Identification of Novel mRNA Isoforms Associated with Acute Heat Stress Response Using RNA Sequencing Data in Sprague Dawley Rats

Simple Summary Global warming and events of heat waves even in temperate climatic zones signify the importance of the genetic makeup of heat stress. Our earlier acute heat stress exposure study (from 30 min to 120 min treatments) in rats led to the designation of aberrant differentially expressed genes (DEGs) and pathways through RNA-sequencing of the three most relevant tissues, i.e., blood, liver, and adrenal glands. However, the causes and mechanisms of differential expression of genes associated with heat stress are still unclear. Using the same RNA-sequencing data, this study identified the differentially expressed mRNA isoforms and narrowed down the most reliable DEG markers and molecular pathways that underlie the thermoregulatory mechanisms. The spatial-temporal differential expression pattern of heat stress-responsive differential mRNA isoforms may be ultimately used in marker-assisted selection for improved thermotolerance. Abstract The molecular mechanisms underlying heat stress tolerance in animals to high temperatures remain unclear. This study identified the differentially expressed mRNA isoforms which narrowed down the most reliable DEG markers and molecular pathways that underlie the mechanisms of thermoregulation. This experiment was performed on Sprague Dawley rats housed at 22 °C (control group; CT), and three acute heat-stressed groups housed at 42 °C for 30 min (H30), 60 min (H60), and 120 min (H120). Earlier, we demonstrated that acute heat stress increased the rectal temperature of rats, caused abnormal changes in the blood biochemical parameters, as well as induced dramatic changes in the expression levels of genes through epigenetics and post-transcriptional regulation. Transcriptomic analysis using RNA-Sequencing (RNA-Seq) data obtained previously from blood (CT and H120), liver (CT, H30, H60, and H120), and adrenal glands (CT, H30, H60, and H120) was performed. The differentially expressed mRNA isoforms (DEIs) were identified and annotated by the CLC Genomics Workbench. Biological process and metabolic pathway analyses were performed using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. A total of 225, 5764, and 4988 DEIs in the blood, liver, and adrenal glands were observed. Furthermore, the number of novel differentially expressed transcript lengths with annotated genes and novel differentially expressed transcript with non-annotated genes were 136 and 8 in blood, 3549 and 120 in the liver, as well as 3078 and 220 in adrenal glands, respectively. About 35 genes were involved in the heat stress response, out of which, Dnaja1, LOC680121, Chordc1, AABR07011951.1, Hsp90aa1, Hspa1b, Cdkn1a, Hmox1, Bag3, and Dnaja4 were commonly identified in the liver and adrenal glands, suggesting that these genes may regulate heat stress response through interactions between the liver and adrenal glands. In conclusion, this study would enhance our understanding of the complex underlying mechanisms of acute heat stress, and the identified mRNA isoforms and genes can be used as potential candidates for thermotolerance selection in mammals.


