Abstract

   The origins of multicellular physiology are tied to evolution of gene
   expression. Genes can shift expression as organisms evolve, but how
   ancestral expression influences altered descendant expression is not
   well understood. To examine this, we amalgamate 1,903 RNA-seq datasets
   from 182 research projects, including 6 organs in 21 vertebrate
   species. Quality control eliminates project-specific biases, and
   expression shifts are reconstructed using gene-family-wise phylogenetic
   Ornstein–Uhlenbeck models. Expression shifts following gene duplication
   result in more drastic changes in expression properties than shifts
   without gene duplication. The expression properties are tightly coupled
   with protein evolutionary rate, depending on whether and how gene
   duplication occurred. Fluxes in expression patterns among organs are
   nonrandom, forming modular connections that are reshaped by gene
   duplication. Thus, if expression shifts, ancestral expression in some
   organs induces a strong propensity for expression in particular organs
   in descendants. Regardless of whether the shifts are adaptive or not,
   this supports a major role for what might be termed preadaptive
   pathways of gene expression evolution.

   Subject terms: Gene regulatory networks, Molecular evolution,
   Evolutionary biology, Gene expression
     __________________________________________________________________

   Multicellularity requires complex coordinated gene expression.
   Fukushima and Pollock find that gene expression in different organs is
   likely to constrain future patterns of gene expression evolution,
   particularly following gene duplication.

Introduction

   Vertebrate organs organize physiological activities, and the diverse
   expression patterns of thousands of genes determines organ identities
   and functions. Because of this, the evolution of gene expression
   patterns plays a central role in organismal evolution. The degree of
   organ expression specificity correlates to how fast amino acids
   substitute^[27]1, how rapidly they change expression levels^[28]2, and
   patterns of histone modifications^[29]3. Major organ-altering
   evolutionary events such as development of the hominoid brain are also
   associated with gene expression shifts^[30]4–[31]7. However, although
   gene duplication is well-known to play an important role in expression
   pattern shifts (see e.g., the ortholog conjecture^[32]8–[33]11), the
   evolutionary dynamics of expression patterns with and without gene
   duplication remain poorly understood. An important question is whether
   long-term expression in one organ predisposes genes to be subsequently
   utilized in other organs.

   A possible theoretical basis for such predisposition is the idea that
   certain preexisting adapted states are more conducive to evolution of
   specific new traits than other preexisting states. This is known as
   preadaptation, and when a trait makes such a shift it is referred to as
   exaptation^[34]12. Evidence for preadaptation was long ago found in
   phenotypic traits^[35]13, and recently in molecular traits such as
   protein sequences during de novo gene birth^[36]14 or during functional
   innovations^[37]15. Protein sequence evolution generally involves
   highly epistatic interactions and context-dependent
   changes^[38]15,[39]16 that affect preadaptation, but the modular nature
   of expression regulation^[40]17 makes it unclear whether preexisting
   expression patterns constrain evolutionary outcomes.

   Evolution of gene expression has been studied at genome-wide scales
   mainly using two distinct approaches: phylogenetic and pairwise
   analyses. Phylogenetic approaches model gene expression dynamics and
   infer ancestral expression patterns in the context of gene phylogenies.
   For example, Brownian motion models embody purely neutral expression
   evolution^[41]18, whereas Ornstein–Uhlenbeck (OU) models are designed
   to detect purifying selection and adaptive evolution along with neutral
   fluctuation^[42]19–[43]21. Although each gene family has a distinct
   evolutionary history, a species phylogeny is often used for the sake of
   simplicity. Because such approximations cannot be applied to gene
   families with lineage-specific gene duplications and losses, its
   application has mostly been limited to single-copy genes. In contrast,
   pairwise analysis compares gene expression between paralogs in single
   species^[44]22,[45]23 or between orthologs or paralogs in pairs of
   species^[46]9,[47]11,[48]24–[49]26. Although pairwise approaches can
   evaluate the effect of gene duplications, ancestral expression cannot
   be inferred.

   To infer the adaptive evolution of gene expression in diverse gene
   families, we apply OU models for complex gene family phylogenies
   containing gene duplications and losses, without assuming species
   phylogeny. We also develop a curation pipeline to amalgamate large
   amounts of transcriptome data from many studies for a better
   phylogenetic resolution.

   The results of this study, using these methods and genome-scale
   datasets, show how gene duplication affects evolution of expression. As
   genes evolve, their patterns of expression occasionally shift from
   primarily one organ to another, forming modular connections. Our main
   conclusion is that these shifts are not random. When a shift occurs,
   the organ of primary expression for the ancestral gene strongly
   predicts the organ of primary expression for the descendant gene. We
   conclude that this supports a major role for what can be described as
   preadaptive pathways of gene expression evolution, by which we mean
   that adaptation of a gene for expression and presumably functional
   utility in one organ predisposes it to be more readily utilized for
   primary expression in another organ. A further result of this study is
   that expression shifts are larger and more frequent following gene
   duplication than in its absence. Each shift in gene expression may or
   may not be itself adaptive, but especially after gene duplication they
   are often accompanied by accelerated or decelerated rates of protein
   evolution. We conclude that this and the larger expression shifts
   observed following gene duplication support the idea that gene
   duplication tends to free genes up for regulatory or structural
   functional divergence, and sometimes both.

