Transcriptomic Analysis Provides Novel Insights into Heat Stress Responses in Sheep

Simple Summary The general increase in global temperatures has meant that heat stress has become an increasingly significant problem for sheep. This has both direct and indirect impact on their physiological functions, productivity, and health of sheep. Sheep generally live in high-temperature environments; however, the genes and pathways that play regulatory roles in the heat stress responses of sheep remain unclear. In this study, we applied RNA-Seq technology to analyze liver tissues of sheep from heat-stressed and control groups, and screened genes and pathways related to sheep heat stress. This work provides a theoretical foundation for the breeding and production of heat-resistant sheep. Abstract With the intensified and large-scale development of sheep husbandry and global warming, sheep heat stress has become an increasingly important issue. However, little is known about the molecular mechanisms related to sheep responses to heat stress. In this study, transcriptomic analysis of liver tissues of sheep in the presence and absence of heat stress was conducted, with the goal of identifying genes and pathways related to regulation when under such stress. After a comparison with the sheep reference genome, 440,226,436 clean reads were obtained from eight libraries. A p-value ≤ 0.05 and fold change ≥ 2 were taken as thresholds for categorizing differentially expressed genes, of which 1137 were identified. The accuracy and reliability of the RNA-Seq results were confirmed by qRT-PCR. The identified differentially expressed genes were significantly associated with 419 GO terms and 51 KEGG pathways, which suggested their participation in biological processes such as response to stress, immunoreaction, and fat metabolism. This study’s results provide a comprehensive overview of sheep heat stress-induced transcriptional expression patterns, laying a foundation for further analysis of the molecular mechanisms of sheep heat stress.


Introduction
Since its foundation in 1988, the Intergovernmental Panel on Climate Change has released five comprehensive climate reports and highlighted the indisputable fact that the world is warming. In the wild, air temperature has a major bearing on the ability of animals to thrive. Long-term high temperatures and excessively warm weather generate serious heat stress for animals. Heat stress directly or indirectly affects animal physiology, production, and health, among others [1][2][3], which is

RNA Extraction and Quality Inspection
Total RNA in the liver tissues of the eight sheep was extracted using TRlzol reagent (Invitrogen, Carlsbad, CA, USA), in accordance with the manufacturer's instructions. The integrity and concentration of RNA were determined by agarose gel electrophoresis and a Nanodrop 2000 (Thermo, Waltham, MA, USA). Finally, the RNA concentration was accurately quantified using Qubit 2.0 (Invitrogen, Carlsbad, CA, USA) and stored at −80 • C until subsequent analysis.

cDNA Library Construction and Sequencing
RNA samples that qualified through quality inspection were used for library construction in accordance with the instructions provided by Illumina (Illumina, San Diego, CA, USA). The specific steps were as follows: (1) Oligo-carrying magnetic beads were used for mRNA enrichment. (2) Fragmentation buffer was used to randomly break enriched mRNA into fragments of about 200 nucleotides in length.
(3) The fragmented mRNA was used as a template; one-chain cDNA was synthesized through inverse transcription using a random primer. In the synthesis of second-chain cDNA, dTTP in dNTPs was replaced by dUTP. (4) AMPure XP beads were used to purify double-chain cDNA, and End Repair Mix was used for end filling and the addition of an A-tail and sequencing linkers. (5) USER enzyme was used to digest double-stranded cDNA so that the library only contained single-stranded cDNA. (6) PCR enrichment of the cDNA fulfilling the above criteria was conducted. (7) Qubit 2.0 and Agilent 2010 (Agilent, Palo Alto, CA, USA) were used for quality inspection of the concentration and fragment size of the established library. The Illumina Hiseq 2500 platform was used to conduct double-end sequencing of the eight established libraries.

Raw Data Preprocessing and Alignment
Raw reads were obtained after Illumina sequencing ended. Skewer software was used to dynamically remove the 3' ends, linker sequences, and low-mass sequences of raw reads, after which clean reads were obtained. FastQC software was used to conduct quality control analysis of preprocessed data and determine indexes such as GC content, Q20, and Q30. STAR software was used to align obtained clean reads onto reference genomes of sheep (Oar_v4.0), where alignment parameters were set as -twopassMode Basic, -outSAMstrandField intronMotif, -alignSJstitchMismatchNmax 5, −1, 5, 5, and other parameters were set as the defaults. Those reads aligned to more than one (>1) genome were removed and the others were kept for subsequent analysis.

