Phenotypic plasticity as a long-term memory easing readaptations to ancestral environments – Science Advances

Life & Non-Humans

Abstract

Phenotypic plasticity refers to environment-induced phenotypic changes without mutation and is present in all organisms. The role of phenotypic plasticity in organismal adaptations to novel environments has attracted much attention, but its role in readaptations to ancestral environments is understudied. To address this question, we use the reciprocal transplant approach to investigate the multitissue transcriptomes of chickens adapted to the Tibetan Plateau and adjacent lowland. While many genetic transcriptomic changes had occurred in the forward adaptation to the highland, plastic changes largely transform the transcriptomes to the preferred state when Tibetan chickens are brought back to the lowland. The same trend holds for egg hatchability, a key component of the chicken fitness. These findings, along with highly similar patterns in comparable experiments of guppies and Escherichia coli, demonstrate that organisms generally “remember” their ancestral environments via phenotypic plasticity and reveal a mechanism by which past experience affects future evolution.

RESULTS

Lowland and Tibetan chickens

Studies of animal and plant domestication have been important in evolutionary biology because of the availability of historical records of evolution and many drastically selected phenotypic traits (1316). Chickens were domesticated from the red jungle fowl in South and Southeast Asia at least 4000 to 4500 years ago (15) and were brought to the Tibetan Plateau by ~1200 years ago (17). Tibetan chickens are now well adapted to their high-altitude environment. For instance, they have more red blood cells (18), a smaller mean blood cell volume (18), and a higher efficiency of oxygen dissociation from hemoglobin in embryos (19) when compared with lowland chickens.

Changes in gene expression and regulation play important roles in environmental adaptations (20), including high-altitude adaptations of humans and other animals (21). To study transcriptomic changes in chicken’s adaptation from the lowland to the highland, we performed mRNA sequencing (RNA-seq) of five tissues (brain, heart, liver, lung, and muscle) from lowland chickens raised in Ya’an (670 m above sea level, Fig. 1A; square 1 in Fig. 1B; “lowland” hereinafter) and Tibetan chickens raised in A’ba (3300 m, Fig. 1A; circle 3 in Fig. 1B; “highland” hereinafter), localities on the opposite sides of the east edge of the Tibet Plateau in Sichuan Province, China. In addition, we conducted reciprocal transplant experiments to investigate whether any detected gene expression difference is a plastic response to a different environment or has genetic basis. Specifically, we hatched and raised lowland chickens in the highland (square 2 in Fig. 1B) and Tibetan chickens in the lowland (circle 4 in Fig. 1B), followed by transcriptome profiling of the same five tissues. In the forward adaptation of lowland chickens to the highland, the three samples marked 1, 2, and 3 in Fig. 1B respectively represent the original adapted stage in the lowland (O, Fig. 1C), the plastic stage in the highland (P, Fig. 1C), and the adapted stage in the highland (A, Fig. 1C). Phenotypic differences between stages O and P are plastic, whereas those between P and A are genetic. When studying the potential reverse adaptation of Tibetan chickens to the ancestral lowland environment, we can view the samples marked 3, 4, and 1 in Fig. 1B as representatives of the original stage in the highland (O, Fig. 1D), plastic stage in the lowland (P, Fig. 1D), and adapted stage in the lowland (A, Fig. 1D). Hence, our experiment allows comparing the relative roles of genetic and plastic phenotypic changes in the adaptation to a new environment and potential readaptation to an ancestral environment. Each of the four samples included eight to nine chickens for each tissue except for the muscle, which had five to six chickens (table S1).

Fig. 1 Reciprocal transplant experiments reveal plastic and genetic differences in gene expression levels between Tibetan and lowland chickens.

