Expression Plasticity of Transposable Elements Is Highly Associated with Organismal Re-adaptation to Ancestral Environments.

Understanding the roles of phenotypic plasticity in adaptive evolution has gained recognition for decades. Studies involving multiple taxa have shown that gene expression plasticity serves as "long-term memory" to facilitate re-adaptations to ancestral environments. Nevertheless, the general pattern and the underlying genetic basis of expression plasticity remain unclear. The transposable elements (TEs) play crucial roles in gene expression regulation and are widely distributed within the genome. Given this, we re-analyzed the transcriptomic data of chicken (Gallus gallus) generated from a reciprocal transplant experiment to examine whether expression shifts of TEs are involved in the re-adaptation process. Similar to the protein-coding genes, the plastic changes of TEs overwhelmingly exceed the genetic changes in the re-adaptation process. Further, the associated TEs co-expressed with diverse genes to perform a regulatory activity. Thus, our study supports the general function of phenotypic plasticity in adaptive evolution, and suggests a regulatory functions of TEs in this process.

Organismal adaptation to a new environment due to the plastic and genetic changes; however, the relative role of the two processes in the adaptation is not clear. A previous study presented that compared to genetic changes, phenotypic plasticity played a great role when organisms return to the ancestral environment. Our study suggests that transposable elements (TEs) may play an important role in this adaptative process. It would improve our understanding of the regulation of TEs and their function in the phenotypic plasticity during the organismal adaptive evolution.

Organismal evolution depends on the accumulation of different patterns of adaptations to the changing environments (Ho et al. 2020). Understanding how organisms adapt to altered environments has become a long-standing issue in evolutionary biology. Generally, organismal adaptations to the altered environments can result from phenotypic plasticity and/or genetic variation (Phillimore et al. 2010; Ho and Zhang 2018). On the one hand, phenotypic plasticity illustrates an organismal response to changing environments by altering its phenotypes without the involvement of gene mutations (Duncanm et al. 2020). On the other hand, genetic variation is associated with a genomic mutation(s) driven mainly by natural selection (Ho and Zhang 2018, 2019). Given the important role of adaptive evolution in an organism’s survival and persistence in the environment, deciphering the relative role of the two phases (plastic vs. genetic changes) has attracted attention in the past decades (Pfennig et al. 2010; Rodriguez-Verdugo et al. 2016; Sun et al. 2020; Zhang et al. 2021).

To date, two hypotheses have been advanced to explain the roles of plasticity in organismal adaptations to novel environments. First, a traditional view based on a few phenotypic traits argued that plasticity serves as a “stepping stone” that allows an organism to rapidly acquire the optimum phenotype before observable genetic and phenotypic changes in the new environment (Levis and Pfennig 2016). On this basis, Jourdan-Pineau et al. (2012) observed that the Parsley frog used phenotypic plasticity rather than genetic variations to adapt to two distinct temporal niches. Nevertheless, studies have shown that phenotypic plasticity varies vastly among species and populations, with the widely distributed species/populations having higher morphological plasticity than others with narrowly distributed niches (Svanback and Schluter 2012). Second, the current hypothesis based on analyzing the expression data of thousands of genes proposed that the plastic expression changes are generally reversed, rather than reinforced, by subsequent genetic expression changes during environmental adaptation; thus, plastic changes are non-adaptive (Ghalambor et al. 2015; Ho and Zhang 2018). Whether or how plasticity precedes and facilitates organismal adaptation to novel environments remains unknown.

A recent study (Ho et al. 2020) based on reciprocal transplant studies of Tibetan chicken from low- and high-elevation showed that the plastic gene expression changes mainly serve as a “long-term memory” that could ease organismal re-adaptations to their ancestral (rather than the novel) environments. This study revealed more genetic changes in the “forward adaptation” from lowland to highland. In contrast, more plastic changes were commonly acquired when Tibetan chickens were brought back to the lowland (reverse adaptation). Although this study highlighted the relationship between plastic and genetic changes as well as proposed the central role of phenotypic plasticity in organismal evolution, the genetic bases, especially whether specific genomic components contribute to the plastic gene expression changes during the reverse adaptations and whether the noncoding elements also display a similar trend, are still unclear.