Results

Duplication-permissive genome-wide analysis of gene families

   To allow evolutionary expression analysis on a broad set of genes, we
   used a phylogenetic approach that deals with the complex history of
   gene family trees with duplications and losses, and applied it to 21
   tetrapod genomes (Supplementary Fig. [50]1). A major challenge in using
   gene trees was divergence time estimation, a prerequisite for applying
   phylogenetic comparative methods. We overcame this problem by
   incorporating phylogeny reconciliation in estimating divergence time of
   gene trees. Gene divergence nodes were constrained by the corresponding
   divergence times in a known species tree, and duplication nodes were
   constrained by ancestral and descendant speciation events (see
   “Methods” for details). Because we estimated individual gene
   phylogenies rather than using a single species phylogeny, we could
   analyze gene families that included many lineage-specific gene
   duplications and losses, making our study less biased toward conserved
   genes with slow gene turnovers^[51]24. Use of gene family trees also
   allowed us to include many organ-specific genes that are enriched in
   lineage-specific and young duplications^[52]24. There were only 1377
   single-copy orthologs for which the species phylogeny was applicable,
   but we were able to include 15,475 genes per species on average
   (including 20,873 human genes, merged into 15,280 gene families). This
   approach eliminates problems with pairwise analyses that ignore
   phylogenetic tree structure, and allowed us to infer expression at
   ancestral nodes in the tree.

Transcriptome amalgamation

   To attain high resolution in our analyses, we amalgamated 1903 RNA-seq
   experiments from 182 research projects (i.e., 182 BioProject IDs in the
   NCBI SRA database) and generated a dataset covering six organs from 21
   vertebrate species without missing data (Supplementary Data [53]1 and
   [54]2 and Supplementary Fig. [55]1). In comparison, other recent
   comparative transcriptomic analyses of
   vertebrates^[56]20,[57]24,[58]26–[59]32 often used the same dataset
   containing 131 RNA-seq experiments from six organs and ten
   species^[60]19, with some additional data in different studies. RNA-seq
   reads were first mapped to corresponding reference genomes and then the
   expression level was quantified by two metrics: transcripts per million
   (TPM) and fragments per kilobase million normalized by trimmed mean of
   M-values^[61]33 (TMM-FPKM). To reduce the among-species variation, the
   TMM normalization was applied across all 1903 samples using the 1377
   single-copy orthologs.

   To allow rapid integrated analysis of datasets, we employed automated
   multi-aspect quality controls, including metadata curation
   (Supplementary Data [62]3), sequence read filtering (Supplementary
   Fig. [63]2), and iterative removal of anomalous RNA-seq samples by
   monitoring correlations between and within data categories (Fig. [64]1a
   and Supplementary Fig. [65]3). The metadata curation enabled us to
   select appropriate samples from the NCBI SRA database. Data that were
   not compatible with those from other research projects (defined by
   BioProject ID) were removed in the correlation analysis by implementing
   a majority rule (Supplementary Data [66]3), resulting in a cleaned
   dataset. This filtering step was designed to fulfill the assumption
   that any samples from the same organ should correlate better than
   samples from different organs within species. When anomalous data were
   detected, all samples belonging to the same research projects (i.e.,
   the same BioProject ID) were discarded.