(A) Map showing the locations and altitudes of the two experimental sites in Sichuan Province, China. Adapted with permission from (33). (B) Design of the reciprocal transplant experiment, in which each of the two breeds are phenotyped in its native environment and the other breed’s native environment. Different symbols indicate different breeds, and different colors indicate different environments where chickens are hatched, reared, and phenotyped. The four samples of chickens phenotyped are numbered. The diagram shows the scenario in which the gene expression level is higher in sample 3 than in sample 1, but the same principle applies if the opposite is true. (C) Samples 1, 2, and 3 in (B) respectively represent the original (O), plastic (P), and adapted (A) stages during the forward adaptation from the lowland to highland. (D) Samples 3, 4, and 1 in (B) respectively represent the O, P, and A stages during the reverse adaptation from the highland to the lowland. (E) Venn diagram showing the number (and percentage) of differentially expressed genes (DEGs) in each tissue between samples 1 and 3 in (B). (F) Numbers of DEGs in each tissue that are plastic-change-only (PO) or genetic-change-needed (GN) during the forward (F) or reverse (R) adaptation. (G) Number of DEGs that are PO divided by the number of DEGs that are GN in F or R adaptation in each tissue and all tissues combined. P values are based on a G test of independence.

After mapping RNA-seq reads (data files S1 and S2), we first identified differentially expressed genes (DEGs) between samples 1 and 3 in Fig. 1B [false discovery rate (FDR) < 5%; see Materials and Methods], which inform us about gene expression differences between lowland and Tibetan chickens in their respective native environments. We found 0.3, 0.9, 1.5, 4.6, and 3.1% of annotated genes as DEGs in the brain, heart, liver, lung, and muscle, respectively (Fig. 1E). We found that the DEGs are most significantly enriched with gene ontology terms of (i) metabolic processes of organonitrogen compound, phosphate-containing compound, amino acids, and protein; (ii) regulation of hydrolase [especially guanosine triphosphatase (GTPase)] activity; (iii) regulation of protein phosphorylation; (iv) regulation of protein secretion; (v) actin cytoskeleton organization and phagocytosis; (vi) stress response, intracellular signal transduction, and apoptosis; (vii) immune processes, lymphocyte differentiation, and leukocyte cell-cell adhesion; and (viii) negative regulation of p38 mitogen-activated protein kinases cascade (data file S3).

For each DEG, we evaluated the relative importance of plastic and genetic changes in its expression difference by the following steps. First, we considered the forward adaptation (Fig. 1C) and defined the expression levels at stages O, P, and A as EO, EP, and EA, respectively. Second, if EO and EP differ significantly but EP and EA are not significantly different, then the difference between EO and EA can be entirely explained by a plastic change so is considered “plastic-change-only” (PO). Third, if EP and EA differ significantly, then the expression variation of the gene is “genetic-change-needed” (GN). Note that these two categories are mutually exclusive, but they together do not include all DEGs because of the situation where neither EO nor EA differs significantly from EP despite a significant difference between EO and EA. Fourth, we similarly classified DEGs to the categories of PO and GN in the reverse adaptation to the lowland (Fig. 1D).

Compared with the forward adaptation, the reverse adaptation shows fewer GN genes but more PO genes in four of the five tissues (Fig. 1F). The ratio (RPO/GN) of the number of PO genes to that of GN genes in the reverse adaptation is 2.9 to 106 times the ratio in the forward adaptation (P < 10−8, G test of independence; Fig. 1G) in these four cases. The exception is the brain, which has the fewest DEGs among the five tissues. For the brain, RPO/GN in the reverse adaptation is 1.06 times that in the forward adaptation, which is not significantly different from 1 (P = 0.9; Fig. 1G). When all five tissues are considered together, RPO/GN in the reverse adaptation is 5.3 times that in the forward adaptation (P < 10−129; Fig. 1G). The elevation in RPO/GN during the reverse adaptation is general rather than specific to certain gene ontology terms (fig. S1).