TEs are DNA sequences that can replicate or move from one locus to others in the genome (Bourque et al. 2018). The TEs occupy a substantial portion of the eukaryotic genome, which are commonly classified into two major classes according to their mechanism of transposition: (1) DNA transposons which can move their locations in the genome through a “cut-and-paste” mechanism (class I) and (2) retrotransposon which can replicate themselves through a “copy-and-paste” way (class II). Although most TEs can lose the capacity of transposition over long-term evolution and could be fixed within the genome, there are yet growing pieces of evidence that the inserted TEs can be transcribed from their locations (Bourque et al. 2018) and can further regulate host gene expression by providing multiple cis-regulatory sequences for the host genome (Sundaram and Wysocka 2020). Hence, the diverse inserted TEs could lead to phenotypic changes and adaptation by modifying the host transcriptional networks (Brawand et al. 2011; Ghalambor et al. 2015; Gilbert et al. 2021). In addition, the landscapes of the TEs could reflect the organism’s ecological niche (Gilbert et al. 2021). Based on these, we ask whether diverse TEs show dynamic expressions and are implicated in the process of reverse and/or forward adaptation.

To address this, we re-analyzed the transcriptomic data of chicken (Gallus gallus) from the reciprocal transplant studies of Ho et al. (2020). This data set contains the gene expression profiles of five tissues (brain, heart, liver, lung, and muscle) of chicken from forward (lowland to highland) and reverse (highland to lowland) adaptation processes. Similar to the study of Ho et al. (2020), we defined three groups of samples as follows: (O) chicken living in the original environment, (P) chicken moved to an altered environment from their original locations, and (A) chicken that has adapted to the new environment. Expression changes between stages O and P are considered plastic, whereas those between P and A are genetic (Ho et al. 2020). For each adaptation process, we mapped the RNA sequencing reads onto the chicken reference genome (version galGal6 available from UCSC Genome Browser) and calculated the expression levels of TEs (as well as the protein-coding genes) at each stage with the SQulRE pipeline (Yang et al. 2019). This pipeline allowed us to calculate and compare the expression values of the same locus between the highland and lowland populations We then determined the plastic-change-only (PO) genes/TEs which showed significantly different expressions between the stages O and P (but without a significant difference between P and A), as well as the genetic-change-needed (GN) genes/TEs whose expressions were significantly different between the stages P and A (fig. 1) (Ho et al. 2020). Given this, we estimated the ratio of PO/GN for both forward and reverse adaptation processes. Interestingly, the ratio of PO/GN of the reverse adaptation was much higher than that of the forward adaptation in almost all tissues except the brain (fig. 1). The same pattern has been observed for the protein-coding genes (Ho et al. 2020). This indicates that the plastic expression changes of TEs might be highly involved in the formation of phenotypic plasticity, especially in the reverse adaptation. This result supports the general trend that phenotypic plasticity plays an essential role in the reverse adaptation to ancestral environments, as reflected in protein-coding genes and non-coding elements (like TEs).

Fig. 1.

The plastic and genetic differences in TEs expression levels and data analysis of DETEs between highland and lowland chickens in forward and reverse adaptation. (A) The diagram of forward and reverse adaptation. O, P, and A represent different stages, respectively. (B) Ratio—the number of PO DETEs (differentially expressed TEs) divided by the number of GN DETEs in forward adaptation and reverse adaptation. (C) The number of different types of DETEs in the GN process. (D) The number of different types of DETEs in the PO process. (E) The number of different LINE types in the GN process. (F) The number of different LINE types in the PO process. Note that snRNA, rRNA, tRNA, and scRNA were classified into the other type of TE due to their extremely small number.