Functional Enrichment Analysis of Differentially Expressed Genes
The DAVID database was used in Gene Ontology (GO) functional enrichment analysis of the differentially expressed genes. All genes were placed on the background list; the differentially expressed genes were included in the candidate list, and hypergeometric distribution was used to calculate p values before multiple testing and Benjamini-Hochberg correction. GO items with a false discovery rate ≤ 0.05 were regarded as significantly enriched. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used in the functional annotation and classification of genes associated with differential peaks.

Verification of Differentially Expressed Genes
The RNA-Seq results were verified using qRT-PCR. A cDNA Synthesis Kit (TransStart Green, Beijing, China) was used for RNA reverse transcription. Oligo 7 software was used for cross-exon design of all primers and specific detection in NCBI. The qPCR SuperMix (TransStart Green, Beijing, China) and LightCycler 480 apparatus (Roche, Basel, Switzerland) were used for quantitative real-time PCR (qRT-PCR). Each sample was measured in triplicate to ensure the accuracy of the quantification. The relative expression of target genes was calculated by the 2 -∆∆Ct method using β-actin (ACTB) as an internal reference gene. Information on the primers used in this study is listed in Supplementary Table S1.

Statistical Analysis
All data reported are expressed as mean ± SE. Student's t-test was carried out using SPSS software for statistical analysis of the data. A p value of < 0.05 was considered to be statistically significant.

Influence of Heat Stress on Serum T 3 and T 4 Levels of Sheep
During heat stress, THI was ≥ 26.18, indicating that the sheep were exposed to extremely severe heat stress; in the control group, THI ≤ 11.82, indicating that the sheep experienced no heat stress (Supplementary Figure S1) [14]. T 3 and T 4 act in regulating the body's metabolic heat production. In this study, the levels of T3 and T 4 in the serum of heat-stressed and control sheep were determined and they were found to be extremely highly expressed in the control group sheep (Table 1).

Summary of Sheep Pen THI and Sequencing Data
The Illumina Hiseq 2500 platform was used to conduct paired-end sequencing of eight established libraries; the results of data quality control are shown in Table 2. A total of 513,487,656 raw reads were generated in the eight libraries, 440,226,436 clean reads were left after quality control, and the average proportion of clean reads was 99.06%. The GC contents of the eight samples were within 49.82%-51.25%, and were thus in accordance with base composition rules. Q20 ≥ 98.10% and Q30 ≥ 94.40%. Clean reads were aligned to the sheep reference genome using STAR software, and over 93.25% of the reads could be accurately aligned with a high matching rate. Approximately 7.05%-12.37% of these clean reads had Animals 2019, 9, 387 5 of 12 multiple aligned positions, 81.99%-86.73% of them had single aligned positions, and 3.66%-6.76% of them were not aligned to the sheep reference genome. Those reads that aligned on the sheep reference genome at a single position were used for further bioinformatic analysis. CG1 and CG2 were males in the control group and CM1 and CM2 were females in the control group, while HG1 and HG2 were heat-stressed males and HM1 and HM2 were heat-stressed females.

Analysis of Overall Gene Expression Levels
Through a comparison of the FPKM values of all genes, the overall gene expression levels of each sample could be evaluated. As shown in Figure 1, the overall gene expression levels of the eight samples were similar. To further evaluate the overall gene expression levels, the genes were divided into eight intervals (0.1-1, 1-5, 5-10, 10-20, 20-30, 30-40, 40-50, and ≥50) according to their FPKM values, and the number of transcripts in each interval and its percentage relative to the total number of transcripts were calculated. There were small percentage changes of various expression quantities in the eight samples, and the number of transcripts gradually decreased as the FPKM value increased (Supplementary Table S2). Overall, most transcripts were expressed at low levels in sheep liver in the heat-stressed group and control group (Supplementary Table S2). aligned positions, and 3.66%-6.76% of them were not aligned to the sheep reference genome. Those reads that aligned on the sheep reference genome at a single position were used for further bioinformatic analysis. CG1 and CG2 were males in the control group and CM1 and CM2 were females in the control group, while HG1 and HG2 were heat-stressed males and HM1 and HM2 were heat-stressed females.