Because at least some gene expression changes are expected to affect fitness, the overall preponderance of GN expression changes in the forward adaptation and that of PO expression changes in the reverse adaptation prompted us to investigate whether similar trends exist for life-history traits that are crucial to fitness. We found that Tibetan chicken eggs have a drastically higher hatchability than lowland chicken eggs when both are incubated in the highland (Fig. 2A), indicating GN in hatchability during the forward adaptation to the highland. When incubated in the lowland, however, Tibetan and lowland chicken eggs show no significant difference in hatchability (Fig. 2A). This observation, together with a significant plastic change in hatchability of Tibetan chicken eggs between the two localities (Fig. 2A), indicates PO in hatchability during the potential reverse adaptation to the lowland. Hence, the variation of egg hatchability also follows the trends observed from the transcriptomic data. We further respectively examined the survival rate of chicken embryos in the first (Fig. 2B), second (Fig. 2C), and third (Fig. 2D) week of hatching, representing temporal components of egg hatchability. Overall, similar trends were observed between these traits and the final hatchability, with minor variations, apparently due to the relatively large standard errors of the survival rate estimates that render some statistical tests underpowered (Fig. 2).

Fig. 2 Reciprocal transplant experiments reveal plastic and genetic differences in egg hatchability between Tibetan and lowland chickens.

(A) Fraction of fertilized eggs that are hatched. (B) Survival rate of chicken embryos in the first week of hatching, measured by the number of viable embryos on day 7 divided by the corresponding number on day 0. (C) Survival rate of chicken embryos in the second week of hatching, measured by the number of viable embryos on day 14 divided by the corresponding number on day 7. (D) Survival rate of chicken embryos in the third week of hatching, measured by the number of viable embryos on day 21 divided by the corresponding number on day 14. Breeds are indicated by symbols, while environments are indicated by colors. Error bars show one standard error based on the binomial distribution. P values determined by a G test of independence are shown as follows: ns, P > 0.06; o, 0.05 < P < 0.06; *, 0.01 < P < 0.05; ** 0.001 < P < 0.01; and ***, P < 0.001.

Other organisms

Notably, studies of potential reverse adaptations of sparrows and deer mice from high-altitude environments to their ancestral low-altitude habitats also suggested a prominent role of plastic transcriptomic changes (22, 23). However, because these earlier studies did not transplant organisms from low to high elevations, they could not compare the roles of plasticity between forward and reverse adaptations. To examine the generality of our findings, especially outside altitude adaptations and vertebrates, we analyzed recently released transcriptomic data from a large study of Escherichia coli experimental evolution (24). Specifically, E. coli cells were first adapted to a benign environment (M9 medium) for ~90 generations, followed by 350 to 800 generations of adaptation to one of 11 harsh environments (M9 medium supplemented with 1 of 11 stressors). Cells adapted to each harsh environment were then single-cloned. Transcriptomic data were collected for the benign environment–adapted strain in the benign environment and each harsh environment, as well as each stress-adapted strain in the benign environment and its respective harsh environment. These data allowed investigating plastic and genetic expression changes in 11 forward adaptations to different stresses and 11 potential reverse adaptations to the benign environment, in the same manner as in Fig. 1B. We found that in all 11 pairs of forward and reverse adaptations, more genes exhibit GN changes in the forward than reverse adaptations (P < 5 × 10−4, one-tailed binomial test; Fig. 3A), and more genes exhibit PO changes in the reverse than forward adaptations (P < 5 × 10−4; Fig. 3B). Combined across the 11 experiments, RPO/GN is 0.63 in forward adaptations but 1.77 in reverse adaptations (Fig. 3C).

Fig. 3 Reciprocal transplant experiments reveal plastic and genetic differences in gene expression levels between E. coli strains adapted to a benign environment and those adapted to 1 of 11 harsh environments.

The experimental design is analogous to that in Fig. 1B. (A) Numbers of DEGs that are GN in the forward adaptation from the benign environment to one of the harsh environments and the reverse adaptation from each harsh environment to the benign environment. Each line represents the result of one pair of forward and reverse adaptations. (B) Number of DEGs that are PO in F and R adaptations. (C) Number of DEGs that are PO divided by the number of DEGs that are GN in F and R adaptations for each harsh environment and all 11 environments combined. P values are based on a G test of independence. The benign environment is the M9 medium with glucose (5 g/liter), whereas the harsh environments numbered 1 to 11 respectively contain cobalt chloride (16 μM), sodium carbonate (32.5 mM), methylglyoxal (350 μM), cetylpyridinium chloride (4.8 μM), crotonate (50 mM), n-butanol (1.25%), methacrylate (8.75 mM), potassium chloride (210 mM), l-lactate (40 mM), sodium chloride (400 mM), and l-malate (30 mM).