To date, numerous TEs have been identified in vertebrate genomes. Thus, we asked whether some specific types of TEs play important roles in reverse adaptation. To address this, we classified the differently expressed TEs (DETEs) according to the chicken genome annotation and compared the two adaptation processes. As shown in fig. 1, we found that LINE (Long Interspersed Nuclear Element) occupied 61.46–72.47% of the total DETEs in the GN process and 56.2–73.91% in the PO process (fig. 1) of the five tissues. Among the differently expressed LINEs, the family—CR1 (chicken repeat 1)—accounts for 98.71–100% in the GN process (fig. 1) and 98.45–100% in the PO process (fig. 1). The CR1 has been previously reported in the G. gallus (Stumph et al. 1981; Hillier et al. 2004) with a total length of 4.5 Kb and two open reading frames (Bertocchi et al. 2018). The CR1 encompasses the largest number of copies (>200,000) in the chicken genome, accounting for more than 98% of the LINE and 64% of the whole genome (Hillier et al. 2004). Such a high proportion of the content of CR1/LINE may be the reason for the highest number of the CR1-type DETEs; however, the CR1 may still be the most important material resource for the host genome (especially the chicken) to adapt to changing environments, which needs to be further investigated in the future.

Given the regulatory capability of TEs (Cowley and Oakey 2013; Lopez-Maury et al. 2009; Slotkin and Martienssen 2007), we further evaluated if the DETEs are involved in the regulation of gene expression. To address this, we first tested whether the DETEs are located around the differently expressed genes (DEGs). To do this, we selected a threshold of 5,000 bp upstream of the DETEs and counted the numbers of DEGs in these regions. As shown in table 1, the muscle tissue in the forward PO process had the least number of DEGs (n = 154) and DETEs (n = 105). At the same time, the heart and lung tissues in the forward GN process had the largest number of DEGs (n = 1,588) and DETEs (n = 1,426), respectively. The analyses additionally showed that an average of 4.5% (1.9–8.5%) of the DEGs are located within 5,000 bp upstream of the TEs (table 1). We further checked whether these intersections are randomly occurred. For each tissue, we randomly sampled the same number of TEs as the DETEs from the background TEs and then determined their overlapping DEGs in the same way. We repeated this analysis 1,000 times. Finally, our results showed that although there was only a small number of intersections between the DETEs and DEGs, this number is still significantly higher than random expectation in almost all tissues and processes (supplementary table S1, Supplementary Material online), indicating some DETEs can be directly involved in regulating nearby gene expression during adaptation.

Table 1

Numbers of DEGs, DETEs, and Genes Within Upstream 5,000 bp Near DETEs

			No. of DEGs	No. of DETEs	No. of Genes Upstream of TEs	No. of Overlap	Ratio (%)[c]
Brain	F[a]	PO	911	569	103	24	2.63
		GN	384	396	84	10	2.6
	R b	PO	622	515	122	12	1.93
		GN	351	560	115	14	3.99
Heart	F	PO	202	137	36	4	1.98
		GN	1,588	1175	226	110	6.93
	R	PO	1,525	1139	248	104	6.82
		GN	430	576	82	20	4.65
Liver	F	PO	398	254	57	15	3.77
		GN	974	810	171	52	5.34
	R	PO	353	263	72	10	2.83
		GN	699	351	102	41	5.87
Lung	F	PO	647	630	154	48	7.42
		GN	1,277	1426	339	109	8.54
	R	PO	634	964	261	37	5.84
		GN	751	675	156	41	5.46
Muscle	F	PO	154	105	23	5	3.25
		GN	1,352	672	137	64	4.73
	R	PO	342	115	30	8	2.34
		GN	391	204	51	9	2.3

F means the forward adaptation.

R means the reverse adaptation.

The ratio equals the number of overlapping genes divided by the number of DEGs.