Analysis of Overall Gene Expression Levels
Through a comparison of the FPKM values of all genes, the overall gene expression levels of each sample could be evaluated. As shown in Figure 1, the overall gene expression levels of the eight samples were similar. To further evaluate the overall gene expression levels, the genes were divided into eight intervals (0.1-1, 1-5, 5-10, 10-20, 20-30, 30-40, 40-50, and ≥50) according to their FPKM values, and the number of transcripts in each interval and its percentage relative to the total number Figure 1. Distribution of gene expression levels. CG1 and CG2 were males in the control group and CM1 and CM2 were females in the control group, while HG1 and HG2 were heat-stressed males and HM1 and HM2 were heat-stressed females.

Screening and Clustering Analysis of Differentially Expressed Genes
DESeq2 software was used to screen out sheep genes differentially expressed between the heat-stressed group and the control group; 1,137 genes fulfilling both criteria of p-value ≤ 0.05 and fold Animals 2019, 9, 387 6 of 12 change ≥ 2 were identified (Figure 2A and Supplementary Table S3). Compared with the levels in the control group, genes were downregulated ( Figure 2B). To obtain a deeper understanding of the gene expression patterns of sheep in the heat-stressed group and control group, clustering analysis of these differentially expressed genes was conducted. The results indicated that the differentially expressed genes of sheep in the heat-stressed group and control group were clustered into a single class ( Figure 2C). In addition, genes with upregulated differential expression and those with downregulated differential expression were clustered into one class. We also analyzed the effects of heat stress on gene expression in sheep of different genders. 4,203 genes were found to be differentially expressed between heat-stressed rams and rams from the control group

Functional Enrichment Analysis of Differentially Expressed Genes
Sheep genes that were differentially expressed between the heat-stressed group and the control group, were significantly enriched in 419 GO terms, including 357 biological processes, 27 cellular components, and 35 molecular functions (Supplementary Table S6). These GO terms are associated with stress response, immunoreaction, and fat metabolism, among others ( Figure 3A). KEGG pathway analysis found that these differentially expressed genes were significantly enriched in 51 pathways (Supplementary Table S7). Through a further analysis, many of these pathways were shown to be related to fat metabolism, such as regulation of lipolysis in adipocytes, MAPK 798 genes in sheep in the heat-stressed group presented upregulated expression, while 339 genes were downregulated ( Figure 2B). To obtain a deeper understanding of the gene expression patterns of sheep in the heat-stressed group and control group, clustering analysis of these differentially expressed genes was conducted. The results indicated that the differentially expressed genes of sheep in the heat-stressed group and control group were clustered into a single class ( Figure 2C). In addition, genes with upregulated differential expression and those with downregulated differential expression were clustered into one class. We also analyzed the effects of heat stress on gene expression in sheep of different genders. 4,203 genes were found to be differentially expressed between heat-stressed rams and rams from the control group (Supplementary Table S4), while 1,369 were differentially expressed between heat-stressed ewes and control group ewes (Supplementary Table S5).

Functional Enrichment Analysis of Differentially Expressed Genes
Sheep genes that were differentially expressed between the heat-stressed group and the control group, were significantly enriched in 419 GO terms, including 357 biological processes, 27 cellular components, and 35 molecular functions (Supplementary Table S6). These GO terms are associated with stress response, immunoreaction, and fat metabolism, among others ( Figure 3A). KEGG pathway analysis found that these differentially expressed genes were significantly enriched in 51 pathways (Supplementary Table S7). Through a further analysis, many of these pathways were shown to be related to fat metabolism, such as regulation of lipolysis in adipocytes, MAPK (mitogen-activated protein kinase), and PI3K-AKt (phosphatidylinositol 3-kinase and protein kinase B) ( Figure 3B). They were also particularly associated with signaling pathways regulating stress reactions, including TNF (tumor necrosis factor) and Rap1 ( Figure 3B). Besides, genes differentially expressed in heat-stressed rams and control group rams were particularly associated with the process of fat metabolism and stress response (Supplementary Tables S8 and S9). The results for genes differentially expressed between heat-stressed Animals 2019, 9, 387 7 of 12 ewes and control group ewes were also particularly associated with the process of fat metabolism and stress response (Supplementary Tables S10 and S11). (mitogen-activated protein kinase), and PI3K-AKt (phosphatidylinositol 3-kinase and protein kinase B) ( Figure 3B). They were also particularly associated with signaling pathways regulating stress reactions, including TNF (tumor necrosis factor) and Rap1 ( Figure 3B). Besides, genes differentially expressed in heat-stressed rams and control group rams were particularly associated with the process of fat metabolism and stress response (Supplementary Table S8 and Table S9). The results for genes differentially expressed between heat-stressed ewes and control group ewes were also