In the chicken and E. coli cases, the new environments are harsher than the ancestral environments. To examine whether this relationship is necessary for our observations, we took advantage of a study of one natural and two human-introduced populations of Trinidadian guppies (Poecilia reticulata) that independently adapted from high- to low-predation streams (10). Here, the new environment is not harsher than the ancestral environment. In each case, we analyzed the published brain transcriptomes of the four groups of guppies from reciprocal transplant experiments (10), as in Fig. 1B. We found that RPO/GN in the reverse adaptation is larger than that in the forward adaptation in all three cases, including two cases where the difference is statistically significant (Table 1). Thus, our finding holds regardless of the relative stress of the ancestral and new environments.

Table 1 Numbers of DEGs that are PO or GN in adaptations of guppies from high- to low-predation streams and the potential reverse adaptations.

Independent evolution occurred in a natural population (natural) and two human-introduced populations (Intro1 and Intro2).

View this table:

DISCUSSION

In this work, we observed a tendency that, compared with genetic phenotypic changes, plastic phenotypic changes play a more prominent role when organisms return to ancestral environments than when they adapt to new environments. This trend appears quite general, because the pattern we initially observed in chickens also exists in guppies and E. coli. The same trend holds in all types of adaptations analyzed, varying from an environmental change in predation to elevation, nutrient, and toxin, regardless of whether the ancestral environment is less or more stressful than the new environment. The observation holds not only for gene expression levels but for life-history traits such as egg hatchability. Why does returning to ancestral environments yield a higher RPO/GN than adaptations to new environments? Below, we evaluate several potential explanations.

First, although a plastic change occurs without mutation, an organism’s ability to make a plastic change has genetic basis. Hence, natural selection can drive the spread of genotypes capable of making certain plastic changes that are beneficial. Consider, for example, the scenario that chickens are frequently exposed to both lowland and highland environments for many generations or the scenario that a strong bidirectional gene flow exists between the lowland and highland populations. It is likely that gene regulatory programs will emerge that turn the chicken transcriptome to the highland-preferred state when in the highland and to the lowland-preferred state when in the lowland, provided that mutations creating such regulations exist and do not have harmful pleiotropic effects (25, 26). This model, however, predicts similar numbers of PO changes in the two environments, so it cannot explain the asymmetry in the relative prevalence of PO changes observed. Furthermore, it is clear, at least in the E. coli and guppy experimental evolution, that no repeated environmental swaps or gene flows occurred.

Second, when a population faces a new environment for a relatively short time, ancestral alleles that are fitter in the ancestral environment may not have been completely replaced by new alleles that are fitter in the new environment. Consequently, the reverse adaptation of this population to the ancestral environment can use standing genetic variants and need not wait for new beneficial mutations (27, 28). But this model is inapplicable in the present study because of the difference in experimental design. For example, regardless of the underlying genetic architecture, lowland-preferred genotypes should be rare among Tibetan chickens such that it is extremely improbable for all of the Tibetan chickens hatched, raised, and transcriptomically profiled in the lowland to be of such genotypes. The same can be said for the guppies and E. coli.