Introduction
As the global temperature increases, heat stress has become a big challenge for human health and animal survival. An assessment report released by the Intergovernmental Panel on Climate Change (IPCC) predicts that by 2100, the global average surface temperature will rise by 0.3-4.8 • C [1], suggesting that it is more urgent to explore the body's adaptation mechanism to heat stress and develop strategies to resist the adverse effects of heat stress. Heat stress is a major stressor affecting the health of animals, and body homeostasis is disturbed while employing thermoregulation [2]. The heat stress response is carried out through the hypothalamic-pituitary-adrenal (HPA) axis, causing hormonal changes including an increase in adrenocorticotropic hormones [3,4]. Homeotherms, upon heat stress, accelerate the respiratory rate, and sweating ensues, among other physiological responses [5,6]. Additional changes induced by heat stress cause spatiotemporal temperature distribution in various tissues [7]. Moreover, being a public health concern [8], heat stress causes physiological and biochemical changes in the body of rodents [9] and livestock [10][11][12]. These changes end up in the morbidity and decline of production and reproduction [10,13]. The negative effects of heat stress are characterized by high oxidative stress [14], molecular and transcriptional changes [15,16], along with posttranscriptional [17,18] and epigenetic changes [19] at the cellular and tissue levels. However, there is limited information about the molecular mechanisms of heat stress response in mammals.
Heat stress is a complex regulatory process combining the neurohormonal, oxidative, and immune responses [20][21][22]. Next-generation sequencing technologies, such as RNA sequencing (RNA-Seq), allow the large-scale, rapid, and accurate identification of heat stress-responsive genes, contributing to elucidating the molecular mechanisms of heat stress more deeply [23][24][25][26][27]. Moreover, numerous studies have demonstrated that the activity of some heat stress-responsive genes is not only regulated by their post-translational modifications [28][29][30] but also affected by changes in mRNA expression. Alternative splicing profoundly plays a vital role in many physiological processes. A study has shown that 60% of disease-causing point mutations are induced by interference with normal mRNA splicing during transcription [31]. In addition, alternative splicing was found to be involved in the process of biotic and abiotic stress responses [32][33][34]. As early as 1994, Takechi et al. reported that alternatively spliced mRNA with 169 additional nucleotides in the 5 noncoding region was found after heat shock in comparison with mRNA transcribed under non-heat shock conditions [35]. Ju et al. found that in heat-stressed Bama miniature pigs, the second exon of the TLR4 gene was spliced and 167 bp shorter in the alternative splicing variant [36]. In heat-stressed Drosophila, Nobuhiro et al. identified three new isoforms of heat shock transcription factor mRNA, including HSFb, HSFc, and HSFd, and found that the ratio of HSFb increased upon heat stress [28]. Although numerous studies reveal hundreds or thousands of heat stress-related genes in animals [9,[37][38][39][40], studies focusing on large-scale mRNA isoforms in heat stressed-animals are still scarce in the literature.
In our previous study, the mRNA expression profiles of the blood, liver, and adrenal glands of rats in the control group (CT, 22 ± 1 • C, relative humidity 50%) and three heat-stressed groups fed at 42 • C for 30 min (H30), 60 min (H60) and 120 min (H120) (relative humidity 50%) were investigated using RNA-Seq [9,17,41]. In this context, the main hypothesis of the current study is that heat stress may change the number of mRNA isoforms involved in the expression levels of transcripts that are transcribed by genes involved in thermotolerance mechanisms in different tissues under different heat stress durations. Therefore, the objective of the present study was to identify the novel mRNA isoforms differentially expressed in blood, liver, and adrenal gland tissues of rats under H30, Biology 2022, 11, 1740 3 of 24 H60, and H120 conditions using RNA-Seq data obtained previously [9,17]. In addition, the mRNA isoforms identified in this study were integrated into different functional analyses in order to identify active physiological pathways and transcripts involved in the process of the heat stress response.

Animals and Sample Collection
Twenty 8-week-old female Sprague Dawley (SD) rats (Beijing Vital River Laboratory, Animal Technology Co., Ltd., Beijing, China) weighing 205 ± 7.16 g were to be used as subjects. During the whole experiment, water and feed were provided ad libitum. Based on changes in rectal temperature and 17 biochemical indicators of blood housed at H30, H60, and H120 conditions [41], a heat-stressed rat model was established in which 48 samples including blood (n = 4 each in the CT and H120 groups), liver (n = 5 each in the CT, H30, H60, and H120 groups), and adrenal glands (n = 5 each in the CT, H30, H60, and H120 groups) were collected. The blood, liver, and adrenal gland tissues in the same treatment group were collected from the same rat. Briefly, the experimental rats were anesthetized with an injection of 1%, 1.2 mL pentobarbital sodium (0.6 mL/100 g body weight), sacrificed, and quickly dissected by sterile surgical scissors and forceps (Shinva Medical Instrument Co. Ltd., Zibo, China) [42]. About 4 mL of blood was collected using Vacutainer ® tubes (Becton Dickinson, Plymouth, UK) containing Ethylenediaminetetraacetic acid (EDTA) from each rat and immediately placed on ice. Then, the peripheral blood mononuclear cells (PBMC) were collected by RNase-free spear via centrifugation (10 min, 3000 rpm) in a 2 mL centrifuge tube containing 1 mL Trizol (Invitrogen 15596018, Thermo Fisher Scientific Inc., Waltham, MA, USA) and stored at −80 • C for RNA extraction. The liver and adrenal glands were collected, washed in ice-cold phosphate-buffered solution, and snap-frozen immediately in liquid nitrogen.

RNA Isolation and Library Construction
The RNA reagent (HUAYUEYANG Biotechnology Co. Ltd., Beijing, China) was used to extract the total RNA from 48 samples of 20 SD rats according to the manufacturer's protocol. The purity of the RNA samples was assessed at A260/230 nm and A260/280 nm ratios using the NanoDrop 2100 (Thermo Fisher Scientific Inc., Waltham, MA, USA). The integrity of the RNA (RIN) was determined using an RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA), and the RNA concentration and genomic DNA contamination were determined using a Qubit ® 2.0 Fluorometer (Thermo Fisher Scientific Inc., Waltham, USA). The OD 260 /OD 280 ratio varied from 1.8 to 2.0 for all samples, and RIN values > 7.0, indicating good RNA quality [43].
A total of 3 µg RNA from each sample was used to prepare the cDNA libraries using the NEBNext ® UltraTM Directional RNA Library Prep Kit from Illumina ® (NEB, San Diego, CA, USA), following the manufacturer's recommendation. Paired-end (2 × 150 pb) reads of samples in the CT and H120 groups were sequenced using a HiSeq 2000 sequencer (Illumina, San Diego, CA, USA) in a previous study [9], and the reads of H30 and H60 groups were sequenced using the HiSeq 2500 sequencer (Illumina, San Diego, CA, USA) with 2 × 150 pb reads paired-end [17]. Both of the RNA-Seq protocols were performed in the Novogene Technology Co., Ltd., located in Tianjin of China.

Sequence Assembly and Quantification
All the transcriptomic analyses of the reads generated in CT, H30, H60, and H120 were performed using the CLC Genomics Workbench software 12.0 (CLC Bio, Aarhus, Denmark). Quality control analyses were performed following the parameters described in Cánovas et al. [44]. The paired-end sequenced reads were mapped to the annotated reference genome Rattus norvegicus 6.0 (rn6, ftp://ftp.ensembl.org/pub/release-105/ genbank/rattus_norvegicus/, accessed on 11 February 2022) using the "Large Gap Read Mapping" tool. Briefly, this tool allows to map the RNA-Seq reads that span introns without requiring prior transcript annotations, which helps to conduct the transcript discovery analysis [45].
The transcript discovery analysis was performed using the "Transcript Discovery" tool implemented in the CLC Genomics Workbench environment (CLC Bio, Aarhus, Denmark) [45,46], which takes large gap read mapping tracks, uses gene and transcript annotations as an input, and produces optimized gene and transcript track annotations as outputs. Therefore, for each gene region, there was a set of transcript annotations that can explain the observed exons and splice sites in this region (data files available but not attached here to this manuscript) [46,47]. Then, the sequences of the samples were mapped to the created tracks (genes and transcripts tracks) using the rat reference genome as a map via the "RNA-Sequencing analysis tool" implemented in the CLC Genomics Workbench. The reads per kilobase per million mapped reads (RPKM) were used to quantify the expression levels for different transcripts across RNA-Seq libraries and were transformed by log10 [48,49]. Transcripts with RPKM ≥ 0.2 were defined as expressed. The expression levels of transcripts identified in each sample were classified as lowly expressed, moderately expressed, and highly expressed according to the criteria of RPKM ≤ 50, 50 < RPKM < 500, and RPKM ≥ 500 [44]. Principal component analysis (PCA, an unsupervised machine learning technique that seeks to find principal components) was performed on the gene expression by using the prcomp R function. The pair-wise Pearson correlation coefficient (PCC, r) between any two of the four or five biological replicates within each treatment of each comparison was also calculated by correlation procedure in the R package (version 4.2.0, https://cran.rproject.org/src/base/R-3/, accessed on 22 April 2022), where the Pearson's correlation significance was computed using the Hmisc procedure.

Identification of Differential mRNA Isoform Expression
The CLC Genomics Workbench software 12.0 was used to conduct the differential mRNA isoform expression analyses [45,50]. Therefore, a two-stage experiment was carried out, applying empirical statistical analysis that implements an "Exact Test" developed by [51] and incorporated into the edgeR [52]. The statistical test was performed for each set up experiment, using the original count values and two parameters for the estimation of the dispersion: 1-"total count filter cut-off > 5, which specifies which features should be considered when estimating the common dispersion component". 2-estimate tag-wise dispersions, which allow a weighted combination of the tag-wise and common dispersion for each transcript [45]. In addition, the original expression values were transformed based on a logarithm (log10) and normalized using scaling normalization described by [53] to ensure that the samples are comparable and assumptions on the data for analysis are met. The isoforms in the blood that met the thresholds of p < 0.05, |fold change (FC)| ≥ 2, and in the liver and adrenal glands met with p ≤ 0.001, FDR ≤ 0.05, and |FC| ≥ 2 were identified as differentially expressed (DEIs) [54], respectively.

Functional Enrichment Analysis and Gene Annotation
The main category biological process (BP) was performed using gene ontology (GO) enrichment analysis [55]. The analysis was performed using the list of genes associated with the DEIs in each tissue under various heat stress conditions. Functional evidence of the relationship between the significant GO terms (FDR ≤ 0.05) and the target phenotypes (e.g., response to heat stress) was identified. The metabolic pathway analysis was performed by the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg, accessed on 30 April 2022) database. The non-redundant biological terms for large clusters of genes in functionally related groups network were analyzed by ClueGO [56] and visualized by Cytoscape 3.8.2 [57].
To annotate novel mRNA isoforms of unknown genes, their coding sequence was obtained from the Genome Data Viewer tool, which is implemented in the National Center of Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/genome/ gdv/browser/genome/, accessed on 6 May 2022). Thereafter, the 'Nucleotide BLAST' method available on the Basic Local Alignment Search Tool (BLAST; https://blast.ncbi. nlm.nih.gov/Blast.cgi, accessed on 6 May 2022) was used to detect similarities between the nucleotide sequences from the novel mRNA isoforms and the nucleotide (nt) collection available in the database. The nt collection consists of GenBank, EMBL, DDBJ, PDB, and RefSeq sequence data, but it excludes EST, STS, GSS, WGS, TSA, patent sequences, phase 0, 1, and 2 HTGS sequences, and sequences longer than 100 Mb. The database is nonredundant; i.e., identical sequences are merged into one entry (but it keeps the accession, GI, title, and taxonomy information from each entry).
Draw Venn Diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/, accessed on 10 May 2022) was used to analyze the common and specific genes identified under various heat stress conditions and in different tissues. Finally, we searched our previous list of differentially expressed genes [41] for the genes associated with the differentially expressed known transcripts and the novel transcript lengths, defined them as DEIDEGs, and summarized the common DEIDEGs involved in the process of thermoregulation.

Global Landscape of the Rat Transcriptome
A total of 45 billion reads of 150 bp paired-end RNA-Seq data were produced, corresponding to an average of~66 million reads per sample (Table S1). High-quality clean reads were mapped to the rn6 genome. On average, 54,042, 86,680, and 102,927 transcripts were identified in the blood, liver, and adrenal glands, mapping to 34,239, 37,043, and 38,689 genes, respectively. The PCA results between any treatment samples within each tissue showed that the samples were divided mainly by various treatments ( Figure S1A). Then, the pair-wise PCC analysis yielded 6 or 10 pair-wise r values per sample group. For example, the mean r value in the CT of the blood was 0.9978 (ranging from 0.9971-0.9996), and the mean r value in the H120 group was 0.9995 (ranging from 0.9989-0.9998) ( Figure S1B). The results of all the PCCs indicated a high level of measurement consistency among the biological replicates ( Figure S1B Statistics of the number of transcripts identified in blood, liver, and adrenal glands under various heat stress conditions. High, medium, and low mean the high, moderate, and low expressed transcripts. The statistic numbers marked in white, orange, and red means the ratio of transcripts with different expression levels to all transcripts identified in blood, liver, and adrenal glands under various heat stress conditions. B, L, and A mean blood, liver, and adrenal gland tissues, respectively. CT means rats housed at 22 ± 1 °C and relative 50% humidity; the H30, H60, and, H120 mean the rats exposed to 42 °C for 30 min, 60 min, and 120 min with a relative humidity of 50%. different expression levels to all transcripts identified in blood, liver, and adrenal glands under various heat stress conditions. B, L, and A mean blood, liver, and adrenal gland tissues, respectively. CT means rats housed at 22 ± 1 • C and relative 50% humidity; the H30, H60, and, H120 mean the rats exposed to 42 • C for 30 min, 60 min, and 120 min with a relative humidity of 50%.

Identification of Various Types of the Differentially Expressed mRNA Isoforms
An overview of the DEIs between the treatments in each tissue is shown in Table 1. A total of 225,2072,2086,1606,1130,1994, and 1864 DEIs were detected in the CT vs. H120 of blood, CT vs. H30, CT vs. H60, and CT vs. H120 comparisons of the liver and adrenal glands, respectively. A larger number of DEIs was identified in the liver compared with the adrenal glands and blood. Over all of the samples, the total number of novel differentially expressed transcript lengths of the annotated genes (6763) in all tissues was 1.75 and 19.43 times higher than the known DEIs (3866) and novel differentially expressed transcripts of non-annotated genes (348). Note: CT means control group, rats housed at 22 ± 1 • C and relative humidity of 50%; H30, H60, and H120 mean 42 • C heat stress for 30 min (H30), 60 min (H60), and 120 min (H120) with a relative humidity of 50%; DEIs means differentially expressed mRNA isoforms; Up means the expression of the mRNA isoforms was up-regulated in the heat treatment group compared to the CT; Down means the expression of the mRNA isoforms was down-regulated in the heat treatment group compared to the CT. B, L, and A refer to blood, liver, and adrenal gland tissues.
3.2.1. Differentially Expressed mRNA Isoforms Annotated in the Rat Genome (Rattus Norvegicus 6.0) A total of 81 known mRNA isoforms were differentially expressed in the blood in the CT vs. H120 comparison; 70.37% of them were up-regulated, and 29.63% were downregulated. The Rpl11-201 (ribosomal protein L11) was the top one annotated DEI with the highest FC of 206.84 (p = 8.34 × 10 5 ), followed by Rpl36al-201 (ribosomal protein L36alike) (p = 1.75 × 10 4 ) (Table S2). In proliferating cells, the ribosomal stress response can be triggered by both increased and decreased ribosomal biogenesis, accompanied by the activation of the p53 pathway [58]. Rpl11 may suppress cell growth via p53-dependent and/or independent mechanisms by interacting with other proteins and RNAs [59]. There is no research related to the role of Rpl11 in heat stress. In addition, previous studies have reported that the transcription-independent p53 can trigger heat stress-induced apoptosis [60]. Therefore, it is necessary to further study whether Rpl11 regulates cell activity during heat stress.
In the liver (Table S3), an average of 698 known DEIs were detected in each comparison, and about 69.12% (482) of them were significantly down-regulated (p ≤ 0.001, FDR ≤ 0.05 and |FC| ≥ 2). Among all the comparisons of the liver, the most significantly changed known DEI was Dbi-202 (diazepam-binding inhibitor, acyl-CoA-binding protein), with FCs of −6671.54 (FDR = 8.37 × 10 87 ), −4386.08 (FDR = 1.54 × 10 85 ), and −4386.08 (FDR = 1.54 × 10 85 ) in CT vs. H30, CT vs. H60, and CT vs. H120 comparisons, respectively. The Dbi gene plays a role in the regulation of mitochondrial steroidogenesis and acyl-CoA metabolism. Previous studies have indicated that the DBI-related peptide can protect neurons and astrocytes from oxidative stress-induced apoptosis through its metabotropic receptor [61]. It is well known that oxidative stress occurs as a consequence of the imbalance between antioxidant defense and reactive oxygen species (ROS) production, and heat stress can stimulate the production of ROS [20], which suggests that Dbi may play an indirect role in the heat stress response. The Venn diagram analysis of all DEIs identified in three comparisons of the liver has shown that 97 isoforms were shared among all comparisons, including 83 significantly up-regulated isoforms and 14 significantly down-regulated isoforms ( Figure 2A).
In the liver (Table S3), an average of 698 known DEIs were detected in each comparison, and about 69.12% (482) of them were significantly down-regulated (p ≤ 0.001, FDR ≤ 0.05 and |FC| ≥ 2). Among all the comparisons of the liver, the most significantly changed known DEI was Dbi-202 (diazepam-binding inhibitor, acyl-CoA-binding protein), with FCs of −6671.54 (FDR = 8.37 × 10 87 ), −4386.08 (FDR = 1.54 × 10 85 ), and −4386.08 (FDR = 1.54 × 10 85 ) in CT vs. H30, CT vs. H60, and CT vs. H120 comparisons, respectively. The Dbi gene plays a role in the regulation of mitochondrial steroidogenesis and acyl-CoA metabolism. Previous studies have indicated that the DBI-related peptide can protect neurons and astrocytes from oxidative stress-induced apoptosis through its metabotropic receptor [61]. It is well known that oxidative stress occurs as a consequence of the imbalance between antioxidant defense and reactive oxygen species (ROS) production, and heat stress can stimulate the production of ROS [20], which suggests that Dbi may play an indirect role in the heat stress response. The Venn diagram analysis of all DEIs identified in three comparisons of the liver has shown that 97 isoforms were shared among all comparisons, including 83 significantly up-regulated isoforms and 14 significantly down-regulated isoforms ( Figure 2A).  In adrenal glands, a total of 1690 isoforms were differentially expressed; 70.23%, 71.08%, and 65.83% of them were significantly up-regulated, which were 2.36, 2.46, and 1.93 times higher than down-regulated isoforms identified in CT vs. H30, CT vs. H60, and CT vs. H120 comparisons, respectively (Table 1). The Rps15a-201 (ribosomal protein S15a) was the top one isoform that was differentially expressed (FDR = 2.43 × 10 60 ) when H30 was compared to CT, with an FC of 3225.49 (Table S4). The primary structure of Rps15a was firstly investigated in 1994, and it was shown that Rps15a has 129 amino acids [62]. Previous studies reported that Rps15a is engaged in the metabolism of RNA and proteins, rRNA processing, and translation [63][64][65]. In the comparison of CT vs. H60 (Table S4), the most significantly changed isoform was Hspe1-201 (heat shock protein family E (Hsp10) member 1), which is a stress-induced mitochondrial matrix protein and molecular chaperone [66]. However, no study was conducted on the function of Hspe1 during the heat stress response. The Giot1 (gonadotropin inducible ovarian transcription factor 1) was the top one significantly changed isoform in CT vs. H120 which was up-regulated (FC = 547.48, FDR = 2.96 × 10 30 ). No shared annotated DEIs were detected among the three comparisons in the adrenal glands ( Figure 2B). most (FC = −115.88) and was significantly down-regulated at H120 (Table S2). Furthermore, a large number of novel transcript lengths annotated for known genes that were identified as differentially expressed in blood were involved in the immune response, such as RT1-Bb (RT1 class II, locus Bb), RT1-Da, and Crip1 (cysteine-rich protein 1).
In adrenal glands (Table S6), a total of 3078 novel transcript lengths annotated for known genes were differentially expressed depending on the thresholds of p ≤ 0.001, FDR ≤ 0.05, and |FC| ≥ 2, and 156 genes were commonly identified in three comparisons, including 63 genes that were up-regulated and 93 genes that were down-regulated when heat stress occurred. We also found that most shared genes were engaged in the response to hormones, the response to heat, and metabolic processes, such as Btg1 (BTG anti-proliferation factor 1), Lrp5 (LDL receptor-related protein 5), Cyp11b3 (cytochrome P450, family 11, subfamily b, polypeptide 3), Nr4a3 (nuclear receptor subfamily 4, group A, member 3), Dnaja1 (DnaJ heat shock protein family (Hsp40) member A1), Dnaja4 (DnaJ heat shock protein family (Hsp40) member A4), Hspd1 and Hsf1 (heat shock transcription factor 1), etc.

Annotation of Differentially Expressed Novel Transcript Lengths Associated with
Non-Annotated Genes in the Rat Genome Reference (Rattus Norvegicus 6.0) The differentially expressed novel transcripts associated with non-annotated genes are shown in Tables 2-4, for blood (Table 2), liver (Table 3), and adrenal glands (Table 4), respectively. In total, 6, 48, and 66 differentially expressed novel transcripts with non-annotated genes were detected in the blood, liver, and adrenal glands, respectively. Furthermore, it was observed that those transcripts presented lengths between 248-16045 bp. Only the sequences of novel mRNA isoforms that had 100% similarity to the nucleotides in the NCBI database were retained.     Among the differentially expressed novel transcripts in the blood, Gene_338_1 and Gene_786_1 were down-regulated and presented 100% similarity with Calm2 and Gypa (Table 2). Calm2 is a kind of Ca 2+ -binding protein; it can modulate cell survival, apoptosis [67], and autophagy [68] via different pathways. A study has proved that Calm2 is involved in cell stress regulation, especially endoplasmic reticulum stress [69]. Gypa is a major erythrocyte membrane sialoglycoprotein with a potential contribution to the development of diabetic complications [70]. In the liver (Table 3), 27, 15, and 6 differentially expressed novel transcripts associated with non-annotated genes were annotated in the rat genome (Rattus norvegicus 6.0), when H30, H60, and H120 were compared to the CT groups, respectively. The Gene_541, located at the genome position "3:163817141-16381790", had the greatest degree of change among the H30 (FC = −46.82), H60 (FC = −44.03), and H120 (FC = −43.04) treatment groups. A total of 35 differentially expressed novel transcripts were obtained in the H120 conditions in the adrenal glands, which is 1.46 and 5 times more than those transcripts in the H60 and H30 conditions ( Table 4). Under H120, the gene encoding the heat shock protein 86 (Hsp86), also known as Hsp90aa1, was annotated by Gene_4216. Numerous studies have revealed the Hsp90aa1 gene can be significantly upregulated in a tissue-specific and time-dependent manner after heat stress [71]. Moreover, our previous study also confirmed the important role of HSP70AA1 in dairy cattle [72]. Thus, the Hsp90aa1 plays a pivotal role in the heat stress response.

Functional Enrichment Analysis
In order to annotate the DEIs further, GO and metabolic pathway analyses were performed on all of the annotated genes associated with DEIs (Tables S7 and S8). In summary, 185 genes were obtained in the blood (CT vs. H120) and 3189 genes were annotated in the liver, including 178, 1672, and 1339 genes in the CT vs. H30, CT vs. H60, and CT vs. H120 comparisons. In addition, a total of 3876 genes were annotated in the adrenal glands when the three heat treatment groups were compared to the CT.

GO Analysis
Based on FDR ≤ 0.05, 77 BP terms were found in the blood for the CT vs. H120 comparison, and most of them were related to the immune response, such as the regulation of adaptive immune response (GO:0002819) and the regulation of T cell-mediated immunity (GO:0002709) (Table S7). Furthermore, 48 genes were significantly enriched in terms of the response to stress (GO: 0006950) (FDR = 1.91 × 10 2 ), including Fn1, Rpl11, Dnajb2, etc., which also have the functions of the positive regulation of T cell migration, the negative regulation of the apoptotic process, as well as the negative regulation of I-kappaB kinase/NF-kappaB (NF-κB) signaling. Heat stress is well documented to have a negative impact on the immune system through humoral immune responses and cell-mediated responses [73,74]. In the current research, all of the 77 significantly enriched BPs are directly or indirectly related to the immune response (Table S7), suggesting that blood regulates heat stress by activating the immune response.
In the liver, totals of 139, 693, and 646 significantly (FDR ≤ 0.05) BP terms were detected in the CT vs. H30, CT vs. H60, and CT vs. H120 comparisons (Table S7), respectively. Eighty of them were shared in all comparisons in the liver. The most significantly enriched term was the small molecule metabolic process (GO: 0044281). There were 395 BP terms significantly enriched by genes that were in common in at least two comparisons. Among them, the process of the response to stress (GO: 0006950) was significantly enriched by 339 and 279 genes in the CT vs. H60 and CT vs. H120 comparisons (Table S7); 24 and 20 genes of them were also significantly (FDR ≤ 0.05) enriched in the process of the response to heat. Venn diagram analysis indicated that 17 genes were commonly identified in the H60 and H120 groups ( Figure 3A), and they also play critical roles (p ≤ 0.05) in the response to arsenic-containing substances ( Figure 3C). In addition, the positive regulation of the coldinduced thermogenesis (GO: 0120162) term was significantly (FDR = 3.91 × 10 2 ) clustered by 14 genes that were specifically identified in the CT vs. H120 comparison ( Figure 3B), and the other functions of these 14 genes are shown in Figure 3D.
HSP90 by binding to the heat response element [76]. This study also reported that HSF-1 not only mediated the expression of HSP during the period of heat stress but also after rewarming to 37 °C. The TRP channels are often activated by a variety of stressors (e.g., heat, cold, pH). TRPV1-4 and TRPM2 are activated by heat stress, and TRPA1 and TRPM8 are activated by cold stress [77]. In the present study, HSF-1 and Trpv2 were found to be involved in the process of the response to heat-and cold-induced thermogenesis at H60 and H120, suggesting that both HSF-1 and Trpv2 played important roles in thermoregulation, but further validation experiments are needed to confirm their functions. The Trpv2 gene can be considered as a novel target for a subsequent heat stress study.  (C) The functional annotation of genes that were significantly clustered in response to heat and were shared between the CT vs. H60 and CT vs. H120 comparisons. (D) The functional annotation of genes in the CT vs. H120 comparison of the liver that were significantly clustered in the positive regulation of cold-induced thermogenesis. Only the GO-BP terms with p ≤ 0.05 are shown in colored words.
In adrenal glands, 114, 649, and 567 BP terms were significantly (FDR ≤ 0.05) detected in the CT vs. H30, CT vs. H60, and CT vs. H120 comparisons, respectively (Table S7). Among them, 353 terms were commonly detected in these three comparisons (marked in green color in Table S7). In addition, the significant terms (FDR ≤ 0.05) related to the response to heat (GO: 0009408, Figure 4A) and cold-induced thermogenesis (GO: 0106106, Figure 4B) were enriched, and six genes including Igfbp7, Dnaja1, LOC680121 (similar to heat shock protein 8), Bag3, Hsp90aa1, and Hspd1 were shared in the CT vs. H30, CT vs. H60, and CT vs. H120 comparisons ( Figure 4A). There were 19, 19, and 17 genes of the H30, H60, and H120 groups clustered in the process of cold-induced thermogenesis, and 5 genes were commonly identified among the CT vs. H30, CT vs. H60, and CT vs. H120 comparisons ( Figure 4B). Furthermore, the Hsf1 and Trpv2 genes were significantly enriched both in the response to heat-and cold-induced thermogenesis processes at H60 and H120 ( Figure 4A,B). Further functional annotation analysis was performed on genes that engaged in response to heat (n = 30) and cold-induced thermogenesis (n = 36). The results show that genes clustered in response to heat also play significant (p ≤ 0.01) roles in protein folding and the positive regulation of reactive oxygen species metabolic processes ( Figure 4C). Genes clustered in the cold-induced thermogenesis process were also significantly (p ≤ 0.01) involved in the regulation of transcription from the RNA polymerase H promoter in response to stress, sterol homeostasis, and the regulation of DNA-template transcription in response to stress ( Figure 4D).
In the liver, 31, 46, and 51 pathways were significantly enriched in the CT vs. H30 CT vs. H60, and CT vs. H120 comparisons, respectively (Table S8). Eleven pathways were commonly shared in the three comparisons, including fatty acid metabolism (ko01212) and the PPAR signaling pathway (ko03320) ( Table S8). Fatty acids can regulate immune homeostasis by the activation of fatty acid receptors involved in inflammation and oxidative stress during heat stress [79,80]. Heat stress increases the expression levels of certain genes involved in promoted fatty acid metabolism, such as PPARγ, CCAAT-enhancerbinding protein α (CEBPα), and fatty acid synthase (FAS) [80]. Importantly, the PPAR signaling pathway activated by fatty acids regulates inflammatory responses and lipid metabolism [79,81,82]. The fatty acid metabolism-related genes and pathways identified in this study show the complexity of heat stress response and its possible implications for the welfare and production traits of livestock in summer. There were also 11 and 67 genes in the H30 and H60 groups significantly clustered in the thermogenesis pathway (ko04714), which help animals adapt to environmental changes. In adrenal glands, 47, 20 and 19 pathways were significantly detected by 1234, 1416, and 1226 genes in the CT vs H30, CT vs. H60, and CT vs. H120 comparisons (Table S8). Only two pathways, including It is widely known that changes in gene expression may be triggered by thermal stress; in turn, the expression of genes may help cells to defend themselves against thermal stress. Some genes were both activated by heat and cold stress, such as the HSP70 (HSP72), HSP90, HSF-1, and TRP channels [75]. A previous study showed that the trimerization of HSF-1 was induced by cold shock at 4 • C, and then increased the expression of HSP70 and HSP90 by binding to the heat response element [76]. This study also reported that HSF-1 not only mediated the expression of HSP during the period of heat stress but also after rewarming to 37 • C. The TRP channels are often activated by a variety of stressors (e.g., heat, cold, pH). TRPV1-4 and TRPM2 are activated by heat stress, and TRPA1 and TRPM8 are activated by cold stress [77]. In the present study, HSF-1 and Trpv2 were found to be involved in the process of the response to heat-and cold-induced thermogenesis at H60 and H120, suggesting that both HSF-1 and Trpv2 played important roles in thermoregulation, but further validation experiments are needed to confirm their functions. The Trpv2 gene can be considered as a novel target for a subsequent heat stress study.
In the liver, 31, 46, and 51 pathways were significantly enriched in the CT vs. H30, CT vs. H60, and CT vs. H120 comparisons, respectively (Table S8). Eleven pathways were commonly shared in the three comparisons, including fatty acid metabolism (ko01212) and the PPAR signaling pathway (ko03320) (Table S8). Fatty acids can regulate immune homeostasis by the activation of fatty acid receptors involved in inflammation and oxidative stress during heat stress [79,80]. Heat stress increases the expression levels of certain genes involved in promoted fatty acid metabolism, such as PPARγ, CCAAT-enhancerbinding protein α (CEBPα), and fatty acid synthase (FAS) [80]. Importantly, the PPAR signaling pathway activated by fatty acids regulates inflammatory responses and lipid metabolism [79,81,82]. The fatty acid metabolism-related genes and pathways identified in this study show the complexity of heat stress response and its possible implications for the welfare and production traits of livestock in summer. There were also 11 and 67 genes in the H30 and H60 groups significantly clustered in the thermogenesis pathway (ko04714), which help animals adapt to environmental changes. In adrenal glands, 47, 20, and 19 pathways were significantly detected by 1234, 1416, and 1226 genes in the CT vs. H30, CT vs. H60, and CT vs. H120 comparisons (Table S8). Only two pathways, including protein processing in the endoplasmic reticulum (ko04141) and ferroptosis (ko04216), were commonly detected by genes in all three comparisons. The pathway of thermogenesis (ko04714) was detected by 38 and 35 genes in the H30 and H60 groups, respectively. (Table S8).

Summary of Differential Isoform Corresponding Genes That Were also Differentially Expressed in Tissues
A total of 43, 786, and 1149 DEIDEGs were found in blood, liver, and adrenal gland tissues under heat stress, respectively (Table 5). Genes enriched in the heat stress process and thermogenesis are shown in Table 6. A total of 35 genes were enriched in the liver and adrenal glands, including 10 shared genes, such as Dnaja1, LOC680121, Chordc1, AABR07011951.1, Hsp90aa1, Hspa1b, Cdkn1a, Hmox1, Bag3, and Dnaja4, suggesting that these genes may regulate the heat stress response through crosstalk between the liver and adrenal glands. Combining with the expression profiles from the RNA-Seq and the validation results of the RT-qPCR experiments performed in our previous study [9], which demonstrated the high reliability of the RNA-Seq technology, we preferentially suggest that these 35 genes could be valuable candidates for underlying mechanisms to cope with heat stress.
Understanding mechanisms of heat stress response and discovering indicators to quantify the degree of heat stress are challenging for researchers [82,83]. Numerous studies have evidenced that gene expression analysis is a great strategy to identify positional and functional candidate genes associated with an organism's response to heat stress [45]. Furthermore, the development of RNA-Seq technology has greatly improved the accuracy and comprehensiveness of gene expression analyses and accelerated the progress of heat stress-related research [84,85]. In the current study, RNA-Seq analysis was performed on three tissues (blood, liver, and adrenal glands) and four treatment groups (CT, H30, H60, and H120) and was not only able to identify the known mRNA isoforms annotated to known genes but also detected some novel transcripts associated with non-annotated genes. Our study also indicated that the expression of heat stress-responsive genes has spatial and temporal differences [17]. Moreover, the novel mRNA isoforms identified in this study would contribute to further studies, providing candidate mechanisms related to heat stress response, and to enriching the rat transcriptome map. If these novel mRNA isoforms are common to different mammalians, these regions could be used to identify universal markers that could be used for breeding for heat-resistant livestock mammals. Note: DEIDEGs refer to the differential isoform corresponding genes that were also differentially expressed. B, L, and A mean blood, liver, and adrenal glands.  Note: FC means fold change, FDR means the false discovery rate, and L and A mean the liver and adrenal glands.

Conclusions
This research described the transcriptome profiles in the blood, liver, and adrenal glands of Sprague Dawley rats under various heat stress conditions, showing that the expression of heat stress-responsive genes in rats has spatiotemporal differences. In total, 35 genes were involved in the process of thermal stress regulation annotated by DEIs and were also differentially expressed in the liver and adrenal gland tissues, which could be preferentially considered as candidate markers for heat stress in rats. Furthermore, the identification of novel mRNA isoforms related to heat stress response allows for new insights into a better understanding of transcripts functionally playing important roles in heat stress, which may in the future be used to select for thermotolerance. The evident role of the inflammation-metabolism nexus appears to be central in the acute-heat-stressed rats through mRNA isoform identification, and hence this study indirectly emphasizes this intervention nexus for heat stress management. The magnitude of transcriptional changes uncovered in this study points towards the need for pathway-level biological investigation, along with further studies on the essential crosstalk between different biological pathways. Future GWAS studies focusing on the thermotolerance breeding of livestock can benefit from this study in the selection of genes, and a careful approach taking into account the spatiotemporal and post-transcriptional changes can be adopted.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology11121740/s1, Figure S1. Principal components analysis (PCA) (A) and Pearson correlation analysis (B) between the CT and heat stress groups based on the expression levels of genes in blood, liver, and adrenal gland tissues. The sample groups were manually preset in the R code for plotting PCA charts with such designed clusters. No unsupervised clustering method was used on the PCA results. Table S1. Overview of the mapping information in blood, liver, and adrenal gland tissues under various heat stress conditions. Table S2. The differentially expressed mRNA isoforms (p < 0.05, |FC| > 2) identified in blood in CT vs. H120. Table S3. The known differentially expressed mRNA isoforms (p < 0.001, FDR < 0.05, and |FC| > 2) identified in the liver under various heat stress conditions. Table S4. The known differentially expressed mRNA isoforms (p < 0.001, FDR < 0.05, and |FC| > 2) identified in adrenal glands under various heat stress conditions. Table S5. The novel differentially expressed transcripts (p < 0.001, FDR < 0.05 and |FC| > 2) length associated with annotated genes identified in the liver under various heat stress conditions. Table S6. The novel differentially expressed transcripts (p < 0.001, FDR < 0.05 and |FC| > 2) length associated with annotated genes identified in adrenal glands under various heat stress conditions. Table S7. Significantly enriched biological process (BP) terms from differentially expressed transcripts corresponding to known genes identified in blood, liver, and adrenal glands. Table S8. Significantly enriched pathways from differentially expressed transcripts corresponding to known genes identified in blood, liver, and adrenal glands. Table S9. Statistics of differentially expressed genes that were also annotated by differential mRNA isoforms (DEIDEGs) in blood, liver, and adrenal glands under heat stress conditions.

Data Availability Statement:
The RNA-Sequencing datasets of the CT and H120 groups (blood, liver, and adrenal glands) and datasets of the H30 and H60 groups (liver and adrenal glands) used in this work are available in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) with BioProject IDs PRJNA589869 and PRJNA690189, respectively.