In addition to directly affecting the expression of nearby genes, the DETEs may form expression networks with other protein-coding genes, thereby playing certain biological functions. To further explore the potential roles of TEs in reverse adaptation, we selected all the expressed genes and DETEs of all tissues and determined the correlation between genes and DETEs using the weighted gene co-expression network analysis (WGCNA). This analysis revealed multiple interesting modules for each tissue, with module sizes larger than 100 (106–1,274). Based on the map of module-trait relationships, we selected one representative module (co-expression network) for each tissue (P ≤ 0.05; supplementary fig. S1, Supplementary Material online) and visualized them using the Cytoscape (Shannon et al. 2003). Top 20 hub genes/TEs that could play an important role (Das et al. 2017) were identified within the five interest modules based on gene connection degrees in the co-expression network using the cytoHubba plugin (fig. 2). We found that diverse DETEs served as hub TEs and participated in the co-expression network with the other genes. These results showed that diverse DETEs might indirectly participate in the co-expression network to ease the phenotypic plasticity and, ultimately, adaptation.

Fig. 2.

Top 20 hub genes/DETEs in the interest modules. A node represents a gene/TE, and its size is directly proportional to the number of other gene/TE connections in the module. The orange color represents TEs, whereas other colors represent genes. (A) Top 20 hub genes/TEs in the salmon module that is positively correlated with PO of brain in forward adaptation. (B) Top 20 hub genes/TEs in the blue module that is positively correlated with PO of heart in forward adaptation. (C) Top 20 hub genes/TEs in the brown module that is positively correlated with PO of liver in forward adaptation. (D) Top 20 hub genes/TEs in the green module that is positively correlated with PO of lung in forward adaptation. (E) Top 20 hub genes/TEs in the magenta module that is positively correlated with PO of muscle in forward adaptation.

To investigate the similarity in the function of genes with DETEs within interest modules, we performed a functional enrichment analysis of all genes in the modules with the DAVID annotation method (Huang et al. 2009). In total, about 101 biological processes (BPs) of gene ontology (GO) terms were significantly enriched in the five modules (table 2). Interestingly, all the GO terms of the tissues, except the muscles, refer to the regulatory process. Analyses showed that the regulatory process mainly include negative regulation of transcription, DNA templated; positive regulation of cartilage development; regulation of transcription, DNA templated; and negative regulation of gene expression, positive regulation of the apoptotic process (table 2). These results further indicate that the DETEs play a regulatory role in prompting the phenotypic plasticity when an organism returns to the ancestral environment.

Table 2

The GO Terms Enriched by the Interested Module Genes

Tissue	Module	Term ID	Terms	P-Value
Brain	Salmon	GO:0055003	Cardiac myofibril assembly	0.03149
		GO:0055010	Ventricular cardiac muscle tissue Morphogenesis	0.04904
		GO:0032781	Positive regulation of ATP (Adenosine Triphosphate)ase activity	0.04904
Heart	Blue	GO:0045669	Positive regulation of osteoblast differentiation	5.45E−04
		GO:0030308	Negative regulation of cell growth	6.29E−04
		GO:0032922	Circadian regulation of gene expression	6.94E−04
		GO:0006351	Transcription, DNA templated	0.00114
		GO:0001938	Positive regulation of endothelial cell proliferation	0.00137
		GO:0043065	Positive regulation of apoptotic process	0.00408
		GO:0010629	Negative regulation of gene expression	0.00476
		GO:0001516	Prostaglandin biosynthetic process	0.00485
		GO:0061036	Positive regulation of cartilage development	0.00485
		GO:0045892	Negative regulation of transcription, DNA templated	0.00504
Liver	Brown	GO:0060324	Face development	3.52E−04
		GO:0034097	Response to cytokine	9.93E−04
		GO:0045766	Positive regulation of angiogenesis	0.00137
		GO:0061036	Positive regulation of cartilage development	0.0024
		GO:0006814	Sodium ion transport	0.0024
		GO:0009408	Response to heat	0.00305
		GO:0000165	MAPK (Mitogen-Activated Protein Kinase) cascade	0.00305
		GO:0043032	Positive regulation of macrophage activation	0.00384
		GO:0045198	Establishment of epithelial cell apical/basal polarity	0.00384
		GO:0045603	Positive regulation of endothelial cell differentiation	0.00629