Because the above two simple models cannot explain our finding, we now consider the possibility that our observation reflects a property of the genetic mechanism by which adaptation to a new environment occurs. For example, let us assume that lowland chickens express a particular gene at a basal level when in the lowland (Fig. 4A) and that the preferred expression level of the gene is higher in the highland. However, the preferred higher expression cannot be realized by plasticity alone (Fig. 4B), because the realization requires the evolutionary acquisition of a regulatory motif that binds to a highland-specific transcriptional activator (Fig. 4C). However, when highland-adapted chickens are moved back to the lowland, the highland-specific transcriptional activator is turned off, causing the expression of the gene to return to the basal level without the need for any genetic change (Fig. 4D). In other words, the reaction norm of highland-adapted chickens becomes different from that of the lowland-adapted chickens, as a byproduct of the adaptation to the highland. The same model involving a transcriptional repressor instead of activator applies if the preferred expression level is lower in the highland than in the lowland. Thus, under this model, a forward adaptation to a new environment requires genetic changes that are de novo or standing in the population, whereas returning to the ancestral environment needs only plastic changes, the asymmetry being caused by the establishment of a basal expression program before the emergence of a condition-specific transcriptional regulation in the new environment. More broadly speaking, our observations are explainable by the presence of a sizable fraction of mutations acquired in new environments that have phenotypic effects in the new but not ancestral environments. We propose that this model reflects a general genetic logic of environmental adaptations, and future studies should test this proposition.

Fig. 4 A potential mechanism responsible for long-term memories of ancestral environments in the form of phenotypic plasticity by gene expression regulation.

(A) Lowland expression of a hypothetical gene in lowland-adapted chickens. (B) Highland expression of the gene in lowland-adapted chickens. (C) Highland expression of the gene in highland-adapted chickens. (D) Lowland expression of the gene in highland-adapted chickens. The gray rectangle shows a regulatory motif that binds to a basal transcription factor (gray oval) to allow transcription at the basal level. Green diamonds indicate the mRNA molecules produced. The orange rectangle shows a regulatory motif that binds to a highland-specific transcription factor (orange oval) to increase the transcription level. From (A) to (B), a plastic change is unable to raise the expression level because of the lack of the orange motif, whereas from (C) to (D), a plastic change can lower the gene expression level because of the turnoff of the highland-specific transcription factor in the lowland. A similar model in which the highland-specific transcription factor is a repressor can explain should the optimal gene expression be lower in the highland than in the lowland.

Our study has several caveats worth discussing. First, in interpreting the transplant experiments, we ignored the possibility of maternal effects and environmental effects before the transplant. If these factors are the causes of the asymmetric observations from the reciprocal transplants, then these factors must be different between the lowland and Tibetan chickens (and the relevant populations of guppies and E. coli), which is highly unlikely. Second, while the forward adaptations actually took place in the cases analyzed, the reverse adaptations are inferred under the assumption that the optimal phenotype in an environment is the same for conspecifics of different genetic backgrounds. Although this assumption may not hold for every trait, previous studies suggested that it is true for expression traits relevant to the environment (26, 29). Furthermore, plasticity renders Tibetan chickens as adapted to the lowland as are lowland chickens at least in terms of egg hatchability, suggesting that plasticity eases readaptations to ancestral environments. Third, not every gene expression change is adaptive, nor does our analysis assume so. Hence, a higher RPO/GN in reverse than forward adaptations means that reverse adaptations are accompanied by relatively more plastic expression changes than forward adaptations. Because at least some expression changes are likely adaptive, it is plausible that plastic expression changes contribute more to reverse than forward adaptations. Note that, although there is no one-to-one correspondence between genetic phenotypic changes and the underlying genetic changes, everything else being equal, our results do suggest that the number of underlying genetic changes for the observed phenotypic changes is greater in the forward than reverse adaptations. Fourth, while it is clear that Tibetan chickens originated from lowland chickens, the specific lowland breed analyzed here may not represent the lowland sister lineage of Tibetan chickens. That is, other than a relatively low altitude, Ya’an may not represent the ancestral environment of Tibetan chickens. This potential environmental mismatch should make our conclusion about chickens conservative, because potentially only the fraction of genes related to the altitude adaptation show an increase in RPO/GN in the reverse adaptation compared with the forward adaptation. Potential environmental mismatches do not exist for the bacteria and guppies analyzed.

Our findings raise two critical questions. First, do the patterns observed depend on how long ago the population left the ancestral environment? The human-introduced guppies had been in the new environment for only three to four generations, while the E. coli experiment had hundreds of generations in the new environments and the Tibetan chickens have been on the highland for at least 2400 generations. Under the model proposed, the plastic change in the ancestral environment (Fig. 4D) is expected regardless of the time since the population left the ancestral environment, as long as the adaptation to the new environment has been established. Second, how many ancestral environments will organisms “remember”? Do they remember only the immediate past or multiple previous environments? Clearly, more studies are needed to answer these questions.