Validation of RNA-Seq data by qRT-PCR
To verify the sequencing results of RNA-Seq, five genes with high expression and five with low expression in the heat-stressed group in this study were selected for confirmatory analysis by qRT-PCR. As shown in Figure 4, although their expression levels differed to a certain degree, their expression patterns were identical, thus verifying that the RNA-Seq results were accurate and reliable.

Validation of RNA-Seq Data by qRT-PCR
To verify the sequencing results of RNA-Seq, five genes with high expression and five with low expression in the heat-stressed group in this study were selected for confirmatory analysis by qRT-PCR. As shown in Figure 4, although their expression levels differed to a certain degree, their expression patterns were identical, thus verifying that the RNA-Seq results were accurate and reliable.

Discussion
Heat stress impacts animal production and endangers animal welfare, causing an annual economic loss of over US$1.2 billion in animal husbandry globally [15]. In conventional breeding methods, the resistance of animals to heat can only be measured by determining some key production indexes, but the genetic changes of animals in association with environmental changes

Discussion
Heat stress impacts animal production and endangers animal welfare, causing an annual economic loss of over US$1.2 billion in animal husbandry globally [15]. In conventional breeding methods, the resistance of animals to heat can only be measured by determining some key production indexes, but the genetic changes of animals in association with environmental changes cannot be determined. In contrast, high-throughput sequencing can be used to shed light on animal heat resistance and genes associated with sensitivity to heat, and it has been proven to be a feasible approach in studies of animals such as cow, pig, and chicken [16][17][18]. Therefore, starting from transcriptomics, this study was aimed at identifying sheep heat resistance-related genes and pathways. Four sheep under heat stress and four without such stress were selected in this study, from which 1,137 differentially expressed genes were screened out using RNA-Seq. qRT-PCR results revealed expression patterns of 10 genes that were identical to those in the RNA-Seq results, proving that the sequencing results were accurate and reliable. These differentially expressed genes were mainly related to stress response, energy metabolism, and immunity. Further functional enrichment analysis of these differentially expressed genes was carried out, which showed that most genes were particularly associated with signaling pathways such as stress response, fat metabolism, and immunity.

T 3 and T 4
T 3 and T 4 are hormones secreted from the thyroid gland and its peripheral tissues. T 3 is a type of thyroid hormone that works on anabolism and catabolism. The reduction in secretion of T 3 is considered to reflect an animal's efforts to reduce its body heat production and accumulation as well as to maintain its body heat balance [1,19]. T 4 also plays an important role in the adaptation of animals to environmental changes. It stimulates intracellular oxygen consumption and heat production, thereby increasing basal metabolic rate, increasing glucose utilization, and altering lipid metabolism [19,20]. In this study, it was found that the concentrations of T 3 and T 4 in the blood of heat-stressed sheep were significantly lower than those in the sheep of the control group, which was consistent with the results in studies of other sheep breeds [7]. A decrease in the concentrations of T 3 and T 4 indicates that heat stress reduces the metabolic activity of the sheep, so as to inhibit the body from producing more heat. In addition, the decrease in thyroxine concentration may also be related to a decrease in blood glucose concentration, which is closely related to energy metabolism during heat stress.

Regulation of Body Temperature
5-Hydroxytryptamine receptor 4 (HTR4) presented significantly high expression in the heat-stressed sheep in this study, while HTR1B was significantly more highly expressed in sheep in the control group. It has been indicated that 5-HT plays a significant role in the process of regulating animal body temperature [21]. When animals are in a high-temperature environment, the 5-HT neuroendocrine system is activated and combines with 5-HTR to activate downstream cyclic adenosine monophosphate (cAMP) and cyclic guanosine monophosphate (cGMP) signaling pathways, which enhance neural activity sensitive to heat so as to inhibit heat production [22]. Fewer studies have been performed on HTR4 with regard to body temperature regulation, but it has been shown in animals to exert important regulatory effects in gastrointestinal sensitivity, memory, and food intake, among others [23]. Under heat stress, sheep exhibited reduced food intake and rumination activities, which further reduced the amount of heat generated by metabolism so as to maintain a stable body temperature [6,9]. During this period, HTR4 might play a significant role. HTR1B is also involved in body temperature regulation, and was shown to be activated to reduce the body temperature of rats [24]. Here, HTR1B presented high expression in sheep in the control group, suggesting its importance in regulating body temperature.