Fig. 1. Transcriptome amalgamation to integrate heterogeneous RNA-seq
samples.

   [67]Fig. 1
   [68]Open in a new tab

   a A simplified flow chart of the transcriptome amalgamation. The full
   chart is available in Supplementary Fig. [69]1. b–d Transcriptome
   curation within species. Data from Monodelphis domestica with
   SVA-log-TMM-FPKM metrics are shown as an example. The heatmaps show
   Pearson’s correlation coefficients among RNA-seq samples (b). Each row
   and column corresponds to one RNA-seq sample. The expression levels of
   all genes were used to calculate the correlation coefficients. Note
   that anomalous samples contaminated in the curated metadata (low
   correlation samples in 1) are successfully removed, and that
   project-specific correlations visible in the uncorrected data (marked
   2) are absent in the corrected data (marked 3). The boxplots show
   distinct distributions of the correlation coefficients depending on
   whether a pair of samples are the same organ or whether they are from
   the same research project (c). The numbers of comparisons are provided
   in the plot. The correlation coefficients are largely improved in
   within-organ comparisons when surrogate variables are removed, while
   within-project biases are attenuated. In this species, nine surrogate
   variables were detected against 52 RNA-seq data from eight projects
   (d). Analysis of those variables by linear regression highlights the
   BioProject feature as the strongest source of removed biases. For full
   description of predictors, see Supplementary Fig. [70]4. e A principal
   component analysis using expression levels of 1377 single-copy
   orthologs from 21 species. Points correspond to RNA-seq samples. Curves
   show the estimated kernel density. Explained variations in percentages
   are indicated in each axis. f Estimated organ-wise expression levels of
   a housekeeping gene. Since data from relatively many BioProjects are
   available, glyceraldehyde-3-phosphate dehydrogenase gene (GAPDH,
   Ensembl gene ID: ENSGALG00000014442) in Gallus gallus is shown as an
   example. Points correspond to the average expression level calculated
   by random subsampling. All data points and the median value (bar),
   rather than a boxplot, are shown if the number of subsampled BioProject
   combinations is <10. Boxplot elements are defined as follows: center
   line, median; box limits, upper and lower quartiles; whiskers, 1.5 ×
   interquartile range; points, outliers.

   Finally, we applied surrogate variable analysis (SVA)^[71]34 to detect
   and correct hidden biases likely originating from heterogeneous
   sampling conditions and sequencing procedures among experiments in both
   log[2]-scaled TPM and TMM-FPKM (SVA-log-TPM and SVA-log-TMM-FPKM,
   respectively; Supplementary Fig. [72]4a, b). This correction greatly
   improved the correlation of expression levels within organs from the
   same species, even when data were derived from different research
   projects (Fig. [73]1b, c and Supplementary Fig. [74]4c–f). Among
   surrogate variables, BioProject IDs tend to show a high predictive
   power, suggesting project-specific sources of bias (Fig. [75]1d and
   Supplementary Fig. [76]4g–h). Although the inclusion of many species
   from phylogenetically diverse lineages makes it difficult to extract
   organ-wise characteristics from the limited number of single-copy
   orthologs, a principal component analysis produced moderate organ-wise
   segregation in the multispecies comparison (Fig. [77]1e and
   Supplementary Fig. [78]4i–k), further indicating that the curated
   dataset is sufficiently reliable for use in cross-species expression
   pattern analyses. The previously-reported uniqueness of testis
   transcriptomes^[79]19 was partly resolved as the third principal
   component (Supplementary Fig. [80]4k).

   To further evaluate the validity of amalgamated transcriptomes, we
   analyzed the expression of community-curated cell-type-specific marker
   genes associated with organs in PanglaoDB^[81]35, which organizes a
   number of single-cell RNA-seq experiments in human and mouse. We
   compared median values of log-transformed expression levels of >100
   marker genes in each organ (Supplementary Fig. [82]5). After SVA
   correction, all RNA-seq samples in the both species showed the
   corresponding marker expression values higher than those from the other
   organs, suggesting our amalgamated transcriptomes preserve the
   organ-specific gene expression. In the cell-type-wise analysis, a few
   cases, such as juxtaglomerular cells in kidney and hepatic stellate
   cells in liver, could not resolve our organ-wise transcriptomes
   (Supplementary Dataset 10.17632/3vcstwdbrn.1). However, such low
   performance was seen in all samples rather than subsets associated with
   particular BioProject IDs, suggesting that the dissection decisions
   have negligible effects to cell-type compositions in the organs.

   In addition to the better phylogenetic coverage, greater accuracy of
   estimated expression levels is another possible advantage of
   integrating many RNA-seq datasets. This idea is supported by
   subsampling analysis on a housekeeping gene glyceraldehyde-3-phosphate
   dehydrogenase, where, as more data are used, estimated expression
   levels in different organs tend to quickly converge to a similar range
   of values (Fig. [83]1f, ca. 15 SVA-log-TMM-FPKM; Supplementary
   Fig. [84]4l, ca. 11.5 SVA-log-TPM).

Modeling expression evolution

   We next used the amalgamated transcriptomes to evaluate how expression
   evolved along 15,280 maximum-likelihood gene family phylogenies
   (Supplementary Data [85]4), employing multi-optima OU models^[86]36 to
   allow for possible adaptive shifts of optimal expression levels and
   neutral fluctuations^[87]19–[88]21. This modeling identified
   statistically supported expression regime shifts^[89]36,[90]37 on each
   gene tree (Fig. [91]2a), which were then analyzed in the context of
   preceding duplication events. Speciation events (S node; Fig. [92]2b)
   with no duplication were considered the baseline mode of expression
   evolution because regulatory environments and expression patterns are
   more preserved among orthologous genes in comparison with paralogous
   genes produced by gene duplication^[93]8–[94]11. Because OU shift
   detection has been applied for gene expression by assuming species tree
   phylogeny in single-copy genes, shifts in S branches are equivalent to
   those characterized previously^[95]19,[96]20 but also include many more
   speciation events in duplication-prone gene families. Gene tree nodes
   associated with preceding duplication events were categorized as
   DNA-based duplication or retrotransposition (D or R nodes,
   respectively) depending on complete intron losses (Fig. [97]2b).

