Abstract
Buergeria frog tadpoles exhibit high thermal tolerance and are occasionally found in water pools that temporarily exceed 40°C. With the aim of understanding how they can cope with the severe heat stress, we performed RNA-seq of three heat-treated (38°C) and three control (25°C) tadpoles and compared their transcriptomic profiles. We identified 382 differentially expressed transcripts. A protein-protein interaction (PPI) network analysis of these transcripts further identified hub proteins involved in protein degradation, stress granule assembly, and global suppression of DNA transcription and mRNA translation. Along with the avoidance behavior against high water temperature, these endurance mechanisms potentially support tadpoles to survive in high temperatures for short periods of time. Similar mechanisms may exist in many other amphibian species whose habitats are prone to high temperatures.
Introduction
In their natural environment, organisms face various environmental stressors such as ultraviolet light, low and high pH, osmotic pressure, freezing, and heat. For survival, organisms exhibit various responses to such stressors at multiple biological levels, including changes in morphology, physiology, behavior, and life history (; ; ). Most of these responses are initially caused by intracellular regulations mediated by the neural and hormonal system (; ; ; ). Therefore, these stress responses become the paradigm of how molecular changes influence individual survival under environmental challenges (). Thus, a better understanding of how organisms cope with extreme stressors can reveal important molecular mechanisms enabling them to survive in such environments. In particular, thermal conditions pose significant challenges to ectotherms. Without thermoregulatory mechanisms, the survival and distribution of organisms are largely determined by temperature, whereby an appropriate response to thermal stress is essential for their survival.
Tree frog species of the genus Buergeria are suitable for studying how organisms cope with extreme thermal stress because several populations live in extreme hot environments. Buergeria frogs are distributed in a wide range of habitats, from mountain streams to coastal areas across continental islands in Eastern and Southeastern Asia (; ). Currently, six species have been identified, and several populations of three species, B. otai, B. choui, and B. japonica, are further found in geothermal hot springs (; ; ; ).
In geothermal hot springs, the water temperature is non-uniform and varies temporally and spatially due to rainfall, water depth, air temperature, and distance from the hot spring source. Field and laboratory observations revealed that Buergeria tadpoles and adults from hot spring populations tended to avoid inhabiting and reproducing in water temperatures above ∼35°C (; ). Avoidance behaviors against heat stress are probably essential for populations to survive in hot spring environments. However, due to the nature of hot springs, the water temperature can fluctuate, and tadpoles must face inevitable temperature spikes in this fully aquatic life stage. Indeed, Buergeria frogs living in hot springs can tolerate hot water: a previous study measured 30 Buergeria tadpoles collected from a hot spring habitat and revealed tolerance of temperatures up to 46°C, whereas above this temperature, tadpoles lost balance and swimming ability or died (). Such high thermal tolerance has been observed not only in hot spring populations but also in freshwater Buergeria populations across their distribution (). The high thermal tolerance necessary for survival at these temperatures demands protective transcriptional and subsequent cellular responses to such stress. These responses, together with behavioral strategies, are expected to play essential roles in Buergeria frogs inhabiting hot spring environments.
Transcriptional and proteomic responses of anurans to cold or freezing temperatures have been previously characterized, and genes and proteins related to heat shock, ubiquitin-proteasome, antioxidation, and glycogenolysis were identified in response to low-temperature stress (; ; ). However, to the best of our knowledge, transcriptome analysis related to heat stress has not been reported, presumably due to the lack of well-known heat-tolerant frog species. Therefore, Buergeria frogs, which are known to inhabit hot springs and exhibit high thermal tolerance, are suitable subjects for studying the transcriptional response to explicit heat stress.
Given their avoidance of high-temperature waters (; ), Buergeria frogs may be less exposed to long-term heat stress. Rather, it may be more important for these frogs to endure temporary heat stress while moving away from the stressor. In this study, to further understand the molecular mechanisms and strategies that function in response to heat stress, we performed a comparative transcriptome analysis of heat-treated and untreated B. otai tadpoles.
Materials and Methods
Sample Preparations
We collected B. otai adult females and males in Hualien, Taiwan in August 2013. The adult individuals were kept in a laboratory where they then laid eggs. It is unclear whether the eggs were laid by a single female-male pair or not. Eggs were transferred to approximately 25°C water tanks (∼50 L) with filtration systems and 9:15 light-dark cycle. The water temperature was chosen according to previous studies (; ). Hatched tadpoles were kept under the same condition until use in thermal experiments. The water was not replaced. Tadpoles were fed boiled spinach once a day.
Six tadpoles at Gosner stages 26–28 () were moved into two 1000-mL glass beakers containing 500-mL of 25°C water and were maintained for a day in incubators. After an additional 12 h fasting, three tadpoles were directly moved into a glass beaker with 38°C water which is slightly cooler than their possible thermal limit (∼40°C) (). We confirmed that all three tadpoles did not exhibit any abnormal behavior after a sudden temperature change. The other three tadpoles were moved into a beaker with 25°C water. We maintained the water temperatures of each beaker for approximately 4 h using incubators, and then tadpoles were euthanized by pithing immediately after they were taken out from the water. Bodies were transferred into a RNAlater Stabilization Solution (Thermo Fisher Scientific, MA, United States).
Total RNA was extracted from the whole-body using a RNeasy mini kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. RNA samples were shipped to BGI (Beijing, China), where RNA-seq was performed using Illumina Hi-Seq 2000 (Illumina, CA, United States).
Since amphibian tadpoles were not subject to ethical review in Hiroshima University where the experiment was conducted at that time, no institutional ethical approval was obtained. This study was carried out in compliance with the ARRIVE (Animal Research: Reporting of In Vivo Experiments) guidelines.1 In the experimental procedures using living tadpoles, we made an effort to minimize stress of individuals following Guidelines for Live Amphibians and Reptiles in Field and Laboratory Research developed by the American Society of Ichthyologists and Herpetologists in 2004.
Transcriptome Analyses
Raw RNA-seq data went through read filtering based on base quality and read length using Trim Galore v0.6.4 with default settings (phred score threshold of 20 and length threshold of 20 bp). Then, de novo assembly was performed using Trinity v2.8.3 (), and summary statistics of the assembly were calculated by SeqKit v0.13.2 (). Completeness of the assembly was assessed using BUSCO version 5.0.0 () specifying tetrapoda_obd10 as a reference. Pseudoalignment and transcript quantification was performed using kallisto v0.45.1 () with 100 bootstraps. Output files of kallisto were imported into the R environment using the R package “tximport” v1.10.1 (), and counts of reads mapped to each contig were normalized by the trimmed means of M (TMM) values method. A differential expression analysis was then performed using the R package “edgeR” v3.24.3 () in R 3.6.3 (R Foundation). A likelihood ratio test based on TMM values was applied in the analysis. Transcripts with an FDR-corrected p (q-value) ≤ 0.05 and an absolute value of log2 fold change ≥ 3 were defined as differentially expressed transcripts.
Annotation of Differentially Expressed Transcripts
To annotate differentially expressed transcripts, we first searched candidate coding regions in each contig of differentially expressed transcripts using Transdecoder v5.5.0 with default settings and the “–single_best_only” option so that only the single best candidate region per transcript was considered (). We then scanned all predicted protein sequences for homology in the Ensembl reference protein sequence of Xenopus tropicalis (v9.1) using local BLASTP v2.10.1 + (). We used the default settings and set the e-value threshold to 10–5, and only the top hits for each contig were considered in the downstream analyses.
Functional Analyses
For an overview of the transcriptomic change in response to the heat stress, we performed a protein-protein interaction (PPI) analysis and functional enrichment analyses based on transcripts that were differentially expressed and annotated by BLASTP search.
A web database STRING () was used for the PPI analysis, specifying X. tropicalis as the target organism with default settings. Annotations of both upregulated and downregulated transcripts obtained from BLASTP annotation using X. tropicalis as a reference were analyzed together in a single PPI analysis. The obtained pairwise PPIs were further analyzed with the R package “igraph” v1.2.5 () to detect hub proteins in the PPI networks. We adopted two different criteria of measurements, including degree centrality and betweenness centrality, to identify hub proteins as applied in and . The degree centrality method measures the number of interactions each protein participates in, and betweenness centrality measures how often the shortest paths connecting two other proteins pass through the given protein. Proteins with a higher score of degree centrality or betweenness centrality were regarded as more important hubs in the networks.
We conducted enrichment analyses implemented in STRING database to identify gene ontologies, pathways, keywords, domains, and features overrepresented by differentially expressed and BLASTP-annotated transcripts. The available frameworks in the STRING analysis were GO (Gene Ontology) annotations, KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways, Uniprot keywords, Reactome pathways, Pfam protein domains, SMART (Simple Modular Architecture Research Tool) protein domains, and InterPro protein features. We performed three independent analyses: upregulated transcripts, downregulated transcripts, and all differentially expressed transcripts (upregulated and downregulated transcripts). In the analyses, strength and Fisher’s exact test p-values were calculated, and the significance threshold was set to FDR-corrected p-value < 0.05. Strength is the log10-transformed number of proteins belonging to a certain category divided by the expected number of proteins belonging to the category when proteins were randomly extracted from the database and indicates the degree of enrichment.
Results
De novo Assembly and Expression Analysis
In total, 106,079,672 reads were obtained from RNA-seq of six tadpoles’ whole body and were applied for de novo assembly (Figure 1 and Supplementary Table 1). We assembled 299,588 contigs with an average length of 979.6 bp (range: 165–60,421 bp) and N50 of 2,272 bp. Estimated completeness of the assembly was 92.1% (40.7% single-copy and 51.4% duplicated BUSCOs). Fragmented and missing BUSCOs were 2.2 and 5.7%, respectively (Supplementary Table 1).
FIGURE 1
In the quantification of transcripts, 17,679,945 reads, on average, were processed per sample, 92% of which were pseudoaligned onto any contigs. Based on the log-likelihood ratio test q-value and log2 fold change thresholds from differential expression analyses, we found 147 and 235 significantly upregulated and downregulated transcripts in response to the heat stress, respectively (Figure 2 and Supplementary Table 2). There were 131 and 190 coding regions predicted by Transdecoder from 147 upregulated and 235 downregulated transcript sequences, respectively.
FIGURE 2
Annotations and Protein-Protein Interactions
The local BLASTP using X. tropicalis as a reference protein sequence annotated 86.7% (113 upregulated and 165 downregulated transcripts) of Transdecoder-predicted coding sequences. There were overlapped annotations among transcripts, and the unique number of annotations was 232.
In the STRING PPI analysis, among 232 unique genes with differentially expressed levels and BLASTP annotations, 64 had protein interactions (Figure 3). The number of protein interactions was significantly larger than expected according to the STRING network stats (PPI enrichment p = 2.2 × 10–3). In total, 13 networks were formed. Members of proteins with the highest degree centrality and betweenness centrality were similar to each other; GAK, DDX3X, and KMT2B exhibited degree centrality and betweenness centrality scores within the top five (Supplementary Table 3) and belonged to the largest network. The second-largest PPI network consisted of 10 proteins, and FBXW4 served as a hub protein.
FIGURE 3
Function of Hub and Related Proteins
Integrating the gene functions, six subgroups were identified from the PPI networks: heat shock protein (HSP), mitochondrial fission, translation regulation (×2), methyltransferase, and ubiquitination.
Heat shock protein subgroup contained three proteins of GAK, HSPA1B, DNAJA1. GAK, a hub protein in the network, is a member of the DnaJ/HSP40 family. DNAJA1, closely located to GAK on the network, is also a member of the family. HSPA1B, which had PPI interactions with GAK and DNAJA1, belongs to the HSP70. mRNA expression levels of these proteins were downregulated in response to the heat stress.
DNM1L, which had an interaction with GAK, mediates mitochondrial fission. Upregulation of this gene was also reported in a cow and mouse in response to heat stress (; ). Under the heat stress, DNM1L-induced mitochondrial fission does not play a protective role for the cell but leads to apoptosis ().
The translation regulation subgroups were found in two non-directly linked regions. One subgroup contained DDX3X, ATF4, eIF4G3, and eIF4ENIF1. DDX3X which was identified as a hub protein of this subgroup is an important component for stress granule assembly and molecular dynamics (). Stress granules are dense aggregations of RNAs and RNA-binding proteins. It is formed in the cytosol when the cell is exposed to stress and is expected to have a protective function for the cell through mRNA translation regulation including, sequestration, degradation, and translation re-initiation of mRNA (). ATF4 is a master transcription factor during the stress response. For the mRNA expression of ATF4, the binding of DDX3X to the eIF4F complex is required (). The downregulated expression level of ATF4 detected in this study thus can be involved in the possible suppression of the eIF4F complex formation under the heat stress. eIF4F (eukaryotic initiation factor 4F) is a translation initiation factor, and consists of eIF4E, eIF4G, and eIF4A. eIF4G3 is an isoform of eIF4G, and its expression level was downregulated in B. otai tadpoles. eIF4ENIF1, a transporter of eIF4E, also exhibited a downregulated expression level in response to the heat stress. It was demonstrated that eIF4F formation is disrupted, and subsequently, mRNA translation is inhibited in mice under heat stress so that protein synthesis and subsequent possible misfolding can be prevented (). The other translation regulation subgroup contained CHD9 and NCOR2, which were closely located to the KMT2 family in the PPI network, are known as transcriptional regulators.
In the methyltransferase subgroup, three proteins (KMT2A–C) of the lysine methyltransferase 2 (KMT2) family were included. This family mainly functions as a histone methyltransferase and is thus known as a transcriptional regulator. In response to the heat stress, KMT2A and KMT2B were upregulated whereas KMT2C was downregulated.
Ubiquitination subgroups including FBXW4, UBR4, RNF126, and HERC4 formed an independent PPI network. They are ubiquitin-related proteins and known to promote ubiquitination and degradation of proteins (; ). All of these proteins are E3 ubiquitin-protein ligases, however, expression patterns were different: FBXW4 and UBR4 were upregulated and HERC4 and RNF126 were downregulated in response to the heat stress.
Enrichment Analyses
Enrichment analyses performed on upregulated proteins demonstrated that three Pfam, four InterPro, and four SMART categories were enriched, which are involved in methylation and other expression regulations, and no GO term and KEGG pathway were enriched (Supplementary Table 4). Methyltransferase genes, such as KMT2A/B, were involved in the enriched categories. In the downregulated protein analyses, only two UniProt categories (coiled coil and kelch repeat), which were involved in expression regulation and ubiquitination, were enriched (Supplementary Table 4). Three proteins, ENC1.2, KLHL22, and KLHL24, and 32 proteins (e.g., EIF4G3, KMT2C, and HSPA1B), were involved in the coiled-coil and kelch repeat, respectively. There was no overlapping category between the upregulated and downregulated protein analyses. The enrichment analyses performed on all differentially expressed proteins demonstrated 16 overrepresented pathway, keywords, domains, and features while no enriched GO term was detected (Table 1). All of them, except for the coiled-coil and kelch repeat, can be directly or indirectly linked to methyltransferase activity.
TABLE 1
| Framework | ID | Description | Strength | q-value |
| KEGG pathways | xtr00310 | Lysine degradation | 1.04 | 0.0339 |
| Annotated keywords (UniProt) | KW-0489 | Methyltransferase | 0.90 | 0.0233 |
| KW-0949 | S-adenosyl-L-methionine | 0.84 | 0.0304 | |
| KW-0175 | Coiled coil | 0.24 | 0.0048 | |
| Protein Domains (Pfam) | PF05964 | F/Y-rich N-terminus | 1.93 | 0.0041 |
| PF05965 | F/Y rich C-terminus | 1.93 | 0.0041 | |
| PF00856 | SET domain | 1.27 | 0.0075 | |
| Protein domains and features (InterPro) | IPR003888 | FY-rich, N-terminal | 1.93 | 0.0070 |
| IPR003889 | FY-rich, C-terminal | 1.93 | 0.0070 | |
| IPR003616 | Post-SET domain | 1.49 | 0.0269 | |
| IPR034732 | Extended PHD (ePHD) domain | 1.31 | 0.0463 | |
| IPR001214 | SET domain | 1.16 | 0.0269 | |
| Protein domains (SMART) | SM00542 | FY-rich domain, C-terminal region | 2.03 | 0.0010 |
| SM00541 | FY-rich domain, N-terminal region | 1.93 | 0.0010 | |
| SM00508 | Cysteine-rich motif following a subset of SET domains | 1.49 | 0.0048 | |
| SM00317 | SET [Su(var)3–9, Enhancer-of-zeste, Trithorax] domain | 1.25 | 0.0034 |
Results of STRING functional enrichment analyses based on combined dataset of upregulated and downregulated proteins.
Strength represents the degree of enrichment, and the q-value is the FDR-corrected p-value.
Discussion
Through a comparative transcriptomic analysis, we identified hundreds of transcripts differentially expressed between heat-treated and control tadpoles of B. otai. Knowledge-based analyses imply that our transcriptome analyses detected sets of genes with important functions in coping with the heat stress. Specifically, the PPI network analysis identified hub proteins with functions such as HSPs, stress granule assembly, methylation, and ubiquitination.
Transcriptional and Translational Regulation Mechanisms
First, the transcriptional response to the environmental stimuli is led by epigenetic regulation mechanisms. This study of comparative transcriptional analyses identified genes, pathways, domains, and features involved in different epigenetic mechanisms like DNA and histone methylation. Methyltransferase itself does not have a specific target of genes but is involved in the global change of DNA and histone methylation status. Subsequently, it can cause the transcriptional change of genes. It is not clear how many differentially expressed genes identified in this study are regulated by the methyltransferase activity. However, expression levels of all genes mentioned in the results section, except for DDX3X on the X chromosome, are associated with DNA methylation levels of neighboring DNA methylation sites in humans ().
Expression levels of KMT2A and KMT2B were upregulated in response to heat stress, whereas that of KMT2C was downregulated. The contrasting expression patterns may be explained by the structural and functional differences between them. For example, the KMT2 enzymes can be distinguished by the number of PHD fingers: KMT2A/B contain four PHD fingers and have a similar overall domain structure. KMT2C, on the other hand, contains eight PHD fingers, which is the largest number among the known KMT2 enzymes (). These structural differences will lead to functional differences. It is known that KMT2A/B and KMT2C have unique subunits, which are known to be associated with differences of binding regions: KMT2A/B are suggested to bind to the promoter and enhancer regions, while KMT2C binds to the enhancer region (). The CXXC domain, unique to KMT2A/B, also contributes to the binding of the enzymes to CpG sites (). Taken together, the contrasting expression pattern of KMT2 genes implies that specific regulatory mechanisms, such as DNA methylation, may be responsive to heat stress.
In addition to methyltransferase activity, we identified differentially expressed genes involved in the activation and suppression of transcription and translation. NCOR2 is a transcriptional co-suppressor, and it was upregulated under heat stress. Conversely, genes involved in forming the eIF4F complex, a translation initiation factor, were downregulated. These results imply that gene transcription and translation are globally suppressed in B. otai tadpoles under heat stress. Presumably, the mechanisms contribute to the prevention of misfolding protein synthesis under heat stress. The upregulation of DDX3X, a stress granule component, is representative of the hypothesis. Stress granule formation is a well-known cellular stress response and prevents RNA from being translated.
CHD9 upregulation is also notable. The gene belongs to the chromodomain helicase DNA-binding (CHD) family. CHD is a chromatin remodeler that contributes to chromatin structural change and thus, affects chromatin accessibility. Different from the above-mentioned proteins that can change the global DNA transcription and RNA translation pattern, CHD9 is suggested to have a specific target of genes related to osteogenesis, such as RUNX2 (). Indeed, it is well documented that heat stimulus promotes osteogenesis (; ), and it is conceivable that CHD9 participates in the process. However, it is unclear whether the CHD9-mediated osteogenesis can be linked to the protective function against heat stress. Since water temperature is known to affect the timing of metamorphosis (; ), the upregulation of CHD9 possibly be related to the body plan change associated with accelerated metamorphosis.
Protein Degradation
We also identified ubiquitination-related proteins that formed an independent PPI network. Proteins damaged by denaturation, such as heat stress are degraded by proteasomes, and the process is usually mediated by ubiquitination (). Transcriptional upregulation of ubiquitination-related genes was also identified in two comparative transcriptome studies in which amphibians were exposed to cold stress (; ). Thus, in combination with above-mentioned protein synthesis suppression, ubiquitination-mediated protein degradation may contribute to reducing dysfunctional proteins.
There were four differentially expressed genes of E3 ubiquitin-protein ligases, although RNF126 and HRC4 were downregulated, while FBXW4 and UBR4 were upregulated in response to the heat stress. There are hundreds of known E3 ligases and each ligase is thought to have specific substrate, however, the relationship is still not well understood, and it is also not known how ligases recognize misfolded proteins. In addition, only a limited number of E3 ligases are known to participate in the misfolded protein ubiquitination (), and various other functions have been proposed for ubiquitination, such as DNA repair, signal transduction, and transcriptional regulation (). Therefore, we cannot provide a clear explanation for the up- and downregulations of each E3 ligase gene in response to the heat stress. However, the contrasting expression patterns strongly suggest that these genes play different roles under heat stress.
Heat Shock Proteins
The HSP family is the most famous protein family produced under heat stress. Three genes belonging to HSP40 and HSP70 were differentially expressed in B. otai tadpoles but were all downregulated in heat-treated individuals. Since HSP40 is a co-chaperone of HSP70, it is probable that they demonstrated the same transcriptional change pattern. However, we have no clear insights to explain why they were downregulated in response to heat stress. Contrary to the current study, several frog species exposed to cold stress exhibited upregulation of HSP70 (; ). Although the upregulation of HSP70 in those species is expected to function in the degradation of proteins misfolded due to the stress, it cannot explain the downregulation of HSP70 under the heat stress, which also leads to protein misfolding. The downregulation of HSP genes in response to heat stress (21°C compared to 11°C) was also observed in a salamander species (). However, upregulation of the other HSP genes were simultaneously observed and the responses of these specific genes have not been discussed in the previous study. HSPs have various functions, several of which change the responsiveness according to the duration of the stress (). Furthermore, they can be controlled by negative feedback. These gene-specific functions and regulations may underlie the suppression of HSP genes in response to heat stress.
Suppressive Transcriptional Response Against Heat Stress
It has been demonstrated that tadpoles prefer cooler water (; ) and the water temperature in hot springs can drop at night due to the lower air temperature; thus, tadpoles can avoid long-term exposure to heat even in hot spring environments. Rather, it probably important for them to tolerate the short-term heat spike. The present study showed that the majority of hub proteins identified in this study had functions related to suppression of protein synthesis and subsequent prevention of protein misfolding. These responses should negatively impact tadpole development if the suppression persists. Therefore, although not verified, the mechanisms are not likely to be sustained under long-term heat stress. That is, Buergeria tadpoles may adopt a combination of suppressive transcriptional responses to tolerate occasional heat stress and avoidance behavioral response to prevent long-term heat exposure.
Suppression of protein synthesis under heat stress is not a mechanism specific to hot spring tadpoles but a widely conserved mechanism. Thus, the mechanisms identified in this study were not evolutionarily acquired in Buergeria lineages. Many anuran species lay eggs in temporally formed shallow water pools, such as paddy fields, whereby water temperatures can reach 36°C (). Therefore, although the present study employed rapid and acute heat stress and did not take into account acclimation effects that may occur in the field, a wide range of anuran species may have a similar transcriptomic response to heat stress as that observed in B. otai tadpoles.
Limitations
It should be noted that the sample size in this study is small. Therefore, it is possible that the set of differentially expressed genes may change when the sample is increased. To reduce the uncertainty arising from the small sample size, we mainly focused on the summarized characteristics of upregulated and downregulated gene sets obtained from enrichment and PPI network analyses.
It is also important to mention that this study only focused on annotated transcripts. This means that even if there are novel heat stress response-related transcripts specific to Buergeria frogs, they cannot be identified in this study design. Since RNA-seq data obtained in this study is publicly available, it will possible to approach this point in future studies.
In addition, RNA was extracted from the whole body. Thus, genes showing tissue- or cell type-specific heat stress response patterns are unlikely to be detected in this study. The differentially expressed genes identified in this study are genes whose expression patterns do not conflict among tissues or cell types.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Statements
Data availability statement
RNA-seq data obtained in this study are deposited in DDBJ Sequence Read Archive (DRA) (BioProject accession number: PRJDB12187). Codes used in this study are available on GitHub (https://github.com/RolyPolyCoily/BuergeriaOtai_2021).
Ethics statement
Ethical review and approval was not required for the animal study because amphibian tadpoles were not subject to ethical review in Hiroshima University where the experiment was conducted at that time; therefore, no institutional ethical approval was obtained.
Author contributions
SK designed the study, carried out field and laboratory works and data analyses, and led the writing of the draft. MM carried out data analyses and participated in the writing of the draft. J-WL, K-HL, Y-PL, and YL participated in sample collection and critically revised the manuscript. S-ML coordinated the field work, participated design of the study, and critically revised the manuscript. TI designed the study, carried out laboratory works and data analyses, and revised the manuscript. All authors gave final approval for publication.
Funding
The study was supported by the JSPS KAKENHI (18K14795 and 25-5065 for SK and 23710282 and 18K06365 for TI).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.706887/full#supplementary-material
Supplementary Table 1Summary of RNA-seq and subsequent analyses.
Supplementary Table 2List of transcripts differentially expressed between heat-treated and control tadpoles with analytical results and annotations.
Supplementary Table 3Scores of degree centrality and betweenness centrality for each protein mapped on STRING protein-protein interaction networks.
Supplementary Table 4Results of STRING enrichment analyses.
Supplementary Text 1Methods and results of additional annotations.
Footnotes
References
1
AdjibadeP.Grenier St-SauveurV.BergemanJ.HuotM. E.KhandjianE. W.MazrouiR. (2017). DDX3 regulates endoplasmic reticulum stress-induced ATF4 expression. Sci. Rep.7, 1–12. 10.1038/s41598-017-14262-7
2
AliM.HomR. A.BlakesleeW.IkenouyeL.KutateladzeT. G. (2014). Diverse functions of PHD fingers of the MLL/KMT2 subfamily.Biochim. Biophys. Acta Mol. Cell Res.1843366–371. 10.1016/j.bbamcr.2013.11.016
3
AltschulS. F.GishW.MillerW.MyersE. W.LipmanD. J. (1990). Basic local alignment search tool.J. Mol. Biol.215403–410. 10.1016/S0022-2836(05)80360-2
4
BoonstraR. (2013). The ecology of stress: a marriage of disciplines.Funct. Ecol.277–10. 10.1111/1365-2435.12048
5
BrayN. L.PimentelH.MelstedP.PachterL. (2016). Near-optimal probabilistic RNA-seq quantification.Nat. Biotechnol.34525–527. 10.1038/nbt.3519
6
ChenJ.ShiZ.-D.JiX.MoralesJ.ZhangJ.KaurN.et al (2013). Enhanced osteogenesis of human mesenchymal stem cells by periodic heat shock in self-assembling peptide hydrogel.Tissue Eng. Part A19716–728. 10.1089/ten.tea.2012.0070
7
ChenK.-L.WangH.-L.JiangL.-Z.QianY.YangC.-X.ChangW.-W.et al (2020). Heat stress induces apoptosis through disruption of dynamic mitochondrial networks in dairy cow mammary epithelial cells.In Vitro Cell Dev. Biol. Anim.56322–331. 10.1007/s11626-020-00446-5
8
ChenS.-J.LiaoD.-L.ChenC.-H.WangT.-Y.ChenK.-C. (2019). Construction and analysis of protein-protein interaction network of heroin use disorder.Sci. Rep.9:4980. 10.1038/s41598-019-41552-z
9
ChenT.-C.KamY.-C.LinY.-S. (2001). Thermal Physiology and reproductive phenology of Buergeria japonica (Rhacophoridae) breeding in a stream and a geothermal hotspring in Taiwan.Zoolog. Sci.18591–596. 10.2108/zsj.18.591
10
ChhanganiD.JanaN. R.MishraA. (2013). Misfolded proteins recognition strategies of E3 ubiquitin ligases and neurodegenerative diseases.Mol. Neurobiol.47302–312. 10.1007/s12035-012-8351-0
11
ClayT. A.SteffenM. A.TregliaM. L.TorresC. D.Trujano-AlvarezA. L.BonettR. M. (2019). Multiple stressors produce differential transcriptomic patterns in a stream-dwelling salamander.BMC Genomics20:482. 10.1186/s12864-019-5814-y
12
CrespiE. J.WilliamsT. D.JessopT. S.DelehantyB. (2013). Life history and the ecology of stress: how do glucocorticoid hormones influence life-history variation in animals?Funct. Ecol.2793–106. 10.1111/1365-2435.12009
13
CsardiG.NepuszT. (2006). The igraph software package for complex network research.InterJournal Complex Syst.16951–9.
14
CuestaR.LaroiaG.SchneiderR. J. (2000). Chaperone hsp27 inhibits translation during heat shock by binding eIF4G and facilitating dissociation of cap-initiation complexes.Genes Dev.141460–1470. 10.1101/gad.14.12.1460
15
do AmaralM. C. F.FrisbieJ.CrumR. J.GoldsteinD. L.KraneC. M. (2020). Hepatic transcriptome of the freeze-tolerant Cope’s gray treefrog, Dryophytes chrysoscelis: responses to cold acclimation and freezing.BMC Genomics21:226. 10.1186/s12864-020-6602-4
16
GoldsteinJ. A.HoffK.vonS.HillyardS. D. (2017). The effect of temperature on development and behaviour of relict leopard frog tadpoles.Conserv. Physiol.5:cow075. 10.1093/conphys/cow075
17
GormallyB. M. G.RomeroL. M. (2020). What are you actually measuring? A review of techniques that integrate the stress response on distinct time-scales.Funct. Ecol.342030–2044. 10.1111/1365-2435.13648
18
GosnerK. L. (1960). A simplified table for staging anuran embryos larvae with notes on identification.Herpetologica16183–190.
19
GrabherrM. G.HaasB. J.YassourM.LevinJ. Z.ThompsonD. A.AmitI.et al (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome.Nat. Biotechnol.29644–652. 10.1038/nbt.1883
20
HaasB. J.PapanicolaouA.YassourM.GrabherrM.BloodP. D.BowdenJ.et al (2013). De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis.Nat. Protoc.81494–1512. 10.1038/nprot.2013.084
21
HaramuraT. (2004). Salinity and other abiotic characteristics of oviposition sites of the Rhacophorid Frog, Buergeria japonica, in coastal habitat.Curr. Herpetol.2381–84. 10.5358/hsj.23.81
22
HershkoA. (1996). The ubiquitin system for protein degradation.Biochem. Soc. Trans.24:627S. 10.1042/bst024627sc
23
HillikerA.GaoZ.JankowskyE.ParkerR. (2011). The DEAD-box protein Ded1 modulates translation by the formation and resolution of an eIF4F-mRNA complex.Mol. Cell43962–972. 10.1016/j.molcel.2011.08.008
24
IkutaK.UrakawaH.KozawaE.HamadaS.OtaT.KatoR.et al (2015). In vivo heat-stimulus-triggered osteogenesis.Int. J. Hyperth.3158–66. 10.3109/02656736.2014.988662
25
KissA. J.MuirT. J.LeeR. E.Jr.CostanzoJ. P. (2011). Seasonal variation in the hepatoproteome of the dehydration- and freeze-tolerant wood frog, Rana sylvatica.Int. J. Mol. Sci.128406–8414. 10.3390/ijms12128406
26
KomakiS.LauQ.IgawaT. (2016b). Living in a Japanese onsen: field observations and physiological measurements of hot spring amphibian tadpoles, Buergeria japonica.Amphibia Reptilia37311–314. 10.1163/15685381-00003052
27
KomakiS.IgawaT.LinS.-M.SumidaM. (2016a). Salinity and thermal tolerance of Japanese stream tree frog (Buergeria japonica) tadpoles from island populations.Herpetol. J.26209–213.
28
KomakiS.ShiwaY.FurukawaR.HachiyaT.OhmomoH.OtomoR.et al (2018). iMETHYL: an integrative database of human DNA methylation, gene expression, and genomic variation.Hum. Genome Var.5:18008. 10.1038/hgv.2018.8
29
KomakiS.SutohY.KobayashiK.SaitoS.SaitoC. T.IgawaT.et al (2020). Hot spring frogs (Buergeria japonica) prefer cooler water to hot water.Ecol. Evol.109466–9473. 10.1002/ece3.6637
30
LakeD.CorrêaS. A. L.MüllerJ. (2016). Negative feedback regulation of the ERK1/2 MAPK pathway.Cell. Mol. Life Sci.734397–4413. 10.1007/s00018-016-2297-8
31
MacielT. A.JuncáF. A. (2009). Effects of temperature and volume of water on the growth and development of tadpoles of Pleurodema diplolister and Rhinella granulosa (Amphibia: Anura).Zoologica26413–418. 10.1590/S1984-46702009000300005
32
MaoY.FisherD. W.YangS.KeszyckiR. M.DongH. (2020). Protein-protein interactions underlying the behavioral and psychological symptoms of dementia (BPSD) and Alzheimer’s disease.PLoS One15:e0226021. 10.1371/journal.pone.0226021
33
MaruyamaA.NemotoM.HamasakiT.IshidaS.KuwagataT. (2017). A water temperature simulation model for rice paddies with variable water depths.Water Resour. Res.5310065–10084. 10.1002/2017WR021019
34
MastrangeloA. M.MaroneD.LaidòG.De LeonardisA. M.De VitaP. (2012). Alternative splicing: enhancing ability to cope with stress via transcriptome plasticity.Plant Sci.185–18640–49. 10.1016/j.plantsci.2011.09.006
35
MatsuiM.MaedaN. (2018). Encycropedia of Japanese Frogs.Tokyo: Bun-ichi Sogo Shuppan.
36
MoritaT. (1988). Terrestrial vertebrate fauna on Kuchinoshima Island of the Tokara Archipelago where wild cows live [transrated from Japanese title].Rikabukaishi3017–19.
37
NagasawaK.TanizakiY.OkuiT.WataraiA.UedaS.KatoT. (2013). Significant modulation of the hepatic proteome induced by exposure to low temperature in Xenopus laevis.Biol. Open21057–1069. 10.1242/bio.20136106
38
NewtonA. H.PaskA. J. (2020). CHD9 upregulates RUNX2 and has a potential role in skeletal evolution.BMC Mol. Cell Biol.21:27. 10.1186/s12860-020-00270-5
39
O’LearyN. A.WrightM. W.BristerJ. R.CiufoS.HaddadD.McVeighR.et al (2016). Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation.Nucleic Acids Res.44D733–D745. 10.1093/nar/gkv1189
40
RaoR. C.DouY. (2015). Hijacked in cancer: the KMT2 (MLL) family of methyltransferases.Nat. Rev. Cancer15334–346. 10.1038/nrc3929
41
RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.Bioinformatics26139–140. 10.1093/bioinformatics/btp616
42
SeppeyM.ManniM.ZdobnovE. M. (2019). BUSCO: assessing genome assembly and annotation completeness.Methods Mol. Biol.1962227–245. 10.1007/978-1-4939-9173-0_14
43
ShenW.LeS.LiY.HuF. (2016). SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation.PLoS One11:e0163962. 10.1371/journal.pone.0163962
44
SonesonC.LoveM. I.RobinsonM. D. (2016). Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences.F1000Research4:1521. 10.12688/f1000research.7563.2
45
SzklarczykD.GableA. L.LyonD.JungeA.WyderS.Huerta-CepasJ.et al (2019). STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.Nucleic Acids Res.47D607–D613. 10.1093/nar/gky1131
46
UniProt Consortium (2019). UniProt: a worldwide hub of protein knowledge.Nucleic Acids Res.47D506–D515. 10.1093/nar/gky1049
47
WilkinsonK. D.VentiiK. H.FriedrichK. L.MullallyJ. E. (2005). The ubiquitin signal: assembly, recognition and termination.EMBO Rep.6815–820. 10.1038/sj.embor.7400506
48
WingfieldJ. C. (2013). Ecological processes and the ecology of stress: the impacts of abiotic environmental factors.Funct. Ecol.2737–44. 10.1111/1365-2435.12039
49
WuC.-S.KamY.-C. (2005). Thermal tolerance and thermoregulation by Taiwanese rhacophorid tadpoles (Buergeria japonica) living in geothermal hot springs and streams.Herpetologica6135–46. 10.1655/04-50
50
YuT.DeusterP.ChenY. (2016). Role of dynamin-related protein 1-mediated mitochondrial fission in resistance of mouse C2C12 myoblasts to heat injury.J. Physiol.5947419–7433. 10.1113/JP272885
51
YuT.FerdjallahI.ElenbergF.ChenS. K.DeusterP.ChenY. (2018). Mitochondrial fission contributes to heat-induced oxidative stress in skeletal muscle but not hyperthermia in mice.Life Sci.2006–14. 10.1016/j.lfs.2018.02.031
Summary
Keywords
amphibia, anura, differentially expressed genes, protein-protein interaction network, RNA-seq, high thermal tolerance, heat stress response
Citation
Komaki S, Matsunami M, Lin J-W, Lee K-H, Lin Y-P, Lee Y, Lin S-M and Igawa T (2021) Transcriptomic Changes in Hot Spring Frog Tadpoles (Buergeria otai) in Response to Heat Stress. Front. Ecol. Evol. 9:706887. doi: 10.3389/fevo.2021.706887
Received
08 May 2021
Accepted
08 September 2021
Published
13 October 2021
Volume
9 - 2021
Edited by
Robert A. Haney, Ball State University, United States
Reviewed by
Saurabh Sarkar, Guskara Mahavidyalaya, India; Ivan Gomez-Mestre, Estación Biológica de Doñana (EBD), Spain
Updates
Copyright
© 2021 Komaki, Matsunami, Lin, Lee, Lin, Lee, Lin and Igawa.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Takeshi Igawa, tigawa@hiroshima-u.ac.jp
This article was submitted to Evolutionary and Population Genetics, a section of the journal Frontiers in Ecology and Evolution
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.