Regulation of Stress Reactions
Under heat stress, animals need to undertake a series of nonspecific responses to maintain homeostasis. The functional enrichment analysis of differentially expressed genes in sheep in the heat-stressed group and the control group showed the significant enrichment of genes in stress-related pathways, such as Rap1, MAPK, and PI3K-Akt. Rap1 is regulated by a wide range of external stimuli as a molecular switch when stress occurs [25]. MAPK is one of the main signaling pathways in mammals and can also be activated by external stimuli, but the activity of the MAPK signaling pathway is stimulated by Rap1 in multiple cells [26,27]. In a similar way, the PI3K-Akt pathway is activated by external stimuli of many cellular types and can involve a cascade reaction with the MAPK signaling pathway [28,29]. In summary, these stress-related signaling pathways can form a complicated cascade reaction to respond to heat stress in sheep, and reduce the harm brought about by such stress.

Regulation of Energy Metabolism
Under heat stress, sheep reduce their food intake in order to reduce metabolic heat production. This study showed that differentially expressed genes in the heat-stressed group were particularly associated with fat metabolism-related pathways. Besides regulating stress reactions, MAPK and PI3K-Akt signaling pathways exert important effects on the fat metabolic process [30,31]. In addition, the cAMP signaling pathway also plays a significant role in the hepatic fat metabolic process and can activate protein kinase A (PKA) to promote β-oxidation of hepatic fat [32]. Genes associated with these signaling pathways were also particularly common among the differentially expressed ones, such as NPR1 (natriuretic peptide receptor 1), ANGPT2 (angiopoietin 2), and SLC13A5 (sodium-dependent citrate transporter member 5); this study indicated that all of these genes participate in the regulation of fat metabolism [33][34][35]. Sheep needed to mobilize a large amount of energy to respond to the heat stress, which enhanced catabolism and reduced anabolism correspondingly. Although genes associated with some glycometabolism-related signaling pathways were also commonly identified in this study, compared with fat, saccharide energy content was low with easy hydration. Therefore, fat mobilization in sheep under heat stress could help them adapt to this condition.

Regulation of Immunoreactions
Moderate heat stress has the effect of promoting immune system activity in animals so that they can better adapt to the environment. However, severe heat stress has a negative influence on the immune systems of animals and results in deterioration of their disease resistance. IL1R (interleukin type 1 receptor) has two subtypes, IL1R1 and IL1R2, both of which exert important effects on immune and inflammatory reactions [36]. HSPA2 belongs to the HSP70 (heat stress protein 70) family. It has been shown that HSP70 can increase the secretion of IL-1β from cells and facilitate the transport of antigens to T cells [37][38][39]. In this study, it was found that IL1R1, IL1R2, and HSPA2 presented high expression in heat-stressed sheep, suggesting that these genes play significant roles in the immunoregulatory process of sheep.

Conclusions
Upon RNA-Seq analysis of eight libraries of heat-stressed and control groups, 1,137 differentially expressed genes were screened out. These genes are mainly involved in biological processes such as stimulus-response reaction, immunoreaction, and fat metabolism. These candidate genes and pathways should deepen our understanding of the molecular mechanisms involved in heat stress in sheep, and provide a theoretical basis for improved sheep production and breeding of heat-tolerant varieties.
Availability of Supporting Data: All the RNA-Sequencing data used in this study has been deposited in the Sequence Read Archive (SRA) public databases under BioProject (PRJNA531367).

Supplementary Materials:
The following are available online at http://www.mdpi.com/2076-2615/9/6/387/s1, Figure S1: THI of sheep 1 week before slaughter, Table S1: Primers used in this study for qRT-PCR, Table S2: Number of genes in different FPKM intervals, Table S3: Differentially expressed genes identified in heat-stressed and control sheep groups by RNA-Seq, Table S4: Differentially expressed genes identified in heat-stressed rams and control group rams, Table S5: Differentially expressed genes identified in heat-stressed ewes and control group ewes, Table S6: GO analysis of differentially expressed genes in heat-stressed and control sheep groups, Table S7: KEGG analysis of differentially expressed genes in heat-stressed and control sheep groups. Table S8: GO analysis of differentially expressed genes in heat-stressed rams and control group rams, Table S9: KEGG analysis of differentially expressed genes in heat-stressed rams and control group rams, Table S10: GO analysis of differentially expressed genes in heat-stressed ewes and control group ewes, Table S11: KEGG analysis of differentially expressed genes in heat-stressed ewes and control group ewes.