Fig. 2. Expression evolution in a complex history of gene family evolution.

   [98]Fig. 2
   [99]Open in a new tab

   a Modeling expression evolution with multi-optima Ornstein–Uhlenbeck
   process. A phylogenetic simulation is shown. Colors show branches
   belonging to different regimes. Regime shifts (change of color) appear
   as a substantial change in optimal trait values. b Nodes and branches
   of gene family phylogeny were categorized into S, D, or R based on the
   branching events, i.e., speciation, DNA-based duplication, or
   retrotransposition, respectively. c Venn diagrams of expression regime
   shifts. Circles represent the sets of branches where regime shifts were
   detected with SVA-log-TMM-FPKM or SVA-log-TPM. d The gene tree of
   phosphoglycerol kinases (orthogroup ID: OG0002332) is shown as an
   example. This orthogroup shows ortholog-specific expression patterns as
   well as regime shifts after speciation and lineage-specific gene
   duplication. Tips correspond to genes. The colors of branches and tip
   labels indicate expression regimes. Node colors match to the
   categorization in b. The heatmap shows expression levels and
   among-organ expression patterns are shown as a pie chart for each
   regime. To the right, the number of introns and located chromosomes (A,
   autosome; X, X chromosome; Y, Y chromosome) are indicated. For full
   information including complete tip labels and bootstrap branch
   supports, see Supplementary Dataset 10.17632/3vcstwdbrn.1.

   Organ-wise means of the two expression values, SVA-log-TPM and
   SVA-log-TMM-FPKM, were separately used to model expression evolution
   with OU processes. The two analyses resulted in similar numbers and
   characteristics of expression regime shifts (Supplementary
   Fig. [100]6), but the shift locations were sometimes inconsistent. S
   branches were a major source of apparently inconsistent regime shifts,
   whereas branches following duplication (D and R) showed largely
   consistent detection between the two metrics (Fig. [101]2c). Although
   the inclusion of inconsistently detected branches did not change the
   results, we retained only consistently detected regime shifts for all
   downstream analysis to draw a more robust conclusion. While
   SVA-log-TMM-FPKM values were reported in the main text unless otherwise
   mentioned, the comparisons with SVA-log-TPM-based analyses are
   available as Supplementary Information (see below for specific
   citations).

   As an example of this analysis, an orthogroup of phosphoglycerol
   kinases (PGKs), containing all three categories of branching events
   followed by expression shifts (S, D, and R), is shown in Fig. [102]2d.
   This protein catalyzes the first ATP-generating step in the glycolytic
   pathway and is required for most cell types including
   sperms^[103]38,[104]39. PGK1, the original copy on the X chromosome, is
   known to have duplicated independently in eutherians and marsupials to
   produce the autosomal retrocopy PGK2 that compensates the protein
   activity during X-inactivation^[105]40–[106]42. Our automated analysis
   correctly recovered both retrotranspositions as well as the subsequent
   gains of testis-specific expression in eutherian and marsupial
   lineages. This illustrates that our automated genome-wide analysis can
   recover evolutionary trajectories that are compatible with focused
   single gene family analyses (see Supplementary Dataset
   10.17632/3vcstwdbrn.1 for individual gene trees).

Duplication-specific effects in expression evolution

   Across gene trees, per-branch frequencies of expression regime shifts
   were significantly different among S, D, and R branches (P ≈ 0;
   χ^2 = 2.11 × 10^4; χ^2 test). Expression regime shifts were relatively
   rare in S branches, for a probability of 2.2% per branch, and an
   average rate of 2.5 × 10^−4 shifts per MY (million years)
   (Fig. [107]3a, b). In agreement with the idea that gene duplication
   tends to free genes up for functional divergence and enhance long-term
   retention of duplicated copies^[108]22,[109]43, the frequency of regime
   shifts in D branches was four times as much per branch (9.0% per
   branch, at a rate of 1.7 × 10^−3 shifts per MY across all genes and all
   D branches). Thus, although far fewer branches are preceded by
   DNA-based duplication events (65,868 branches) than speciation events
   (542,978 branches), D branches account for over 33% of all regime
   shifts consistently detected by the two expression measures. We note
   that this result reinforces previous results on the ortholog
   conjecture, the idea that duplicated gene copies (paralogs) are more
   prone to expression shifts than orthologs^[110]8–[111]11. R branches
   were far more likely to result in expression regime shifts (37.3% per
   branch, at a rate of 7.0 × 10^−3 shifts per MY), but with only 2963 R
   branches, this resulted in only 1106 shifts (5.6% of the total).
   Translocated genes are more likely to shift expression than those that
   do not (Supplementary Fig. [112]7), in line with previous observations
   from the human genome^[113]23. While the expression shift frequency in
   S and D branches varies slightly across the phylogeny, R branches
   showed much stronger among-lineage heterogeneity, and had a
   particularly high frequency in the mammalian lineage (Fig. [114]3b).
   These retrotransposition-related expression changes may be related to
   the variation in the retrotransposition rate itself, which is known to
   vary across lineages^[115]44,[116]45. Among-species heterogeneity in
   gene prediction quality may also be attributed to this pattern because
   early-diverging species tended to show higher percentages of missing
   single-copy orthologs than those in mammalian species (Supplementary
   Fig. [117]8; but see Danio rerio, Oreochromis niloticus, and
   Oryctolagus cuniculus as counterexamples). In the absence of regime
   shifts, expression levels varied most in ovary and testes, which had
   significantly higher average stationary variances than the other four
   organs on the basis of tree-wise stationary variance (Fig. [118]3c).
   This supports previously observed high variation of gene expression in
   testes^[119]19 and extends it to the reproductive organs of both sexes.