In summary, our work uncovers a phenomenon conserved from bacteria to vertebrates that organisms remember their ancestral environments in the form of phenotypic plasticity. Although this observation does not mean that no genetic adaptation is needed when organisms return to their ancestral environment, it does suggest that fewer genetic phenotypic changes are required for readaptations to ancestral environments than adaptations to novel environments. Our finding contributes to the recent debate on the relative roles of plastic and genetic changes in adaptation and reveals the importance of considering whether the environment is changing to a novel or ancestral one. Our results demonstrate the impact of evolutionary histories on future evolutionary dynamics and are therefore a manifestation of historical contingency. It remains to be seen whether this phenomenon can be used to infer ancestral environments of organisms.

MATERIALS AND METHODS

Chicken tissue sample preparation and RNA-seq

We hatched and reared Tibetan chickens and lowland chickens (Pengxian yellow chicken breed) in Ya’an, as well as Tibetan chickens and lowland chickens in A’ba. For each of the four groups of chickens, we collected five tissues (brain, liver, lung, heart, and muscle) from five to nine healthy males of similar body weights at the age of 120 days (table S1). Total RNA was isolated from these tissues by the standard TRIzol method (Invitrogen). RNA sequencing was performed on an Illumina HiSeq 2000 or NovaSeq 6000 sequencing system, and a total of 1384.50 Gb of sequences from 165 libraries were generated. Read mapping was performed using TopHat2 (30) against the reference genome Galgal5.0 from Ensembl release 92. The raw reads were submitted to the National Center for Biotechnology Information (NCBI) and National Genomics Data Center (NGDC) (https://bigd.big.ac.cn/) with the reference number of GSE85318 and PRJCA001948, respectively.

Analysis of chicken transcriptomes

We used DESeq2 (31) to identify DEGs using the cutoff of an FDR of 0.05. In the forward adaptation, the expression difference of a DEG is considered PO if the expression levels from samples 1 and 2 in Fig. 1B are significantly different but those from samples 2 and 3 are not significantly different based on DESeq2 (FDR cutoff of 0.05 considering tests of all DEGs). In the forward adaptation, the expression difference of a DEG is considered GN if the expression levels from samples 2 and 3 are significantly different. In the reverse adaptation, the expression difference of a DEG is considered PO if the expression levels from samples 3 and 4 are significantly different but those from samples 4 and 1 are not significantly different. In the reverse adaptation, the expression difference of a DEG is considered GN if the expression levels from samples 4 and 1 are significantly different.

Phenotyping of chicken life-history traits

Eggs from Tibetan chickens and lowland chickens were collected daily from day 2 after the first insemination to day 6 after the last insemination. In total, we collected 466 eggs from Tibetan chickens and 534 eggs from lowland chickens in A’ba. We also collected 576 eggs from Tibetan chickens and 512 eggs from lowland chickens in Ya’an. All eggs were incubated at 37.8°C with a relative humidity of 60% in an automatic incubator (Changsha Photosynthetic Solar Energy Co. Ltd.). Embryonic mortality was determined using an egg candler. This animal handling experiment was approved by the Institutional Animal Care and Use Committee of Sichuan Agricultural University under permit number YCS-B20100804.

When comparing the hatching performance of different samples (Fig. 2A), we performed a G test of independence using a 2 × 2 contingency table, where the columns are different samples and the rows are the number of eggs hatched or not hatched. When comparing the hatching performance in the second week (Fig. 2C) and in the third week (Fig. 2D), only embryos surviving the previous week were considered.

Analysis of E. coli transcriptomes

Horinouchi et al. (24) performed experimental evolution of E. coli first in the M9 medium with glucose (5 g/liter) for ~90 generations and then in 11 different harsh environments for 350 to 800 generations. In each harsh environment, five replicate evolution experiments were conducted with the same starting strain. Upon the completion of the experimental evolution, single clones were picked and transcriptomically profiled. One microarray transcriptome was generated for each of the five replicates from each of the 11 harsh environments in its adapted environment (corresponding to sample 3 in Fig. 1B). In addition, three replicate transcriptomes of the ancestral strain in the M9 medium (corresponding to sample 1 in Fig. 1B), one transcriptome from each replicate of each of the 11 adaptations in the M9 medium (sample 4 in Fig. 1B), and one replicate of the ancestral strain in each harsh environment (sample 2 in Fig. 1B) were generated. All transcriptome data were from table S2 in the paper (24).

For each harsh environment, we first tested the equality of gene expression level between sample 1 and sample 3 by a two-sample t test. Genes with FDR < 0.05 in this test were referred to as DEGs. We categorized DEGs to PO and GN on the basis of the same criteria used in the chicken data analysis. Note that before any two-sample t test, we performed an F test to examine the equality of variance. When the nominal P value of the F test is smaller than 0.05, we conducted a two-sample t test assuming unequal variances; otherwise, we conducted a two-sample t test assuming equal variances.

Analysis of guppy transcriptomes

Ghalambor et al. (10) performed experimental evolution of Trinidadian guppies (P. reticulata) previously adapted to the streams with high numbers of cichlid predators [high-predation (HP) environment] in two different cichlid-free streams [low-predation (LP) environment]. After 1 year of experimental evolution (approximately three to four generations), fish were collected from each introduced population (Intro1 and Intro2). In addition, fish originating from an HP environment and naturally adapted to an LP environment were collected. In the laboratory, breeds of Intro1, Intro2, and the natural LP population were respectively reared in an LP environment (corresponding to sample 3 in Fig. 1B) and an HP environment (corresponding to sample 4 in Fig. 1B). Breeds of the ancestral population sampled from the HP stream were also respectively reared in the LP (corresponding to sample 2 in Fig. 1B) and HP environments (corresponding to sample 1 in Fig. 1B). Brain RNAs from these reared fish were extracted, sequenced, and quantified. Numbers of fish individuals meeting quality filters are four for Intro1 and Intro2 reared in each environment, three for natural LP reared in HP environment, two for natural LP reared in LP environment, and five for the ancestral HP reared in each environment. The numbers of mapped reads for 37,493 genes were provided by the authors. DESeq2 was used to identify DEGs at FDR < 0.05. For each of Intro1, Intro2, and natural LP populations, we first tested the equality of gene expression level between sample 1 and sample 3, as in Fig. 1B. We categorized DEGs into PO and GN on the basis of the same criteria used in the chicken data analysis.

Software for statistical analysis

The t test and F test were respectively performed using the function “ttest2” and “vartest2” in MATLAB R2019b. FDR adjustment was performed by the function “mafdr” in MATLAB R2019b or the function “p.adjust” with the “BH” method in R. The G test of independence was performed using the function “G.test” from the package “RVAideMemoire” in R.

The enrichment test of gene ontology terms was performed by the web server g:Profiler (32) using the default setting. Specifically, the tool g:GOSt applies a hypergeometric test for each gene ontology term with its default multiple testing correction tool g:SCS. The test was performed in November 2019, when g:Profiler had been updated to Ensembl 98. Only the genes with at least one annotation were considered in the reference background.

SUPPLEMENTARY MATERIALS

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

Acknowledgments: We thank W. Qian, J. -R. Yang, and members of the Zhang laboratory for valuable comments. Funding: This work was supported by U.S. NIH research grant GM103232 to J.Z. and Sichuan Provincial Department of Science and Technology Program 2015JQO023 to D.L. Author contributions: W.-C.H., D.L., Q.Z., and J.Z. designed the research. W.-C.H. and D.L. performed the experiments and analyzed the data. W.-C.H. and J.Z. wrote the paper with input from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The raw sequencing reads were submitted to NCBI (https://ncbi.nlm.nih.gov/) and NGDC (https://bigd.big.ac.cn/) with the reference number of GSE85318 and PRJCA001948, respectively. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.