Lung	Green	GO:0043484	Regulation of RNA splicing	0.00395
		GO:0043280	Positive regulation of cysteine-type endopeptidase activity involved in apoptotic process	0.00395
		GO:0006351	Transcription, DNA templated	0.00731
		GO:0008380	RNA splicing	0.00774
		GO:0033146	Regulation of intracellular estrogen receptor signaling pathway	0.0092
		GO:0009968	Negative regulation of signal transduction	0.02099
		GO:0006355	Regulation of transcription, DNA templated	0.02905
		GO:0019985	Translesion synthesis	0.03659
		GO:0006397	mRNA processing	0.03889
		GO:0045944	Positive regulation of transcription from RNA polymerase II promoter	0.04724
Muscle	Magenta	GO:0001944	Vasculature development	0.02845
		GO:0007275	Multicellular organism development	0.03621
		GO:0048821	Erythrocyte development	0.04436

Previous studies have shown that organisms can “remember” their ancestral environments (evolutionary memory), thus affecting their future evolution (Desai 2009; Ho et al. 2020). Our study based on the analyses of the TEs further provides evidence that the plastic expression changes, not only in protein-coding genes but also in TEs (maybe also involve other non-coding regions), can be highly involved in the phenotypic plasticity to ease organismal re-adaptations to ancestral environments. Although most transposition events of TEs may be harmful to the host by reducing the genomic stability, multiple TEs insertions can be transcribed and subsequently participated in the gene regulation landscape of host adaptive evolution. Our study, therefore, recommends detailed molecular and cellular experiments to further explore the regulatory activities of TEs.

Materials and Methods

Differential Expression Analysis

We used the SQuIRE (Software for Quantifying Interspersed Repeat Expression; Yang et al. 2019) to measure the TEs and gene expression levels of the PO and GN processes during forward and reverse adaptation. The SQuIRE is a powerful tool that performs differential expression of genes and TEs (Yang et al. 2019). The analysis process consists of three stages: preparation, quantification, and analysis, with a total of procedures: Fetch, Clean, Map, Count, and Call (https://github.com/wyang17/SQuIRE). The DEGs and DETEs were identified using a threshold of a P-value = 0.05 from the SQuIRE_Call results. For each DETE, we evaluated the relative importance (measured as the ratio of PO-DETEs divided by GN-DETEs) of plastic and genetic changes in its expression level.

WGCNA

The WGCNA is a valuable tool for finding modules of highly correlated genes by clustering similar expressions of genes. All genes (n = 6,885) and DETEs in the five tissues (brain, heart, liver, lung, and muscles) were selected for the WGCNA analysis. The general input data include the standardized gene expression matrix and traits matrix (the trait here refers to the different stages), we used the FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) expression matrix and numerical characteristic matrix here. The FPKM expression matrix originates from the SQuIRE_Count result. Due to the input limit, genes/DETEs with low expression (FPKM value ≤0) were filtered out. Sequent analysis was performed for the remaining genes and DETEs. Subsequent analyses were performed with the WGCNA package (Langfelder and Horvath 2008) in R. Based on the criterion of approximate scale-free topology, we chose the proper soft thresholding power for each adaptive process of each tissue using the built-in function “pickSoftThreshold.” Other important parameters were set as default settings: “min module size ≥30” and “merge cut height = 0.25.” The interesting modules for each tissue were selected for sequent analysis, with module sizes larger than 100 and P ≤ 0.05 based on the map of module-trait relationships. Top 20 hub genes/TEs within interest modules were identified by calculating intra-modular connectivity using “cytoHubba plug-in” in the Cytoscape (version. 3.9.0).

Enrichment Analysis

We used the DAVID (Database for Annotation, Visualization and Integrated Discovery) tool to perform enrichment analysis of all interested genes. We chose the significantly differential GO terms with P ≤ 0.05 from the BP.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.