Fig. 3. Characteristics of expression shifts in 15,280 gene trees.

   [120]Fig. 3
   [121]Open in a new tab

   a The species tree showing analyzed genomes and their divergence time.
   A part of animal silhouettes were obtained from PhyloPic
   ([122]http://phylopic.org). The silhouettes of Astyanax mexicanus and
   Oreochromis niloticus are licensed under CC BY-NC-SA 3.0
   ([123]https://creativecommons.org/licenses/by-nc-sa/3.0/) by Milton Tan
   (reproduced with permission), and those of Anolis carolinensis (by
   Sarah Werning), Ornithorhynchus anatinus (by Sarah Werning), and Rattus
   norvegicus (by Rebecca Groom; with modification) are licensed under CC
   BY 3.0 ([124]https://creativecommons.org/licenses/by/3.0/). b Mapping
   of 18,812 expression shifts in the species tree. The number and
   proportion of expression regime shifts in S, D, and R branches are
   shown. Corresponding branches in the species tree are indicated in a. c
   Organ-specific stationary variances (γ) of expression level evolution
   in vertebrates. The distribution of γ between reproductive and
   nonreproductive organs were compared by a two-sided Brunner–Munzel
   test^[125]95. d–f Cumulative frequency of change in organ expression
   specificity (d), change in maximum expression level (e), and expression
   complementarity between sister lineages (f) among detected expression
   shifts. Number of analyzed regime shifts are shown in b. The D
   statistics and P values of pairwise branch category comparisons were
   calculated with two-sided Kolmogorov–Smirnov tests. Boxplot elements
   are defined as follows: center line, median; box limits, upper and
   lower quartiles; whiskers, 1.5 × interquartile range.

   If expression regime shifts are due to functional divergence, it is
   highly relevant to characterize how expression properties changed from
   ancestral to derived regimes. To do this, we examined changes in organ
   expression specificity, maximum expression level, and organ expression
   complementarity. The specificity measure τ ranges from 0 for uniformly
   expressed genes to 1 for genes with highly specific expression^[126]46.
   The distributions of regime shifts in D and R branches are shifted
   toward greater organ specificity compared to shifts in S branches, with
   R branches creating the most specific expression (Fig. [127]3d). To
   characterize the on state transcriptional activity, we analyzed the
   maximum fitted expression levels among the six organs (μ[max]). D and R
   branches appear to be enriched for downregulation compared to S
   branches (Fig. [128]3e). Complementarity of organ expression patterns
   was measured to evaluate the differentiation between a pair of sister
   branches. We used a metric on the fitted organ-wise expression levels
   (μ) called TEC, which ranges from 0 for completely overlapping
   expression to 1 for mutually exclusive patterns^[129]43^. Nearly all
   branches with regime shifts (95%) had complementarity values >0.5,
   indicating that most regime shifts detected involve differentiation of
   expression patterns, rather than overall up- or downregulation. Regime
   shifts in D and R branches often had more complementary expression than
   those in S branches (Fig. [130]3f), further supporting the role of gene
   duplication in functional differentiation. The more drastic effect in R
   branches probably reflects the regular loss of regulatory elements in
   retrotranspositions, whereas DNA-based duplication can more often
   retain regulatory regions^[131]47. Jointly, these results indicate
   that, compared with the baseline from speciation-associated shifts,
   gene duplication tend to produce more organ-specific, more often
   downregulated, and more differentiated expression patterns. Although
   the downregulation may be explainable by a tendency to need less of the
   newly functional expression regime, it may also be explained by either
   recent nonfunctionalization^[132]48,[133]49 or specialized expression
   in organs that were not part of this analysis.

Context-dependent change in the rate of protein evolution

   Change in gene expression can be accompanied by accelerated or
   decelerated protein evolution, which may be detected by change in the
   ratio of nonsynonymous/synonymous substitutions (dN/dS or ω) along
   branches. In D and R branches, median ω values are more than double the
   baseline seen in S branches (Fig. [134]4a; Supplementary Fig. [135]9a),
   again supporting the ortholog conjecture and the idea that gene
   duplication tends to free genes up for functional
   divergence^[136]22,[137]43. Within each of S, D, and R branch
   categories, branches with expression regime shifts accompany an
   increased ω compared to sister branches (Fig. [138]4a). The increased
   rate was quite small in S branches (median ω, 0.096 in branches with
   shifts versus 0.102 in sister branches), much bigger in D branches
   (0.182 versus 0.244), and huge in R branches (0.036 versus 0.394). If
   the changes in expression and the rate of protein evolution are due to
   functional changes, this may indicate that functional divergence is
   sometimes effected by joint changes in expression and accelerated
   protein evolution.

Fig. 4. Context-dependent changes of protein evolution rate coupled with
expression regime shifts.

   [139]Fig. 4
   [140]Open in a new tab

   a Distribution of ω values. A plus (+) indicates branches with
   expression shifts, whereas minus (−) branches are sisters to the plus
   branches. Statistical differences between pairs of distributions were
   tested using a two-sided Brunner–Munzel test^[141]95.
   Non-log-transformed median values are shown above the boxplots. For
   visualization purposes, extreme values exceeding ±10 were clipped.
   Boxplot elements are defined as follows: center line, median; box
   limits, upper and lower quartiles; whiskers, 1.5 × interquartile range.
   b Relationships between protein evolution rate and change in expression
   properties. Points correspond to expression regime shifts. Dashed lines
   indicate no between-branch difference in ω. Solid lines show a linear
   regression. Its slope and number of regime shifts are provided in the
   plot. Regime shifts with negative and positive changes were separately
   analyzed for organ expression specificity (upper) and maximum
   expression level (middle). P values indicate whether the slopes were
   significantly different from zero (two-sided t tests).

   Although 52.3% of all branches with expression regime shifts had higher
   ω relative to sister branches (ω ratio), 27.5% of sister branch pairs
   are relatively undifferentiated, with differences in ω within ±5%, and
   42.5% had lower ω in the branches with expression shifts. This finding
   led us to hypothesize that the direction of the rate changes in protein
   evolution is linked to how, rather than whether, expression is changed.
   Strikingly, we found a context-dependent association between protein
   and expression evolution (Fig. [142]4b and Supplementary Fig. [143]9b).
   Increased ω ratio was linked strongly to increased, rather than
   decreased, organ expression specificity in S and D branches,
   potentially reflecting adaptive evolution coupled with specialized
   expression. However, it was in turn highly associated with decreased
   specificity in R branches, which may be explained by frequent gene
   decay in unsuccessful retro-copies. The change in maximum expression
   level was overall negatively correlated with ω ratio, but this link was
   stronger in downregulation compared with upregulation, except for D
   branches where the ω ratio was smaller when Δμ[max] was larger
   (Fig. [144]4b). It has been reported that high expression slows protein
   evolution^[145]1, and our results suggest that DNA-based duplication
   creates such constraints when accompanied by upregulation. The organ
   expression complementarity between sister lineages was positively
   correlated with ω ratio (Fig. [146]4b), and its association was
   strongest in R branches, suggesting that protein evolution accelerates
   as gene expression patterns differentiate from their ancestral state.
   Collectively, these results suggest that protein evolution rate is
   linked to changes in expression properties through a complicated
   association, which masks their relationships in a global, unstratified
   analysis, and potentially explain a previous report of no strong
   relation^[147]26.

Organ-specific propensity in gene expression evolution

   The preadaptation hypothesis predicts that the ancestral organ
   expression prior to the shifts will affect which organs are likely to
   become the target of newly specific expression. To assess this, we
   tested whether expression shifts are random with respect to change in
   expression from one organ to another, by characterizing the organ in
   which genes are most highly expressed (primary-expressed organ, PEO).

   Across vertebrates, switching from one PEO to another was detected in
   6886, 3586, and 746 regime shifts in S, D, and R branches,
   respectively. The gain/loss ratios are heterogeneous among organs
   (Fig. [148]5a), suggesting that vertebrate organs serve as both sources
   and sinks in expression evolution, but that their relative
   contributions are organ-specific. Although S, D, and R branches shared
   a global trend of relatively abundant testis-related PEO shifts, their
   distributions are largely different (P = 1.52 × 10^−77; χ^2 = 531; χ^2
   test). D branches were moderately similar to both S and R branches
   (Spearman’s ρ ~ 0.6), but S and R branches were dissimilar (ρ = 0.28).
   This pattern, including the abundant shifts related to testis, was
   robust against the correction by the organ-wise numbers of expressed
   genes (Supplementary Fig. [149]10). This result suggests a role for
   gene duplication, including by retrotranspositions, in remodeling the
   among-organ flow of expressed genes.

Fig. 5. Evolutionary dynamics of gene expression.

   [150]Fig. 5
   [151]Open in a new tab

   a Shift distributions of primary-expressed organs (PEOs). Y-axis was
   sorted by abundance in S branches. Spearman’s correlation coefficients
   among S, D, and R branches are shown above the plots. b Preadaptation
   networks in organ expression. Arrows represent transitions from
   ancestral PEOs to derived PEOs, and its color shows statistical
   significance based on 10,000 permutations. The results were obtained
   with SVA-log-TMM-FPKM, and an asterisk (*) indicates the statistical
   significance supported also by SVA-log-TMM-based analysis
   (Supplementary Fig. [152]11). c The global polarity of PEO shifts. The
   global polarity is defined by the scaled sum of differences between two
   opposite PEO shifts. Boxplots show the distribution estimated by 1000
   bootstrap resampling. The D statistics and P values of pairwise branch
   category comparisons were calculated with two-sided Kolmogorov–Smirnov
   tests. Boxplot elements are defined as follows: center line, median;
   box limits, upper and lower quartiles; whiskers, 1.5 × interquartile
   range.

   Controlling the total number of shifts from and to each PEO, some PEO
   shifts are significantly different from the random expectation
   (Fig. [153]5b and Supplementary Data [154]5). There are clear patterns
   of evolutionary transitions that are statistically supported by
   independent OU modeling of the two expression metrics. In S branches,
   the pairs of brain–testis and testis–ovary showed strong connections,
   indicating a solid exchange module. Kidney and liver also donate genes
   to one another, forming a separate module from brain–testis–ovary. D
   and R branches showed a pronounced acceleration of PEO shifts between
   testis and ovary. PEO shifts in S and D branches were moderately
   symmetric in the flow between pairs of organs (Fig. [155]5c), meaning
   two organs tend to donate comparable numbers of expressed genes each
   other. In contrast, R branches were more asymmetric. The lower symmetry
   may have perturbed the evolutionary dynamics of gene expression.

   To check the robustness of our analysis, we analyzed high-confidence
   subsets of expression regime shifts. The most drastic expression
   changes were characterized by introducing a cutoff of organ expression
   specificity (τ > 0.5) to define organ-specific genes. Although some
   previously significant trends were not recovered due to small sample
   sizes, the result was largely consistent with the broader analysis
   (Supplementary Fig. [156]11). Because tree inference errors can bias
   the downstream analysis including OU modeling, we also analyzed
   expression regime shifts found in clades which have a high support in
   tree inference (>99% bootstrap support). Again, the results were
   largely consistent (Supplementary Fig. [157]11). Especially, the
   brain–testis–ovary and kidney–liver modules in S branches and the
   testis–ovary connection in D and R branches were always reproduced in
   the analyses with the above thresholds in combinations with the two
   expression metrics (SVA-log-TMM-FPKM and SVA-log-TPM; Supplementary
   Fig. [158]11). The analysis of shifts in high-support branches also
   reproduced the other main results in this paper (Supplementary Dataset
   10.17632/3vcstwdbrn.1), demonstrating the robustness of our conclusion.

   The effect of gene trees was further examined by replicating the
   analysis with alternative tree topologies inferred by species tree
   reconciliation, which takes into account duplication–loss
   rates^[159]50. This reconciliation step is expected to correct
   erroneous tree topology, while possibly introducing another bias
   derived from over-correction of biological signals such as incomplete
   lineage sorting. With the reconciled trees, the OU modeling with
   SVA-log-TMM-FPKM values resulted in equivalent numbers of expression
   shifts in S and D branches compared with those with non-reconciled
   trees (97% [23,231/23,985] and 104% [9407/9018], respectively)
   (Supplementary Fig. [160]12a). In contrast, the phylogeny
   reconciliation substantially reduced the number of shifts in R branches
   (39% [481/1238]). This could be explained by the correction of
   erroneous tree topology caused by the fast-evolving retro-copies
   (Fig. [161]4a), although the differences in shift numbers did not
   correlate with the topological differences measured by the
   Robinson–Foulds distance^[162]51 (Supplementary Fig. [163]12b).
   Nevertheless, resulting PEO shift distributions were largely similar
   (Supplementary Fig. [164]12c), with the reproduced accelerations in the
   brain–testis–ovary and kidney–liver modules (Supplementary
   Fig. [165]12d) and the asymmetric PEO shifts in R branches
   (Supplementary Fig. [166]12e), suggesting the robustness of detected
   modules against gene tree topology.

   To obtain insight into the biological relevance of the among-organ
   modules, we characterized human genes involved in PEO shifts from all
   branch categories using the Kyoto Encyclopedia of Genes and Genomes
   (KEGG) pathway enrichment analysis. The genes descended from the
   testis–ovary PEO shifts enriched only one KEGG pathway term “cell
   cycle” (Supplementary Data [167]6; adjusted P value < 0.05, Fisher’s
   exact test with the Benjamini–Hochberg correction), likely reflecting
   their function in meiosis. Although the adjusted P value was not
   statistically significant, it is noteworthy that the top-ranked term
   for the brain–testis connection was “endocrine and other
   factor-regulated calcium reabsorption” (unadjusted P
   value = 1.66 × 10^−3; adjusted P value = 0.26) annotated to four genes
   including GNAQ, which has been implicated to tumor formation in
   neuronal tissues^[168]52,[169]53. In the kidney–liver module, 15 terms
   were significantly enriched, many of which appear to be related to the
   functions and diseases in those organs, for example, “bile secretion,”
   “phagosome,” “lysosome,” “ABC transporters,” “sphingolipid metabolism,”
   and “Type-I diabetes mellitus” (Supplementary Data [170]6). These
   results suggest that the among-organ modules in the PEO shifts played a
   role in supplying functionally important genes.

   Because ancestral expression was shown to orient new expression by the
   analysis of PEO shifts, we concluded that there was prevalent
   organ-specific propensity, which supports the presence of preadaptation
   in gene expression.

Discussion

   Our results suggest that the landscape of expression evolution is
   strongly shaped by mechanisms of gene birth. Expression shifts are more
   pronounced following gene duplication in agreement with the results of
   pairwise gene expression analyses^[171]22,[172]24,[173]43, and shifts
   in patterns of PEOs strongly depend on the expression state in the
   ancestral organism. Thus, by analyzing such influences on a genome-wide
   scale for a moderately large number of species, the question whether
   long-term expression in one organ predisposes genes to be subsequently
   utilized in other organs has been answered in the affirmative. There
   are preadaptive propensities in the evolution of vertebrate gene
   expression, and the propensity varies with the presence and type of
   gene duplication. Furthermore, the approach developed in this study,
   using complex gene family phylogenies including gene duplications and
   losses that do not assume perfect match to the species phylogeny, and
   incorporating a curation pipeline to amalgamate large amounts of
   transcriptome data from many studies, was essential to obtain the
   necessary species density and phylogenetic resolution to answer this
   question. The extensibility of this method will allow for more species
   and more organs to be incorporated as further studies come into the
   literature from diverse laboratories.

   The mechanisms responsible for the preadaptive propensities that
   influence expression shifts among organs are, however, unknown. A key
   question in understanding these shifts may be the role of adaptation in
   the shift, and in subsequent evolution. We have been careful so far to
   simply describe the shifts, but adaptive possibilities include
   subfunctionalization, escape from adaptive conflict (EAC), and
   neofunctionalization^[174]54–[175]56. The increased number of shifts
   following duplication suggests that drift alone is not the explanation,
   but subfunctionalization easily could be. Subfunctionalization is the
   idea that, if a gene has multiple functions prior to duplication, they
   may be segregated among the duplicates following gene duplication.
   Thus, the expression shifts may be simply a shift in focus of a
   duplicated copy on a subset of the necessary expression profile needed
   at the organismal level. In this scenario, any accompanying
   acceleration of amino acid substitution would be caused by a loss of
   constraint and reduced purifying selection in one expression
   environment or the other.

   EAC involves more adaptation by adding the simple idea that prior to
   duplication, the multiple functions and expression regimes were at
   least partially in conflict. Such conflicts could clearly occur at the
   amino acid level, but could also occur at the expression level. For
   example, if expression levels were focused on a most-important tissue
   or most-sensitive tissue prior to duplication, but after duplication
   could be more tailored to what is better for the new expression regime.
   Finally, neofunctionalization would occur at the sequence or expression
   level, if the loss of selection on a duplicate allowed mutations that
   were previously harmful to the old function, but now are not, and are
   able to carry out some novel functional aspect that was previously
   prohibited. Neofunctionalization is perhaps the most interesting and
   extraordinary possible cause for the expression regime shifts we see,
   but it requires strong evidence and it is not a necessary explanation
   for what we observe.

   In this context, the patterns of expression regime shifts we observed
   may be explained at different levels of biological organization, from
   the tissues and cells that make up organs, to subcellular compartments,
   chromatin structure, promoter usage, and protein biochemistry. Part of
   the propensity shifts we observed can be explained by the
   out-of-the-testis hypothesis, which posits that accelerated gains of
   testis expression are based on the permissive chromatin state, abundant
   transcriptional machinery, relatively simple promoters required for the
   expression in spermatogenic cells, and following gains of new
   expression patterns^[176]4,[177]57. This theory fits to the accelerated
   testis-related PEO shifts, and could fit with any of the adaptive
   scenarios discussed above, but the other detected patterns (e.g.,
   kidney–liver module) require other explanations.

   One potential mechanism of preferences in expression regime shifts is a