Transposable Element Tissue-Specific Response to Temperature Stress in the Stenothermal Fish Puntius tetrazona

Simple Summary Teleosts are one of the most diversified group of vertebrates, as they colonized different aquatic environments worldwide. Stenothermal fish are characterized by a limited ability to tolerate temperature changes, therefore variations in this abiotic factor represent a threat for their survival. In the present study, we selected the tiger barb Puntius tetrazona, for which RNA-Seq data were available for brain, gill, and liver, to investigate the possible role of the transcriptional activity of transposable elements (TEs) in the rapid adaptation of this species. Our findings highlighted a tissue-specific response of TEs with a remarkable increase at 13 °C recorded in liver, strengthening the view of TEs as source of genetic variability acting positively in species resilience. Abstract Ray-finned fish represent a very interesting group of vertebrates comprising a variety of organisms living in different aquatic environments worldwide. In the case of stenothermal fish, thermal fluctuations are poorly tolerated, thus ambient temperature represents a critical factor. In this paper, we considered the tiger barb Puntius tetrazona, a freshwater fish belonging to the family Cyprinidae, living at 21–28 °C. We analyzed the available RNA-Seq data obtained from specimens exposed at 27 °C and 13 °C to investigate the transcriptional activity of transposable elements (TEs) and genes encoding for proteins involved in their silencing in the brain, gill, and liver. TEs are one of the tools generating genetic variability that underlies biological evolution, useful for organisms to adapt to environmental changes. Our findings highlighted a different response of TEs in the three analyzed tissues. While in the brain and gill, no variation in TE transcriptional activity was observed, a remarkable increase at 13 °C was recorded in the liver. Moreover, the transcriptional analysis of genes encoding proteins involved in TE silencing such as heterochromatin formation, the NuRD complex, and the RISC complex (e.g., AGO and GW182 proteins) highlighted their activity in the hepatic tissue. Overall, our findings suggested that this tissue is a target organ for this kind of stress, since TE activation might regulate the expression of stress-induced genes, leading to a better response of the organism to temperature changes. Therefore, this view corroborates once again the idea of a potential role of TEs in organism rapid adaptation, hence representing a promising molecular tool for species resilience.


Introduction
The evolutionary success of a species is strictly related to their ability to cope with changes in environmental conditions. Phenotypic plasticity is the main mechanism used by organisms to respond to rapid environmental perturbations. Indeed, it can produce well-adapted phenotypes that, in turn, might improve the fitness of organisms [1,2]. Phenotypic plasticity is also the result of environmental sensitivity of the genome that, through modifications in the heterochromatin structure, varies the expression not only of genes first time an activation of TE transcriptional activity in response to salinity variations in the marbled eel Anguilla marmorata [11]. Since most teleosts are poikilotherms, environmental stimuli affect physiological and metabolic activities in these organisms [34][35][36][37]. In particular, ambient temperature might be a critical factor for stenothermal fish in which thermal fluctuations are less tolerated. The tiger barb Puntius tetrazona is a freshwater fish belonging to the family Cyprinidae, living at 21-28 • C. In this paper, we analyzed the RNA-Seq data of tiger barb exposed at 27 • C and 13 • C [38] to investigate the transcriptional activity of TEs and genes encoding for proteins involved in their silencing in the brain, gill, and liver. Our results pointed out a clear response of TEs in the liver consistent with the idea of a potential role of these elements in the rapid adaptation, thus representing a promising molecular tool for species resilience. Furthermore, all TE silencing mechanisms were active in the hepatic tissue and ago genes also in the brain and gill tissues.

