Abstract Phenotypic variation among individuals provides the raw material for evolution, and gene expression is a key mediator between genetic and phenotypic variation. As for phenotypes, the range of gene expression is limited and biased by evolutionary and developmental constraints. Observed expression variability due to biomolecular stochasticity and cell-to-cell heterogeneity has been well-studied in isogenic populations of unicellular organisms. However, for multicellular organisms with a diversity of cells and tissues sharing the same genetic background, the interplay between expression variability, gene and organ function, and gene regulation remains an open question. Here, we use highly multiplexed 3’-end bulk RNA sequencing to generate transcriptome profiles spanning at least nine organs in outbred individuals of three ray-finned fish species: zebrafish, Northern pike, and spotted gar. Per organ, we quantify individual-to-individual gene expression variability independent of mean expression level. Lowly variable genes are enriched in cellular housekeeping functions whereas highly variable genes are enriched in stimulus-response functions. Furthermore, highly variable genes evolve under weaker purifying selection at the protein-coding sequence, indicating that intra-species expression variability predicts inter-species protein sequence divergence. Genes that are broadly expressed across organs are both highly expressed and lowly variable, whereas organ-biased genes are typically highly variable within their top organ. Among organ-biased genes, patterns of expression variance are dependent on the top organ. Specifically, brain- and gonads-biased genes have lowly variable expression across different organs, suggesting stabilizing selection. These patterns suggest that gene regulatory mechanisms evolve under organ-specific selective pressures. Files Bucao_VariabilityAcrossOrgans_Supp.zip contains scripts and results files for all steps of the analysis pipeline, excluding promoter sequence analysis (contained in Supplemental_Code_promoter-sequence-analysis.zip). Directory Structure config/: Contains YAML file indicating package versions for conda environment data/: Contains input data counts/: Contains counts and UMI-deduplicated counts. gene_metadata/: Contains gene biotype information from Ensembl sample_metadata/: Contains sample metadata files for each species selectome/: Contains selection statistics from the Selectome database results/: Contains output files sorted by subfolders labeled after each step of the analysis pipeline. Only R notebook HTML files are available on the Git repository, please check Zenodo for R data files. run_pipeline.Rdata and run_pipeline_revisions.Rdata: Contains all parameters used for each step of the analysis pipeline workflow/: Contains scripts used for the analysis pipeline analysis/: Contains all steps of the analysis pipeline, available as .Rmd files functions/: Contains all functions used for analysis/ renv/: Used for package management in R run_pipeline.R and run_pipeline_revisions.R: Runs all the steps under analysis/ run_go_figure.sh: Runs GO-Figure! 1.0.0 (downloaded separately) demultiplex_brbseq_fastq.sh: Used for demultiplexing BRB-seq fastq files using BRB-seqTools 1.6.1 (downloaded separately) for uploading to NCBI SRA rename_fastq_files.sh: Used for renaming demultiplexed fastq files by mapping each barcode to their corresponding sample name renv.lock: Lockfile for managing R package versions. Run renv::restore() to set up the R environment based on packages specified in the lockfile. All package versions used are also specified in the output HTML files under results/. Species Codes LOC: Lepisosteus oculatus (spotted gar) ELU: Esox lucius (Northern pike) DRE: Danio rerio (zebrafish)