Materials and Methods
Tiger barb RNA-Seq raw data were obtained from the Sequence Read Archive (SRA) (https://www.ncbi.nlm.nih.gov/sra, accessed on 15 June 2022) under the accession number SRP153005 [38]. Data were obtained sequencing the brain, gill, and liver of ninety fish (eight-month-old) divided into two groups, one used as the control group (27 • C) and another for testing the acute cold stress until 13 • C. The temperature was set to 27 ± 0.5 • C on the first day, followed by a decrease of 2 • C every 24 h. Fish were kept at 13 ± 0.5 • C for 24 h, and then, five fish, in triplicate, in each tank, were anesthetized before dissection. Raw paired-end reads were imported in the CLC Genomics Workbench v.12 (Qiagen, Hilden, Germany) (Table S1) and trimmed removing sequencing adapters, low-quality bases, and low-quality read ends using default parameters. Trimmed reads were then de novo assembled with the "De Novo Assembly" of CLC Genomics Workbench v.12 using the default parameters. Completeness of the de novo assembled transcriptome was evaluated through BUSCO v.5 using the Actinopterygii OrthoDB v.10 database as a reference [39].

Estimation of TE Transcriptional Activity
To estimate the TE transcriptional activity, we first identified TEs in the de novo assembled transcriptomes with RepeatMasker v.4.1.0 (http://www.repeatmasker.org/cgibin/WEBRepeatMasker, accessed on 15 July 2022) using a custom library created following the methodology described in our previous work [30]. After the RepeatMasker analysis, the trimmed reads related to 27 • C and 13 • C were mapped against the reference transcriptome to calculate the expression values using the proprietary RNA-Seq Analysis tool included in the CLC Genomics Workbench v.12 and setting the following mapping parameters: length fraction = 0.75 and similarity fraction = 0.98. To remove redundancies, the RepeatMasker output file was filtered, removing entries not classified as TEs and keeping those with the highest score and length values. The overall expression of each TE class was calculated by summing the expression values of each TE type: DNA transposons, LINE, LTR (also including endogenous retroviruses), non-LTR, Retro (retroelements that are not classified in any of the two main subclasses), and SINE. In our results, we also reported the fraction of unclear elements that are referred to repetitive elements without specific features to determine their attribution to a given TE class. We included this fraction to show the low percentage of the unclear fraction, representative of the goodness of our analyses. The expression values were then transformed into a percentage of mapped reads to achieve the comparability between species.

Identification and Expression Enalysis of Genes of Interest
Genes of interest were retrieved through TBLASTN [40] from the assembled transcriptome obtained from the brain, gill, and liver of P. tetrazona. In particular, the search was made for genes involved in heterochromatinization (chromobox homolog 5 (cbx5), chromobox homolog 1 (cbx1), chromobox homolog 3 (cbx3), DNA (cytosine-5-)-methyltransferase 1 (dnmt1), DNA (cytosine-5-)-methyltransferase 3 alpha (dnmt3a), and SET domain bifurcated histone lysine methyltransferase 1 (setdb1)); genes related to the NuRD complex (krab-like, trim33, chromodomain helicase DNA binding protein 3 (chd3), chromodomain helicase DNA binding protein 4 (chd4), histone deacetylase 1 (hdac1), methyl-CpG binding domain protein 2 (mbd2), methyl-CpG binding domain protein 3 (mbd3), metastasis associated 1 (mta1), metastasis associated 1 family, member 2 (mta2), metastasis associated 1 family, member 3 (mta3), GATA zinc finger domain containing 2 (gatad2), and retinoblastoma binding protein 4 (rbbp4), retinoblastoma binding protein 7 (rbbp7)); four genes of the Argonaute subfamily (argonaute RISC component 1 (ago1), argonaute RISC component 2 (ago2), argonaute RISC component 3 (ago3), and argonaute RISC component 4 (ago4)); and three genes of the GW182 family proteins (tnrc6a, tnrc6b, and tnrc6c). Transcripts were translated using the EMBOSS Transeq translation tool (https://www.ebi.ac.uk/Tools/st/emboss_transeq/, accessed on 12 September 2022), and the UTR and CDS regions were identified (Table S2). The aforementioned sequences were deposited in GenBank under the accession numbers listed in Table S2. To ensure the comparison between species, the gene expression values were computed using a scaling factor based on the cumulative expression of a dataset composed of 2124 orthologs derived from the Actinopterygii OrthoDB v.10 database [39]. In detail, the dataset was created as follows: for each transcriptome of the tissue considered in this study, expression levels of the genes attributed to the BUSCO analyses as "complete and single copy" and "fragmented" genes were kept as the number of mapped reads; expression values of genes classified as "duplicated" (most probably derived from transcriptional isoforms) were computed as the sum of each copy of single BUSCO, and the expression levels of "missing" genes were set to 0. This dataset was then used as a calibration set, computing a scaling factor that was applied to the original expression values of the genes of interest, as described in Biscotti et al. (2016) [41]. Transcriptional values of the genes analyzed in this study were reported as Transcripts Per Million (TPM).

Statistics
For each tissue, the data on gene expressions and TE transcriptional activity obtained from the three replicates for each tested condition were expressed as the mean ± standard error, and statistically significant differences were evaluated by one-way ANOVA. The symbol * indicates p-values < 0.05, ** for p-values < 0.01, and *** for p-values < 0.001.

Results
The transcriptional activity of the TEs was evaluated as a percentage of the mapped reads in the brain, gill, and liver of the tiger barb ( Figure 1). The data reported in Figure 1A were referred to the total TE transcriptional contribution at the control (27 • C) and test (13 • C) conditions. Comparing the total TE transcriptional activity at 27 • C between the analyzed tissues, the highest level was showed in the brain, followed by the gill and then liver. At 13 • C, the highest TE transcriptional activity in the brain was confirmed, and the difference between the gill and liver was less remarkable. Moreover, analyzing the data obtained in each tissue, a statistically significant difference between the two conditions considered here, was detectable only in the liver (p-value: 1.65 × 10 −5 ), in which the total TE transcriptional activity was increased at 13 • C. Profiles of the TE relative abundances between the brain and gill samples were similar. Overall, a prevalence of DNA transposons followed by LINE, unclear, Retro, LTR, SINE, and non-LTR retroelements was shown. In the liver, LINE retroelements prevailed at 27 • C, while the TE distribution pattern was similar to that observed in the other two tissues at 13 • C. Despite the absence of an appreciable variation in the total contribution of TEs in the brain and gill, statistically significant differences emerged for the same tissues in single TE types between the two tested conditions. Indeed, in the brain, the LINE, LTR, and SINE retroelements showed a significant decrease passing from 27 • C to 13 • C ( Figure 1B). In the gill, except for unclear elements and SINE retroelements, all TE types experienced significant changes, with DNA transposons and Retro having a higher activity at 13 • C, while LINE, LTR, and non-LTR retroelements have a lower activity in this condition than in the control one ( Figure 1C). In the hepatic tissue, with the exception of LINE retroelements that showed lower values at 13 • C, all the analyzed TE types significantly increased their transcriptional contribution ( Figure 1D). These increases were ascribable to a few specific elements. Indeed, the ten elements having the highest values of variation in transcriptional activity between the two conditions represented 62.27%, 86.0%, and 73.35% for the DNA transposons, SINE, and LTR retroelements, respectively. In particular, two elements, DNA/hAT-Ac and LTR/DIRS, mostly contributed to the upregulation of their belonging TE type (32.67% and 51.49%, respectively). Moreover, in the other tissues, the ten elements having the highest values of variation in transcriptional activity between the two conditions were different from those identified in the liver, except for one LTR and four SINE retroelements.
differences emerged for the same tissues in single TE types between the two tested conditions. Indeed, in the brain, the LINE, LTR, and SINE retroelements showed a significant decrease passing from 27 °C to 13 °C ( Figure 1B). In the gill, except for unclear elements and SINE retroelements, all TE types experienced significant changes, with DNA transposons and Retro having a higher activity at 13 °C, while LINE, LTR, and non-LTR retroelements have a lower activity in this condition than in the control one ( Figure 1C). In the hepatic tissue, with the exception of LINE retroelements that showed lower values at 13 °C, all the analyzed TE types significantly increased their transcriptional contribution ( Figure 1D). These increases were ascribable to a few specific elements. Indeed, the ten elements having the highest values of variation in transcriptional activity between the two conditions represented 62.27%, 86.0%, and 73.35% for the DNA transposons, SINE, and LTR retroelements, respectively. In particular, two elements, DNA/hAT-Ac and LTR/DIRS, mostly contributed to the upregulation of their belonging TE type (32.67% and 51.49%, respectively). Moreover, in the other tissues, the ten elements having the highest values of variation in transcriptional activity between the two conditions were different from those identified in the liver, except for one LTR and four SINE retroelements. The transcriptional activity of TEs is controlled by the host genome through silencing mechanisms that determine a stronger heterochromatinization status through the activity of heterochromatin proteins, DNA methyltransferases, and the NuRD complex, and AGO proteins also involved in mRNA cleavage and translation repression. In this study, a total of 34 genes encoding proteins involved in TE silencing mechanisms were retrieved: ten concerning heterochromatin formation (cbx5, cbx1a, cbx1b, cbx3a, cbx3b, dnmt1, dnmt3aa,  dnmt3ab, setdb1a, and setdb1b); 16 related to the NuRD complex (krab-like, trim33, chd3, The transcriptional activity of TEs is controlled by the host genome through silencing mechanisms that determine a stronger heterochromatinization status through the activity of heterochromatin proteins, DNA methyltransferases, and the NuRD complex, and AGO proteins also involved in mRNA cleavage and translation repression. In this study, a total of 34 genes encoding proteins involved in TE silencing mechanisms were retrieved: ten concerning heterochromatin formation (cbx5, cbx1a, cbx1b, cbx3a, cbx3b, dnmt1, dnmt3aa, dnmt3ab, setdb1a, and setdb1b); 16 related to the NuRD complex (krab-like, trim33, chd3, chd4a, chd4b, hdac1, mbd2, mbd3a, mbd3b, mta1, mta2, mta3, gatad2ab, gatad2b, rbbp4, and rbbp7); four to ago (ago1, ago2, ago3, and ago4); and four belonging to the GW182 family proteins (tnrc6a, tnrc6b, tnrc6c1, and tnrc6c2). Of these, 20 showed a complete CDS, 2 were incomplete, and 12 fragmented gene sequences were replaced with the complete version available in public repositories to further refine the assembled transcriptome (Table S2).
Regarding the transcriptional activity, our findings pointed out that all these mechanisms were active in P. tetrazona (Figures 2-4). In the case of the brain, a significant increase Animals 2023, 13, 1 6 of 11 in the transcription of the genes involved in the heterochromatin formation was detected for cbx3a, as well as the three dnmt (dnmt1, dnmt3aa, and dnmt3ab), while a significant decrease was shown for cbx3b ( Figure 2A). No remarkable change was detected comparing the expression of the genes encoding for the proteins of the NuRD complex between the two considered temperatures ( Figure 2B) differently from the ago genes and those of the GW182 family proteins that recorded a statistically significant increase ( Figure 2C,D). Analyzing the RNA-Seq data of a gill, a general statistically significant decrease in the transcriptional level of the genes involved in heterochromatin formation (cbx1a, cbx3b) and the NuRD complex (krab-like, chd3, chd4b, and mta3) was observed, while, in the latter, a significant increase was observed for mta1 ( Figure 3A,B). Differences in the expression of ago genes and in those encoding the GW182 family proteins were identified also in the gill ( Figure 3C,D). In the hepatic tissue, most of the analyzed genes showed higher values at 13 • C (Figure 4). In particular, cbx5, cbx1b, cbx3b, dnmt1, dnmt3aa, dnmt3ab, trim33, chd4a, hdac1, gatad2ab, gatad2b, rbbp4, the four ago genes, and two of the GW182 family proteins showed significant changes. chd4a, chd4b, hdac1, mbd2, mbd3a, mbd3b, mta1, mta2, mta3, gatad2ab, gatad2b, rbbp4, and  rbbp7); four to ago (ago1, ago2, ago3, and ago4); and four belonging to the GW182 family proteins (tnrc6a, tnrc6b, tnrc6c1, and tnrc6c2). Of these, 20 showed a complete CDS, 2 were incomplete, and 12 fragmented gene sequences were replaced with the complete version available in public repositories to further refine the assembled transcriptome (Table S2).
Regarding the transcriptional activity, our findings pointed out that all these mechanisms were active in P. tetrazona (Figures 2-4). In the case of the brain, a significant increase in the transcription of the genes involved in the heterochromatin formation was detected for cbx3a, as well as the three dnmt (dnmt1, dnmt3aa, and dnmt3ab), while a significant decrease was shown for cbx3b ( Figure 2A). No remarkable change was detected comparing the expression of the genes encoding for the proteins of the NuRD complex between the two considered temperatures ( Figure 2B) differently from the ago genes and those of the GW182 family proteins that recorded a statistically significant increase ( Figure  2C,D). Analyzing the RNA-Seq data of a gill, a general statistically significant decrease in the transcriptional level of the genes involved in heterochromatin formation (cbx1a, cbx3b) and the NuRD complex (krab-like, chd3, chd4b, and mta3) was observed, while, in the latter, a significant increase was observed for mta1 ( Figures 3A,B). Differences in the expression of ago genes and in those encoding the GW182 family proteins were identified also in the gill ( Figure 3C,D). In the hepatic tissue, most of the analyzed genes showed higher values at 13 °C (Figure 4). In particular, cbx5, cbx1b, cbx3b, dnmt1, dnmt3aa, dnmt3ab, trim33, chd4a,  hdac1, gatad2ab, gatad2b, rbbp4, the four ago genes, and two of the GW182 family proteins showed significant changes. Values are reported for both tested temperature conditions (27 °C and 13 °C). Statistically significant differences are presented as * for p < 0.05, ** for p < 0.01, and *** for p < 0.001. Values are reported for both tested temperature conditions (27 • C and 13 • C). Statistically significant differences are presented as * for p < 0.05, ** for p < 0.01, and *** for p < 0.001. Values are reported for both tested temperature conditions (27 °C and 13 °C). Statistically significant differences are presented as * for p < 0.05, ** for p < 0.01, and *** for p < 0.001. Values are reported for both tested temperature conditions (27 • C and 13 • C). Statistically significant differences are presented as * for p < 0.05, ** for p < 0.01, and *** for p < 0.001. Values are reported for both tested temperature conditions (27 °C and 13 °C). Statistically significant differences are presented as * for p < 0.05, ** for p < 0.01, and *** for p < 0.001.

Discussion
The ability of organisms to face changes of biotic and abiotic factors is mainly due to phenotypic plasticity that also requires modifications at the genome level. In particular, environmental changes can cause variations in the epigenetic status, leading to gene activation but also to the impairment of transposon silencing mechanisms inducing the activation of these mobile elements with the possible creation of genetic variability. Therefore, TEs represent a powerful adaptive response to environmental perturbations [7,8,27,42,43]. Although this issue has been largely studied in plants [20][21][22]25,[44][45][46][47], information is still scarce in animals and, in particular, in vertebrates [27,43,48].
In this paper, we investigated the transcriptional TE response and the activity of genes involved in their silencing in the stenothermal fish P. tetrazona, a very popular ornamental freshwater fish native to Southeast Asia in places such as Malaysia and Borneo, places characterized by an equatorial climate; some of them are found in Thailand, Sumatra Island, and Cambodia [49]. This fish, living at 21-28 • C, tolerates small thermal fluctuations and experiences serious histopathological damages at low temperatures [38]. Overall, ambient temperature is an abiotic factor that acts on physiological and metabolic activities of poikilothermic teleosts, and in particular, on those of stenothermal fish [34][35][36][37]. Transcriptional analyses conducted by Liu and colleagues (2020) [38] on RNA-seq data obtained from P. tetrazona exposed at 27 • C (control) and 13 • C (test) have showed a high number of differentially expressed genes: in brain and liver the upregulated genes prevailed, differently from gill, in which the downregulated genes were the most represented. This is in line with the general idea that the genome epigenetic status varies in response to a stress, leading to changes in the expression levels. Our analysis of TE transcriptional contribution, performed on the same datasets, highlighted a different response in the three analyzed tissues. Apparently, the brain and gill showed the same behavior: in both the transcriptional activity of TEs did not vary between the two tested conditions. However, in the cerebral tissue, the activity of genes involved in heterochromatin formation, as well as ago genes and those related to GW182 family proteins, presented higher transcriptional levels at the stress condition. The lack of variation in TE transcriptional contribution in this tissue might be related to the silencing activity of these genes or to the short exposure period. In gill, although most of genes related to TE silencing mechanisms decreased their activity at 13 • C, no variation was observed in total TE transcriptional activity. However, for DNA transposons and Retro retroelements a slight increase of expression levels was appreciable at the stress condition. This might indicate that a chronic exposure is necessary to gain a remarkable TE variation.
Interestingly, a statistically significant increase of TE transcriptional activity was showed only in liver when fish were exposed at 13 • C. Our findings revealed that this increase was mainly due to few specific TEs that seem to undergo a remarkable transcriptional variation only in the hepatic tissue. This is in line with previous papers, reporting a relationship between a kind of stress and a specific TE type response. The mPing DNA transposon has been reported to be activated in response to salt stress in rice [44,46], the ONSEN retrotransposon as consequence of heat stress in Arabidopsis [45,47], and the nuclear import of Tam3 transposon related-transposase has been showed to increase after exposition to temperature decrease in Antirrhinum [20][21][22]25]. In addition, the different TE behavior might be related to their tissue-specific response. Indeed, Hunter and colleagues (2012) [50] have reported a stress-induced activation of particular TE families as ERV intracisternal-A particle (IAP), B2_RN SINE, and L1_RN in the rat hippocampus.
It is well known that the transcriptional and transpositional activity of TEs can have positive effects increasing genetic variability and leading to advantageous phenotypes. On the other side, TEs can provoke negative effects, deleterious for the host genome that in turn activates TE silencing mechanisms. In the tiger barb P. tetrazona, the high TE transcriptional levels at 13 • C in liver were accompanied by an increase in the transcriptional activity of genes related to heterochromatinization, the NuRD complex, ago, and the GW182 family proteins. This picture might suggest that the hepatic tissue has a quicker responsiveness than other tissues. Indeed, it is a target organ for this kind of stress since several metabolic activities, also linked to body thermoregulation, occur in liver. Species respond to thermal conditions by modulating the metabolic rate to satisfy energy demand [51]. In P. tetrazona liver, Liu and colleagues (2020) [38] have reported an upregulation of cold-induced genes belonging to metabolism and biosynthesis pathways such as biosynthesis of steroid and fatty acids and tryptophan metabolism. The role of the hepatic tissue is extremely important since the tiger barb is a poikilothermic organism and also a stenothermal fish, tolerating small thermal fluctuations. The activation of pathways associated with organism response to stress conditions might be related to TE transcriptional activity. It is well-known that these genetic elements contribute in the up-or downregulation of nearby genes through cisor trans-regulatory sequences located within TEs or TE-derived noncoding RNAs [26,52,53]. Alternatively, the decrease of methylation levels of genes involved in stress response might passively affect the heterochromatic state of nearby TEs leading to an increase of their transcriptional activity and consequently genes involved in their silencing are activated to counteract this epigenetic deregulation [54]. On the other side, liver might better tolerate the negative mutagenic effects related to TE activation thanks to its ability to regenerate after cell damage.

Conclusions
In conclusion, the analysis of the transcriptional activity of TEs and genes related to their silencing highlighted a clear response of these mobile elements in liver. Overall, the tissue-specific activation of TEs might represent a promising molecular tool, favoring species adaptation and resilience. Knowledge in this field could be useful in the perspective of wildlife species conservation since climate changes are threatening their biodiversity.