Next Article in Journal
Transcriptomic and Metabolomic Analyses Reveal the Roles of Flavonoids and Auxin on Peanut Nodulation
Next Article in Special Issue
Genome-Wide DNA Methylation and Transcriptome Integration Associates DNA Methylation Changes with Bovine Subclinical Mastitis Caused by Staphylococcus chromogenes
Previous Article in Journal
Investigating the Physiological and Molecular Responses of Solanum lycopersicum hp Mutants to Light of Different Quality for Biotechnological Applications
Previous Article in Special Issue
A Satellite-Free Centromere in Equus przewalskii Chromosome 10
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 24
Issue 12
10.3390/ijms241210150
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analysis of CircRNA Expression in Peripheral Blood of Holstein Cows in Response to Heat Stress

by
Congcong Zhang
1,
Shuhui Wang
1,†,
Lirong Hu
2,
Hao Fang
1,
Gong Chen
1,
Xiaojuan Ma
1,
Ying Yu
2,
Yachun Wang
2,* and
Qing Xu
1,*
1
College of Life Sciences and Bioengineering, Beijing Jiaotong University, Beijing 100044, China
2
National Engineering Laboratory for Animal Breeding, Key Laboratory of Animal Genetics, Breeding and Reproduction, MARA, College of Animal Sciences and Technology, China Agricultural University, Beijing 100193, China
*
Authors to whom correspondence should be addressed.
†
Current address: Laboratoire d’Analyse et de Modélisation de Systèmes pour l’Aide à la Décision, Paris Dauphine University—PSL, Pl. du Maréchal de Lattre de Tassigny, 75016 Paris, France.
Int. J. Mol. Sci. 2023, 24(12), 10150; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms241210150
Submission received: 6 May 2023 / Revised: 6 June 2023 / Accepted: 8 June 2023 / Published: 15 June 2023
(This article belongs to the Special Issue Animal Genomes and Epigenomes)
Browse Figures
Review Reports Versions Notes

Abstract

:
The present study aimed to identify key circRNAs and pathways associated with heat stress in blood samples of Holstein cows, which will provide new insights into the molecular mechanisms driving heat stress in cows. Hence, we evaluated changes in milk yield, rectal temperature, and respiratory rate of experimental cows between heat stress (summer) and non-heat stress (spring) conditions with two comparisons, including Sum1 vs. Spr1 (same lactation stage, different individuals, 15 cows per group) and Sum1 vs. Spr2 (same individual, different lactation stages, 15 cows per group). Compared to both Spr1 and Spr2, cows in the Sum1 group had a significantly lower milk yield, while rectal temperature and respiratory rate were significantly higher (p < 0.05), indicating that cows in the Sum1 group were experiencing heat stress. In each group, five animals were chosen randomly to undergo RNA-seq. The results reveal that 140 and 205 differentially expressed (DE) circRNAs were screened in the first and second comparisons, respectively. According to the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, these DE circRNAs were mainly enriched in five signaling pathways, including choline metabolism, the PI3K/AKT signaling pathway, the HIF-1 signaling pathway, the longevity-regulating pathway, and autophagy. Then, we obtained the top 10 hub source genes of circRNAs according to the protein–protein interaction networks. Among them, ciRNA1282 (HIF1A), circRNA4205 (NR3C1), and circRNA12923 (ROCK1) were enriched in multiple pathways and identified as binding multiple miRNAs. These key circRNAs may play an important role in the heat stress responses of dairy cows. These results provide valuable information on the involvement of key circRNAs and their expression pattern in the heat stress response of cows.

1. Introduction

With the global temperature rise, the dairy industry is confronted with a significant challenge brought about by severely damaging impacts of heat stress. Heat stress increases mortality rates and threatens the health of livestock, especially in highly productive and intolerant species. The increase in ambient temperature–humidity index (THI) in summer will affect the balance of organisms and lead to heat stress in animals. Dairy cows are more sensitive to high THI due to their high metabolism and low heat dissipation. Heat stress in dairy cows has become increasingly severe, causing huge economic losses to the dairy industry. It has been reported that heat stress has caused a loss to the dairy industry of $1.5 billion each year [1]. Therefore, it is necessary to develop strategies to reduce the damage caused by heat stress in animal husbandry. Investigating the effect of heat stress on cow performance and its potential mechanisms is very important to fundamentally reduce the negative effect of heat stress. With the rapid development of sequencing technology, genome-wide association studies (GWAS) [2] and transcriptome [3], metabolome [4,5], and microbiome [6,7] analyses were used to screen biochemical indicators and identify molecular mechanisms associated with heat stress in dairy cows. However, there are few studies about the regulation of post-transcriptional noncoding RNAs (ncRNAs), such as long noncoding RNAs (lncRNA), microRNAs (miRNAs), and circular RNAs (circRNAs), in dairy cows after exposure to heat stress conditions, especially studies focusing on circRNAs.
Recently, circular RNA (circRNA) has attracted significant attention in research due to its potentially rich functions in cells. CircRNAs are a large class of ncRNAs ubiquitously expressed in eukaryotic cells. They have a covalently bonded closed-loop structure, resulting in their stronger stability than linear RNAs [8]. Subsequent research has tried to uncover the functional role of circRNAs in various species and tissues. A handful of circRNAs have been demonstrated to have functions in post-transcriptional regulation, acting as miRNA sponges [9]. An increasing number of studies have been conducted to improve livestock milk production [10], meat quality [11] and skeletal muscle development [12] by describing the mechanisms of action and related pathways of circRNAs. Hao et al. [10] found 18 classes of miRNAs play essential roles in the regulation of lactation and mammary growth through targeting circRNAs in sheep. CircRNA is also associated with the development of livestock embryos, and circTUT7 regulates the expression of HMG20B in a ceRNA mechanism, suggesting their critical functions in embryonic skeletal muscle development in porcine embryonic muscle [13]. Regarding the role of circRNAs in meat quality in cows, circFUT10 was observed to bind miR-133a to regulate myoblast differentiation in bovine skeletal muscle tissue [14]. However, little is known about the role of circular RNAs in the responses of dairy cows undergoing heat stress.
Limited studies on circRNA and heat stress in dairy cows have been reported. Wang et al. [15,16] identified several circRNA-associated ceRNA networks involved in milk fat metabolism under heat stress and found circEZH2 regulated milk fat metabolism through miR-378b sponge activity. Furthermore, copious changed circRNAs were discovered in cow mammary gland tissues [15,17] and hypothalamic and pituitary tissues [18] under heat stress. Among all these papers, different tissue samples of cattle under heat stress were used, but not peripheral blood. It is worth mentioning that other cattle tissue samples are more difficult to obtain than peripheral blood, since it is readily available from cows. Moreover, sampling of peripheral blood is less invasive, more convenient, and safer. Especially as one common cellular model, the role and mechanism of cattle peripheral blood mononuclear cells (PBMCs) have been extensively studied [19,20]. In addition, when designing experiments, we had to consider the effect of various factors on the incidence of heat stress, such as individual genetic background, lactation, and breed. Most studies on heat-stress-related circRNAs have selected cows with the same lactation stage as samples, which may ignore the individual differences [15,17]. Liu et al. [21] described 86 expression patterns of mRNA and miRNA in the blood of the same individual cows under different conditions from heat stress to thermal equilibrium. They suggested individual influence on cows’ response to heat stress. Thereby, we investigated the effect of heat stress on circRNA by two experimental comparisons, the same individuals under different lactation stages in comparison1 (Spr1 and Sum1 groups) and different individuals under the same lactation stage in comparison2 (Spr2 and Sum1 groups).
Here, we examined the changes in milk yield, rectal temperature, and respiratory rate of dairy cows under heat stress. Then, we focused on the effect of heat stress on the key biological processes and candidate circRNAs in blood samples of cows through two experimental comparisons. The obtained results comprehensively elucidate the underlying circRNAs of heat stress response and provide a new insight into breeding of thermotolerant dairy cows.

2. Results

2.1. The Performances of Holstein Cattle under Heat Stress

The experiments were conducted in April and July 2017. Between the Spr (April) and Sum (July) groups, there was obvious difference in the temperature–humidity index (THI) calculated from ambient temperature and humidity (p < 0.05). The average THI in spring was about 55, while it increased to about 81 in summer (Table 1). From THI, it can be seen that no heat stress occurred in 2 Spr groups in April, when an average THI only reached to 55. In contrast, the Sum1 group suffered moderate heat stress during July, since the THI reached a high value: 81. Additionally, from Table 1, it is evident that the seven-day average milk yield (7AMY) and milk yield (MY) were generally higher in the Spr1 group than in the Sum1 group (p < 0.05). Furthermore, both physiological indicators, rectal temperature (RT) and respiration rate (RR) of cow, were significantly increased in the Sum1 group compared with the Spr2 group.

2.2. Identification and Characterization of circRNAs in Holstein Cattle

In order to profile the expression of circRNAs in Holstein cattle, 15 samples from Spr1, Spr2, and Sum1 groups were used to perform RNA-seq and extract the characteristics of circRNAs. After removing the de-duplicated reads, we obtained an average of 74,665,172 reads in each sample, of which 91.45–93.96% of the reads were successfully mapped to the bovine reference genome, and 68.68–78.32% were uniquely mapped to the genome (Table S1). These results show our sequencing data had good quality, which supported its usage for the subsequent analysis. After removing the linear RNA and ribosomal RNA, 24,250, 26,201, and 24,848 circRNAs were identified in blood samples of Holstein cows from Spr1, Spr2, and Sum1 groups, respectively (Figure 1a). The majority of circRNAs were composed of exons (average = 47.93%) and introns (average = 44.32%), while a small proportion of circRNAs (average = 7.75%) originated from intergenic sequences (Figure 1b). Figure 1c mainly shows the distribution of circRNAs on chromosomes. It can be seen that these detected circRNAs in the current study originated from all chromosomes of cows, while the largest number of circRNAs came from chromosome 1. In addition, the number of exons contained in a circRNA is examined. From Figure 1d, it can be observed that although the majority of them consisted of multiple exons, more than 20% of circRNAs were composed of a single exon.

2.3. Analysis of Differential Expression Patterns of circRNAs

Using SRPBM as the metric, the overall expression of circRNAs in each cow is presented in Figure 2a. With the log2 (fold change) ≥ 1 and p < 0.05 as a threshold to screen the DE circRNAs, 140 circRNAs were differentially expressed in Sum1 and Spr1 groups (comparison1, 71 of them were up-regulated and 69 were down-regulated, Figure 2b), while 205 circRNAs were differentially expressed in Sum1 and Spr2 groups (comparison2, 83 of them were up-regulated and 122 were down-regulated, Figure 2c). A statistical summary of RNA sequencing data for each comparison is shown in Tables S2 and S3. Furthermore, 34 circRNAs were differentially expressed in both comparisons (Figure 2d, Table S4). More circRNAs showed differential expression between Sum1 and Spr1 groups, compared to between Sum1 and Spr2 groups. A hierarchical clustering heatmap of the DE circRNAs showed that the expression patterns of the circRNAs were clearly differentiated and aggregated between Sum1 and Spr1 groups, as well as between Sum1 and Spr2 groups (Figure 2e,f). These results suggest that heat stress greatly modulated the expression of circRNAs in cows’ blood. Moreover, as for the two groups with the same cows in different lactation stage in comparison, their DE circRNAs were due to different heat stress and lactation stage, while regarding the two groups with different cows in the same lactation stage in comparison2, the DE circRNAs were mainly due to different heat stress and individuals. Accordingly, 34 overlapped DE circRNAs in 2 comparisons should be more highly correlated with heat stress, since the circRNAs’ abundance was affected by either lactation stage or individual genetic background.

2.4. GO and KEGG Analyses of Source Genes of DE circRNAs

In order to investigate the potential biological functions of source genes of DE circRNA, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were carried out with source genes of 140 and 205 DE circRNAs in the first and second comparisons, respectively. GO categories were assigned to the source genes of circRNAs, and the genes were classified into three GO categories: cellular component, biological process, and molecular function (Figure 3a and Figure 4a). Summary statistics for GO and KEGG entries for each comparison are shown in Tables S5–S8. The current study presents only the top 50 GO terms with significance (p < 0.05) and the top 20 KEGG pathways. Based on the enrichment results of parental genes from DE circRNAs, RNA polymerase ΙΙ regulatory region sequence-specific DNA binding, protein import into nucleus, autophagy, and ion transport GO terms were detected in Sum1 and Spr1 groups (Figure 3a). Meanwhile, the pathways significantly changed by the DE circRNAs in Sum1 vs. Spr1 were the sphingolipid signaling pathway, cAMP signaling pathway, and PI3K-Akt signaling pathway (Figure 3b). Similar to the previous analysis, the GO and KEGG enrichment for Sum1 vs. Spr2 groups is shown in Figure 4. From Figure 4a, we can find that double-strand break repair via homologous recombination, innate immune response, and protein autophosphorylation items are most enriched. For the KEGG pathway, they are most significantly enriched in the autophagy–animal pathway, transcriptional misregulation in cancer, and the thyroid hormone signaling pathway (Figure 4b).
Meanwhile, the significantly changed pathways related to heat stress and their enriched DE circRNAs in two comparisons were visualized by using the OmicStudio tools. As shown in Figure 5, there were 15 significantly changed pathways in Sum1 vs. Spr1 groups and 9 in Sum1 vs. Spr2 groups (p < 0.05). These pathways were enriched for 39 parental genes. Among 15 pathways in Sum1 vs. Spr1 groups, the sphingolipid signaling pathway showed the most significance (p = 0.001) and contained five enriched DE circRNAs: circRNA12069 (PTEN), circRNA12923 (ROCK1), circRNA2560 (NFKB1), circRNA3691 (PIK3CB), and circRNA5055 (PLD1). Furthermore, circRNA2560 (NFKB1) and circRNA3691 (PIK3CB) were enriched in the longevity-regulating pathway, cAMP signaling pathway, and HIF-1 signaling pathway, which played important roles in energy metabolic regulation in cattle under heat stress. Meanwhile, among 9 pathways from Sum1 vs. Spr2 groups, the autophagy–animal pathway was the most significantly enriched functional pathway (p = 0.009), which included 6 DE circRNAs: circRNA1153 (AMBRA1), circRNA11953 (MTMR3), circRNA13493 (RHEB), circRNA5260 (SH3GLB1), circRNA6642 (AKT3), and ciRNA1282 (HIF1A).

2.5. Analyses of the Overlapped DE circRNAs in Two Comparisons

There are 34 circRNAs overlapped in 2 comparisons. The GO and KEGG annotations for these two differentially expressed circRNAs are displayed in Table 2 and Table 3. It was found that most of these differential circRNAs were enriched in pathways related to energy metabolism and RNA processing. In the biological process category, two important GO subcategories, guanyl-nucleotide exchange factor activity and cellular response to transforming growth factor beta stimulus, were annotated. In the molecular function category, the small GTPase-mediated signal transduction and repressing transcription factor binding subcategories were also annotated. The biological interpretations of the circRNA source genes were further analyzed using the KEGG pathway database, and three KEGG pathways were enriched: choline metabolism in cancer, endocytosis, and the longevity-regulating pathway. CircRNA27183 (DOCK7), circRNA13493 (RHEB), CiRNA3301 (WWOX), and circRNA4205 (NR3C1) were enriched in multiple GO and KEGG entries and may exhibit functions in response to heat stress of cattle.

2.6. Construction of the Protein–Protein Interaction Network and Determination of Hub Source Genes of circRNA

Using the protein interaction network, the protein interaction relationships between different genes can be determined. With the cytoNCA plug-in of Cytoscape 3.9.1 online software, we determined which source gene of circRNA was in the hub position in the regulatory network. PPI networks were constructed on the basis of STRING database. A total of 39 KEGG enriched source genes of circRNA in the network (Figure 5) and 34 overlapped source genes of circRNA (Table 3) were selected to map the PPI network using String database (http://string-db.org (accessed on 20 March 2023)) (Figure 6a). In total, 57 nodes and 47 edges were displayed in the PPI network. Finally, 26 hub source genes of circRNA were found by cytoNCA (Cytoscape plug-in) (Figure 6b). The top 10 source genes of circRNA are circRNA12069 (PTEN), ciRNA1282 (HIF1A), circRNA14021 (PTPRC), circRNA13493 (RHEB), circRNA22544 (EZH2), ciRNA292 (RUNX2), circRNA2467 (NR3C1), circRNA4205 (NR3C1), circRNA3691 (PIK3CB), circRNA12923 (ROCK1), and circRNA5055 (PLD1). Detailed information regarding these 10 hub genes is listed in Table 4.

2.7. Interaction of circRNA and miRNA

A myriad of studies have shown that circRNAs can regulate miRNA as competitive endogenous RNAs. To further understand the function of circRNAs, we used both TargetScan and miRanda to predict the interactions between circRNAs and miRNAs. Then, the interaction network data files were generated and imported into OmicStudio tools. The attributes of the target circRNAs were visualized in the network. Based on all 345 differentially expressed circRNAs, 1917 miRNAs were estimated to have a tendency to interact with circRNAs (miRanda_score ≥ 145; miranda_Energy ≤ −20 Kcal/mol, Table S9). According to the above differential expression and GO/KEGG/PPI enrichment analyses, we identified 10 genes of different sources related to heat stress in cattle, corresponding to 11 circRNAs (Table 4). Figure 7 shows the miRNA–circRNA network diagram for all 11 circRNAs, in which circRNA292 interacted with 38 miRNAs, circRNA1282 interacted with 19 miRNAs, and circRNA2467 and 4205 interacted with 8 miRNAs. These results suggest that the identified circRNAs may be involved in heat-stress-related signaling pathways through sponging multiple miRNAs.

2.8. Verification of circRNA by RT-qPCR

In order to confirm the RNA-seq data, six circRNAs were randomly chosen for validation by PCR amplification with divergent primer or with one primer spanning the spliced junction. The gel electrophoresis results reveal that each circRNA had a single band at the expected length location (Figure 8a). With the intention of verifying that the observed bands were truly the corresponding circRNAs, the cDNA sequencing results were obtained. Because of the specificity of primers and the circular character of circRNAs, the back-splicing junction can be observed in circRNA products (Figure 8b). The sequencing results are in agreement with the expected RNA sequence, and the splice sites were spotted. Finally, we successfully confirmed the expressions of these six circRNAs by RT-qPCR. Results were compared with the high-throughput RNA-seq results, which showed that the expression of the six circRNAs was consistent with the trends obtained from RNA-seq data (Figure 8c).

3. Discussion

High ambient temperature has a profound effect on animal welfare. With escalating global temperature, compounded by the increased intensification of production, heat stress significantly impacts the dairy industry. Cows are homoeothermic animals, and their thermal comfort ranges from 5 °C to 25 °C [22]. Cows can maintain physiological stability through their own regulatory mechanisms, but when the external environmental temperature is higher than this range, cows will suffer from terrible heat stress [23,24]. Therefore, determining whether a cow is under heat stress and its response level has become a research focus. Various indicators have been established to assess the status of cows experiencing heat stress, such as the THI. Usually, THI = 72 is used as a threshold for judging the heat stress status of dairy cows. When THI is higher than 72, cows will encounter heat stress [25]. In this study, experimental cows were undergoing heat stress (Summer: THI = 81.08 ± 4.57), triggering a series of thermoregulation responses. Accompanied by reduced feed intake, increased water intake, and respiratory rate, cows exhibit a range of physiological activities [26]. The RT and RR, two important indicators of heat stress [27], were found to significantly increase (p < 0.05) in Sum1 cows, suggesting the involvement of physiological thermoregulation. We also found that heat stress had a significant effect on milk yield. It has been shown that cow milk production in the heat-stressed group decreased by about 10–20% compared to the group under appropriate temperature conditions [28,29], which is consistent with our results. These results reveal that Holstein cows in Beijing were undergoing severe temperature stress in July. Furthermore, heat stress had damaging effects on milk production traits and physiological condition of Holstein cows in the Sum1 group.
CircRNAs are a new class of endogenous non-coding RNAs that are covalently bound in a closed-loop structure and are more stable than linear RNAs [30]. Despite being discovered nearly 40 years ago, the significance of circular RNA has only recently been recognized, and circRNAs are now known to play a key role in gene regulation [31]. Due to their ability to bind to miRNA molecules, circRNAs can reduce the repressive effect of miRNAs on their target genes in the circRNA–miRNA–mRNA regulatory network [32]. Moreover, many reports have confirmed the role of circRNAs both in development and disease by analyzing their expression profiles [33,34]. However, little is known about the role of circular RNAs in cattle undergoing heat stress in peripheral blood. The present work investigated (i) changes in circRNAs that occurred in cattle under heat stress and (ii) key pathways and hub circRNAs involved in the cattle in response to heat stress.
To explore the regulatory effect of circRNAs on dairy cows under heat stress, we profiled the total expression of circRNAs in Holstein cattle with RNA-seq and extracted the characteristics of circRNAs in response to heat stress. According to GO and KEGG pathway analysis results, critical pathways related to the metabolic processes and autophagy in dairy cows under heat stress were identified. Although reduced dry matter intake due to heat stress has long been recognized to be the main cause of reduced milk production, a growing body of research is revealing that significant changes in the metabolism processes of cows in hot environmental conditions may be the key to milk production and milk composition [17,35]. In the present study, we found the key signaling pathways enriched in KEGG metabolic pathways in peripheral blood between the Spr and Sum groups, such as choline metabolism, the PI3K/AKT signaling pathway, the sphingolipid signaling pathway [36], and the cAMP signaling pathway [37], all of which have been previously reported to involved in heat stress response in other studies. Yang et al. [38] reported that treatment of mammary epithelial cells with choline effectively inhibited heat-stress-induced oxidative stress and apoptosis. In our results, overlapped DE in cows’ circRNAs under heat stress was most enriched in choline metabolism, which is related to alleviating heat-stress-induced damages by reducing oxidative stress [39]. Furthermore, PI3K/AKT [40], sphingolipid [36] and cAMP [37] have been proved to be important signaling molecules participating in adaptation to heat stress by affecting multi-level regulatory networks. From our results, the RT and RR significantly increased but the milk yield decreased with the THI in cows (Table 2), which revealed that Holstein cows in Beijing were undergoing severe temperature stress in July. Thus, cows tried their best to mitigate the damage caused by heat stress through a variety of metabolic processes.
Furthermore, heat stress can alter oxidative homeostasis and then cause cellular autophagy and apoptosis. This is consistent with the enriched KEGG pathways presented here: the HIF-1 signaling pathway, the longevity-regulating pathway, autophagy, and the mitophagy pathway. Heat stress was suggested to be responsible for inducing oxidative stress during summer in livestock animals. Later on, oxidative stress will result in the production of reactive oxygen species (ROS), alteration of ATP product, DNA damage, and inhibition of protein synthesis [41]. As the consumption of oxygen is elevated under thermal stress, limitations of oxygen may lead to a transient hypoxia. Additionally, the oxidative stress may activate the HIF-1 pathway [42]. Moreover, the enriched autophagy and mitophagy pathway are also highly related to oxidative stress triggered by heat stress. Baechler et al. [43] found that autophagy and mitophagy are two essential mechanisms regulating mitochondrial oxidative stress and mitochondria-associated cell death events.
Key DE circRNAs were identified as related to the metabolic processes and autophagy in cattle under heat stress. Among them, the expression of ciRNA1282 (HIF1A) was significantly changed in Sum1 compared with Spr2, and it binds to multiple miRNAs. Chen et al. [44] found that Hsp90 relieves heat-stress-induced damage by involvement of autophagic HIF-1α signaling, which suggested that HIF1A might alleviate the damage caused by heat stress in animals. Moreover, ciRNA1282 (HIF1A) might bind to the bta-let-7 family, in which Bta-let-7c has been shown to promote adipocyte proliferation and inhibit adipocyte differentiation by acting as a sponge for circFUT10 [45]. Furthermore, circRNA13493 (RHEB) and circRNA4205 (NR3C1) were differentially and consistently expressed in two groups. RHEB is a small H-Ras-like GTPase that is a direct and essential activator of the lysosomal surface rapamycin complex 1 (mTORC1) machinery, regulates cell growth, and activates autophagy [46,47]. Heat stress may activate the modulatory effects of RHEB on the mTOR pathway, thereby registering proliferation and autophagy in dairy mammary cells. As for NR3C1 (glucocorticoid receptor), it was found to reduce the promoter methylation under heat stress in dairy cattle [48]. Another interesting finding is that the expression of circRNA22544 (EZH2) significantly differed between Sum1 and Spr1 [16]. Wang et al. [16] found that heat-stress-related circEZH2 affected the proliferation, apoptosis, and lipid metabolism of mammary gland epithelial cells, which is consistent with our results.
Usually, in the investigation of heat stress of cattle, with the increase in experimental samples’ quantity, the analysis is statistically more robust and reliable. However, sampling and sequencing require a lot of manpower and expenses. In practice, generally more than three samples per group are selected for analysis [49,50]. In this manner, we randomly selected 5 animals per group of 15 animals for RNA-seq. Moreover, the comparison design can be constructed either by using cows with the same lactation stage in different temperatures, such as Spr2 and Sum1 in this study, or by using the same samples assessed in two seasons, such as Spr1 and Sum1 in this case (apparently, the lactation stage of samples is different). Between Sum1 and Spr1 groups, pathways closely related to energy metabolic and oxidative stress were found, such as the sphingolipid signaling pathway, cAMP signaling pathway, PI3K-Akt signaling pathway and HIF-1 signaling pathway. Liu et al. [21] selected five Holstein cows (female, healthy) and extracted blood samples at the same moment by coccygeal venipuncture in the summer (August) and winter (December) for RNA-seq. Their functional analyses showed that the MAPK signaling pathway, cellular senescence cGMP-PKG signaling pathway, and thermogenesis are related to oxidative damage and cell apoptosis. This result suggested that heat stress might have different impacts on cows in different stages of lactation. This might be due to the fact that heat stress can lead to an aggravation of lipid peroxidation and a decrease in total antioxidant capacity in dairy cows at different lactation stages. Moreover, heat stress negatively impacts feed intake, milk yield, and negative energy balance in early lactation [51]. However, in Sum1 and Spr2 groups, the largest number of DE circRNAs were enriched in heat-stress-related signaling pathways such as autophagy. Luo et al. [2] selected 8 primiparous (4 cows in April, 4 cows in July) cows with DIM ranging from 135 to 144 d (mid-lactation and pregnant) for RNA-seq. In their study, some genes associated with autophagy and apoptosis were identified. It was pointed out by Zhang et al. [52] that exposure of germ cells to hyperthermia resulted in several specific features of the autophagic process. Therefore, heat stress might mainly affect the expression of circRNA through autophagy in dairy cows at the same lactation stage and thus induce oxidative injury in paired samples.
From the above analysis, it can be seen that different comparison designs will lead to different results. It is important for researchers to choose an appropriate comparison design for heat stress study. When dairy cows are in lactation stage, they are more likely to be affected by heat stress. Consequently, their energy consumption is higher and their metabolic activity is more active [53]. When we are more concerned about the effects of heat stress on cows’ lactation performance, the impact on cow lactation cannot be ignored. In this case, comparison1 is probably a better choice. However, supposing that we are more concerned about the effects of heat stress on the cow as an organism, we should weigh the impact of different lactation stages on cows. Moreover, the breed of dairy cattle is also relevant to the choice of experimental samples. Holstein cows display fewer individual differences than local breeds, and thus, they are more appropriate for the design of comparison2. Our results provide the experimental design suggestions and basis to study dairy cows’ responses to heat stress.

4. Materials and Methods

4.1. Animals and Sampling

The experiment was performed in agreement with the Committee on Ethics of Animal Experimentation at the Beijing Jiaotong University, Beijing, China (Code ID: SS-QX-2014-06; 26 June 2014).
A total of 30 healthy primiparous Holstein lactating cows with similar age and body weight from a commercial farm (Beijing, China) were selected as experimental animals. All cows were raised in a loose pen with an exercise yard outside the house, free access to water, and a total mixed ration (TMR) throughout the day. What is more, the diet’s nutrient composition remained consistent throughout the test period, and its nutritional level was in accordance with the standard for dairy cattle feeding (NY/T 34-2004). Table S10 shows the composition and nutritional level of the feed. The temperature and humidity were automatically recorded every half hour. The temperature–humidity index (THI) was calculated according to the equation THI = 0.8AT + RH × (AT − 14.1) + 46.4 (AT and RH are temperature and relative humidity) [54]. Out of 30 cows, the Spr1 group (15 cows) and the Spr2 group (15 cows) were sampled in April 2017 (the thermoneutral period, THI = 55.43), and the Sum1 group (same 15 cows as in Spr1) was sampled in July 2017 (heat stress period, THI > 78). Among them, Spr1 (DIM = 61.6 ± 10.62 d) and Sum1 (DIM = 137.60 ± 10.62 d) were the same individuals, while Spr2 (DIM = 137.53 ± 7.63 d) and Sum1 (DIM = 137.60 ± 10.62 d) were in the same lactation stage.

4.2. Measurement of Phenotypic Data

During the experimental periods, daily milk yield (MY) for all 45 cows were collected and the average milk yield for 7 days (7AMY) was calculated. Rectal temperature (RT) was measured 3 times a day using a digital thermometer with a precision of 0.1 °C (MC-347, OMRON, Tokyo, Japan) by leaving the animal thermometer in the rectum for 10 s in the morning (0700–0900 h), afternoon (1400–1600 h), and evening (2100–2300 h) for 3 continuous days. At the same time, the breathing rate was measured by the flanking movement of the cow’s body for 30 s by visual inspection and multiplied by 2 to determine breaths per minute (breaths/min).

4.3. Blood Sample Collection and Total RNA Extraction

Blood samples (10 mL) for each cow were collected into ETDA2K anticoagulant tubes via tail vein puncture and centrifuged at 1400× g for 15 min at 4 °C. The middle leukocyte was extracted and stored at −80 °C for further RNA extraction. Following the protocol of the manufacturer, Trizol reagent (Invitrogen, Carlsbad, CA, USA) was employed for total RNA isolation and purification. RNA integrity was assessed using 1% agarose gel electrophoresis. The quality and quantity of RNA were assessed through NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA) and Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Only when RNA integrity number was >7.0 and RNA electrophoresis results (28S/18S ≥ 1.0), the RNAs could be used for further analysis.

4.4. Library Construction and RNA Sequencing

In each group, five cows were selected randomly from each group for RNA-seq. Sequencing libraries were generated following the manufacturer’s recommendation, and index codes were added to assign sequences to each sample. Afterwards, in accordance with the instructions of RiboZeroTM rRNA Removal Kit (Illumina, San Diego, CA, USA), ribosomal RNA (rRNA) was depleted. After the depletion of rRNA, the remaining RNAs were fragmented into small pieces and later reverse-transcribed to construct the cDNA library. Implementations include synthesizing U-labeled second-stranded DNAs with E. coli DNA polymerase I, RNase, and dUTP; adding A-base; ligating the adapter; and some other operations. Finally, llumina Hiseq 4000 (LC Bio, Hangzhou, Zhejiang, China) platform was used to perform the pair-end sequencing according to the vendor’s recommended protocol.

4.5. Circular RNA Identification and Differential Expression Analysis

Firstly, raw data obtained were filtered by Cutadapt [55], which could remove the reads that contained adaptor contamination, low-quality bases (Q ≥ 10), and undetermined bases. Then, in order to ensure the accuracy of subsequent analysis, sequence quality was evaluated by using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed on 8 February 2022)). After that, Bowtie2 [56] and Hisat2 [57] were applied to map reads to the Bos taurus reference genome. Later on, the remaining reads (58 unmapped reads) were extracted and mapped to genome using tophat-fusion [58]. With the aim of predicting circRNA more credibly, at the beginning, CIRCExplorer2 [59,60] and CIRI [61] were used to de novo assemble the mapped reads to circular RNAs. Then, tophat-fusion was used again to identify the back-splicing reads. Only the circRNA with at least one back-splicing read can be confirmed. To estimate the relative expression level of a circRNA, its quantity is denoted as spliced reads per billion mapped reads (SRPBM). The DE circRNAs were calculated between different treatment groups with |log2(fold change)| ≥ 1, p < 0.05 by R package edgeR [62].

4.6. GO and Pathway Analyses of Source Genes of DE circRNAs

The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted to explore the potential function of circRNAs. The gene corresponding to the parental mRNA of circRNA was identified with terms in the GO and KEGG database. In this manner, the function of gene was classified by GO pattern and KEGG pathway, and then all differentially expressed circRNAs could be mapped into terms of GO and KEGG. After comparing with the genome background using Fisher’s exact test, the significantly enriched (p < 0.05) GO and KEGG terms are kept for further analysis. OmicStudio tools were used to visualize the significantly changed heat-stress-related pathways and their enriched DE circRNAs in two comparisons.

4.7. Protein–Protein Interaction (PPI) Network Construction of Hub Source Genes of DE circRNAs

The String online database (https://cn.string-db.org/ (accessed on 20 March 2023)) provides some reliable information on PPI [63]. In this study, minimum required interaction score = 0.400 was considered. After extracting information on PPI, CytoNCA in Cytoscape was used for network visualization. Cytoscape is an open source software platform for visualizing complex networks and integrating these with any type of attribute data [64]. CytoNCA is a Cytoscape plugin for centrality analysis and evaluation of protein interaction networks [65] and is used to identify the hub source genes of DE circRNAs obtained from the PPI network.

4.8. Construction of circRNA–miRNA Network

Previous study hypothesized that by acting as miRNA sponge, circRNAs can bind miRNAs through an miRNA response element (MRE) and negatively regulate their activity [31]. The co-expressed circRNA and miRNA have potential as biomarkers and therapeutic targets. To uncover this co-expression network, we used miRanda [66] and Targetscan [67] to determine the interaction between circRNA and miRNA.

4.9. Quantitative Real-Time PCR Validation

Six circRNAs (four up-regulated and two down-regulated) were selected randomly to perform qRT-PCR to confirm the results of RNA-seq. In brief, 500 ng total RNA corresponding to the 15 samples was reverse-transcribed to cDNA. The expression levels of the 6 circRNAs were investigated by qPCR using Eastep R qPCR Master Mix (Promega, Shanghai, China), and the final result was observed by harnessing ABI Prism 7900HT sequence detection system (Thermo Fisher, Waltham, MA, USA). The primers were designed for RT-qPCR using circPrimer and Primer3. The circRNA and their corresponding primer sequences are provided in Table S11. Biological replicates were considered for all the quantitative PCR reactions. Afterward, using the designed primers and cDNA template, a polymerase chain reaction was conducted. The PCR conditions were 94 °C denaturation for 5 min, 40 cycles at 94 °C for 10 s, 54 to 60 °C for 15 s, and 72 °C for 30 s. The relative expression levels of circRNA were calculated using the 2−∆∆Ct method. The GAPDH was used as an internal control. Finally, agarose gel electrophoresis and DNA sequencing by Sanger sequencing were employed to confirm the qPCR products.

4.10. Statistical Analysis

The quantitative data were analyzed using SPSS software (ver. 24.0, SPSS Inc., Chicago, IL, USA), and graphs were generated using GraphPad Prism 7.0 software (ver. 5.03, GraphPad Software, San Diego, CA, USA). Data were expressed as mean ± standard error of the mean (s.e.m). Significant differences were determined by Student’s t-test or one-way ANOVA. p values < 0.05 were considered statistically significant differences unless otherwise stated.

5. Conclusions

In this study, two comparisons consisting of three experimental groups in spring and summer were designed to identify the effect of regulation of circRNA on heat stress. It can be seen that heat stress drastically reduced the overall milk yield and damaged the physiological condition of Holstein cows in Beijing in July. There were significant changes in circRNA expression profiles in peripheral blood samples from cattle under heat stress. Key pathways related to metabolic processes and autophagy such as the HIF-1 signaling pathway and longevity-regulating pathway were identified as participating in the cattle heat stress response. In addition, heat-stress-responsive circRNAs and their targeting miRNAs or corresponding genes were found. This valuable information will enrich our knowledge about involvement of key circRNAs and their expression patterns in cows’ heat stress processes. Moreover, our study can be further utilized in the identification, characterization, and designation of breeding strategies to develop both high-yielding and thermotolerant Holstein dairy cows.

Supplementary Materials

The following supporting information can be downloaded at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/ijms241210150/s1.

Author Contributions

Conceptualization and supervision, Q.X. and Y.W.; methodology and formal analysis, C.Z. and S.W.; writing—original draft preparation, C.Z.; resources and data curation, L.H. and H.F.; writing—review and editing, Y.Y.; visualization, G.C. and X.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Key Research and Development Projects (2022YFE0115700), China Agriculture Research System of MOF and MARA, the Key Research Project of Ningxia Hui Autonomous Region (2022BBF02017), and the Program for Changjiang Scholar and Innovation Research Team in University (IRT_15R62).

Institutional Review Board Statement

The experiment was performed in agreement with the Committee on Ethics of Animal Experimentation from the Beijing Jiaotong University, Beijing, China (Code ID: SS-QX-2014-06; 26 June 2014).

Informed Consent Statement

Not applicable.

Data Availability Statement

All the pertinent data are presented in the manuscript and associated supplementary files. Raw sequencing data can be obtained from the corresponding author.

Acknowledgments

We gratefully thank the reviewers for their careful reading and thoughtful suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Baker, J.S.; Murray, B.C.; McCarl, B.A.; Feng, S.Y.; Johansson, R. Implications of alternative agricultural productivity growth assumptions on land management, greenhouse gas emissions, and mitigation potential. Am. J. Agric. Econ. 2013, 95, 435–441. [Google Scholar] [CrossRef]
  2. Luo, H.; Hu, L.; Brito, L.F.; Dou, J.; Sammad, A.; Chang, Y.; Ma, L.; Guo, G.; Liu, L.; Zhai, L.; et al. Weighted single-step gwas and rna sequencing reveals key candidate genes associated with physiological indicators of heat stress in holstein cattle. J. Anim. Sci. Biotechnol. 2022, 13, 108. [Google Scholar] [CrossRef] [PubMed]
  3. Hu, L.; Sammad, A.; Zhang, C.; Brito, L.F.; Xu, Q.; Wang, Y. Transcriptome analyses reveal essential roles of alternative splicing regulation in heat-stressed holstein cows. Int. J. Mol. Sci. 2022, 23, 664. [Google Scholar] [CrossRef] [PubMed]
  4. Sammad, A.; Hu, L.R.; Luo, H.P.; Abbas, Z.; Umer, S.; Zhao, S.J.; Xu, Q.; Khan, A.; Wang, Y.J.; Zhu, H.B.; et al. Investigation of metabolome underlying the biological mechanisms of acute heat stressed granulosa cells. Int. J. Mol. Sci. 2022, 23, 2146. [Google Scholar] [CrossRef]
  5. Hu, L.R.; Brito, L.F.; Zhang, H.L.; Zhao, M.; Liu, H.Z.; Chai, H.; Wang, D.S.; Wu, H.J.; Cui, J.H.; Liu, A.R.; et al. Metabolome profiling of plasma reveals different metabolic responses to acute cold challenge between inner-mongolia sanhe and holstein cattle. J. Dairy Sci. 2022, 105, 9162–9178. [Google Scholar] [CrossRef] [PubMed]
  6. Czech, B.; Wang, Y.C.; Wang, K.; Luo, H.P.; Hu, L.R.; Szyda, J. Host transcriptome and microbiome interactions in holstein cattle under heat stress condition. Front. Microbiol. 2022, 13, 998093. [Google Scholar] [CrossRef]
  7. Czech, B.; Szyda, J.; Wang, K.; Luo, H.P.; Wang, Y.C. Fecal microbiota and their association with heat stress in bos taurus. BMC Microbiol. 2022, 22, 171. [Google Scholar] [CrossRef]
  8. Barrett, S.P.; Salzman, J. Circular rnas: Analysis, expression and potential functions. Development 2016, 143, 1838–1847. [Google Scholar] [CrossRef] [Green Version]
  9. Lee, E.C.S.; Elhassan, S.A.M.; Lim, G.P.L.; Kok, W.H.; Tan, S.W.; Leong, E.N.; Tan, S.H.; Chan, E.W.L.; Bhattamisra, S.K.; Rajendran, R.; et al. The roles of circular rnas in human development and diseases. Biomed. Pharmacother. 2019, 111, 198–208. [Google Scholar] [CrossRef]
  10. Hao, Z.Y.; Zhou, H.T.; Hickford, J.G.H.; Gong, H.; Wang, J.Q.; Hu, J.Q.; Liu, X.; Li, S.B.; Zhao, M.L.; Luo, Y.Z. Identification and characterization of circular rna in lactating mammary glands from two breeds of sheep with different milk production profiles using rna-seq. Genomics 2020, 112, 2186–2193. [Google Scholar] [CrossRef]
  11. Wang, L.D.; Liang, W.S.; Wang, S.S.; Wang, Z.X.; Bai, H.; Jiang, Y.; Bi, Y.L.; Chen, G.H.; Chang, G.B. Circular rna expression profiling reveals that circ-plxna1 functions in duck adipocyte differentiation. PLoS ONE 2020, 15, e0236069. [Google Scholar] [CrossRef] [PubMed]
  12. Yan, X.M.; Zhang, Z.; Meng, Y.; Li, H.B.; Gao, L.; Luo, D.; Jiang, H.; Gao, Y.; Yuan, B.; Zhang, J.B. Genome-wide identification and analysis of circular rnas differentially expressed in the longissimus dorsi between kazakh cattle and xinjiang brown cattle. PeerJ 2020, 8, 8646. [Google Scholar]
  13. Hong, L.J.; Gu, T.; He, Y.J.; Zhou, C.; Hu, Q.; Wang, X.W.; Zheng, E.Q.; Huang, S.X.; Xu, Z.; Yang, J.; et al. Genome-wide analysis of circular rnas mediated cerna regulation in porcine embryonic muscle development. Front. Cell Dev. Biol. 2019, 7, 289. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Li, H.; Yang, J.M.; Wei, X.F.; Song, C.C.; Dong, D.; Huang, Y.Z.; Lan, X.Y.; Plath, M.; Lei, C.Z.; Ma, Y.; et al. Circfut10 reduces proliferation and facilitates differentiation of myoblasts by sponging mir-133a. J. Cell. Physiol. 2018, 233, 4643–4651. [Google Scholar] [CrossRef]
  15. Wang, D.; Chen, Z.; Zhuang, X.; Luo, J.; Chen, T.; Xi, Q.; Zhang, Y.; Sun, J. Identification of circrna-associated-cerna networks involved in milk fat metabolism under heat stress. Int. J. Mol. Sci. 2020, 21, 4162. [Google Scholar] [CrossRef]
  16. Wang, D.; Zhao, Z.; Shi, Y.; Luo, J.; Chen, T.; Xi, Q.; Zhang, Y.; Sun, J. Circezh2 regulates milk fat metabolism through mir-378b sponge activity. Animal 2022, 12, 718. [Google Scholar] [CrossRef]
  17. Qi, Y.; Zhang, L.; Guo, Y.; Wang, J.; Chu, M.; Zhang, Y.; Guo, J.; Li, Q. Genome-wide identification and functional prediction of circular rnas in response to heat stress in chinese holstein cows. Anim. Biotechnol. 2021, 33, 1170–1180. [Google Scholar] [CrossRef]
  18. Zeng, H.F.; Xia, H.B.; Wang, X.L.; Wang, Y.; Fang, J.; Li, S.J.; Zhai, Y.F.; Han, Z.Y. Comprehensive profiling of cerna (circrna-mirna-mrna) networks in hypothalamic-pituitary-mammary gland axis of dairy cows under heat stress. Int. J. Mol. Sci. 2023, 24, 888. [Google Scholar] [CrossRef]
  19. Joshi, G.; Sharma, R.; Kakker, N.K. Phenotypic and functional characterization of t-cells and in vitro replication of fmdv serotypes in bovine lymphocytes. Vaccine 2009, 27, 6656–6661. [Google Scholar] [CrossRef]
  20. Gao, Y.H.; Li, J.B.; Cai, G.Z.; Wang, Y.J.; Yang, W.J.; Li, Y.Q.; Zhao, X.X.; Li, R.L.; Gao, Y.D.; Tuo, W.B.; et al. Single-cell transcriptomic and chromatin accessibility analyses of dairy cattle peripheral blood mononuclear cells and their responses to lipopolysaccharide. BMC Genom. 2022, 23, 338. [Google Scholar] [CrossRef]
  21. Liu, G.; Liao, Y.; Sun, B.; Guo, Y.; Deng, M.; Li, Y.; Liu, D. Effects of chronic heat stress on mrna and mirna expressions in dairy cows. Gene 2020, 742, 144550. [Google Scholar] [CrossRef]
  22. Dikmen, S.; Hansen, P.J. Is the temperature-humidity index the best indicator of heat stress in lactating dairy cows in a subtropical environment? J. Dairy Sci. 2009, 92, 109–116. [Google Scholar] [CrossRef] [Green Version]
  23. Bernabucci, U.; Lacetera, N.; Baumgard, L.H.; Rhoads, R.P.; Ronchi, B.; Nardone, A. Metabolic and hormonal acclimation to heat stress in domesticated ruminants. Animal 2010, 4, 1167–1183. [Google Scholar] [CrossRef] [Green Version]
  24. Berman, A.; Folman, Y.; Kaim, M.; Mamen, M.; Herz, Z.; Wolfenson, D.; Arieli, A.; Graber, Y. Upper critical-temperatures and forced ventilation effects for high-yielding dairy-cows in a sub-tropical climate. J. Dairy Sci. 1985, 68, 1488–1495. [Google Scholar] [CrossRef]
  25. Min, L.; Cheng, J.-b.; Shi, B.-l.; Yang, H.-j.; Zheng, N.; Wang, J.-q. Effects of heat stress on serum insulin, adipokines, amp-activated protein kinase, and heat shock signal molecules in dairy cows. J. Zhejiang Univ. Sci. B 2015, 16, 541–548. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Fan, C.Y.; Su, D.; Tian, H.; Hu, R.T.; Ran, L.; Yang, Y.; Su, Y.J.; Cheng, J.B. Milk production and composition and metabolic alterations in the mammary gland of heat-stressed lactating dairy cows. J. Integr. Agric. 2019, 18, 2844–2853. [Google Scholar] [CrossRef]
  27. Schutz, K.E.; Rogers, A.R.; Poulouin, Y.A.; Cox, N.R.; Tucker, C.B. The amount of shade influences the behavior and physiology of dairy cattle. J. Dairy Sci. 2010, 93, 125–133. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Tang, C.; Liang, Y.; Guo, J.H.; Wang, M.Q.; Li, M.X.; Zhang, H.M.; Arbab, A.A.I.; Karrow, N.A.; Yang, Z.P.; Mao, Y.J. Effects of seasonal heat stress during late gestation on growth performance, metabolic and immuno-endocrine parameters of calves. Animals 2022, 12, 716. [Google Scholar] [CrossRef]
  29. do Amaral, B.C.; Connor, E.E.; Tao, S.; Hayen, J.; Bubolz, J.; Dahl, G.E. Heat-stress abatement during the dry period: Does cooling improve transition into lactation? J. Dairy Sci. 2009, 92, 5988–5999. [Google Scholar] [CrossRef]
  30. Hsiao, K.Y.; Sun, H.S.; Tsai, S.J. Circular rna—New member of noncoding rna with novel functions. Exp. Biol. Med. 2017, 242, 1136–1141. [Google Scholar] [CrossRef] [Green Version]
  31. Memczak, S.; Jens, M.; Elefsinioti, A.; Torti, F.; Krueger, J.; Rybak, A.; Maier, L.; Mackowiak, S.D.; Gregersen, L.H.; Munschauer, M.; et al. Circular rnas are a large class of animal rnas with regulatory potency. Nature 2013, 495, 333–338. [Google Scholar] [CrossRef] [PubMed]
  32. Sang, M.J.; Wu, M.; Meng, L.J.; Zheng, Y.; Gu, L.N.; Liu, F.; Sang, M.X. Identification of epithelial-mesenchymal transition-related circrna-mirna-mrna cerna regulatory network in breast cancer. Pathol. Res. Pract. 2020, 216, 153088. [Google Scholar] [CrossRef] [PubMed]
  33. Xi, Y.; Fowdur, M.; Liu, Y.; Wu, H.; He, M.; Zhao, J. Differential expression and bioinformatics analysis of circrna in osteosarcoma. Biosci. Rep. 2019, 39, BSR20181514. [Google Scholar] [CrossRef] [Green Version]
  34. Xu, H.C.; Wang, C.Y.; Song, H.Y.; Xu, Y.X.; Ji, G. Rna-seq profiling of circular rnas in human colorectal cancer liver metastasis and the potential biomarkers. Mol. Cancer 2019, 18, 8. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Yao, Q.; He, Y.L.; Wang, N.; Dong, S.S.; Tu He Ta Mi Shi, M.E.; Feng, X.; Chen, H.; Pang, L.J.; Zou, H.; Zhou, W.H.; et al. Identification of potential genomic alterations and the circrna-mirna-mrna regulatory network in primary and recurrent synovial sarcomas. Front. Mol. Biosci. 2021, 8, 707151. [Google Scholar] [CrossRef]
  36. Dickson, R.C.; Lester, R.L. Sphingolipid functions in saccharomyces cerevisiae. Biochim. Biophys. Acta-Mol. Cell Biol. Lipids 2002, 1583, 13–25. [Google Scholar] [CrossRef]
  37. Dickson, R.C.; Nagiec, E.E.; Skrzypek, M.; Tillman, P.; Wells, G.B.; Lester, R.L. Sphingolipids are potential heat stress signals in saccharomyces. J. Biol. Chem. 1997, 272, 30196–30200. [Google Scholar] [CrossRef] [Green Version]
  38. Yang, M.; Kuang, M.Q.; Wang, G.L.; Ali, I.; Tang, Y.J.; Yang, C.X.; Li, Y.A.; Li, L. Choline attenuates heat stress-induced oxidative injury and apoptosis in bovine mammary epithelial cells by modulating perk/nrf-2 signaling pathway. Mol. Immunol. 2021, 135, 388–397. [Google Scholar] [CrossRef]
  39. Zhou, Z.; Trevisi, E.; Luchini, D.N.; Loor, J.J. Differences in liver functionality indexes in peripartal dairy cows fed rumen-protected methionine or choline are associated with performance, oxidative stress status, and plasma amino acid profiles. J. Dairy Sci. 2017, 100, 6720–6732. [Google Scholar] [CrossRef]
  40. Zhu, S.G.; Yang, Z.; Kong, L.L.; Kong, L.J.; Zhang, Y.Z. Arbutin inhibited heat stress-induced apoptosis and promoted proliferation and migration of heat-injured dermal fibroblasts and keratinocytes by activating pi3k/akt signaling pathway. Evid.-Based Complement. Altern. Med. 2022, 2022, 8798861. [Google Scholar] [CrossRef]
  41. Sies, H.; Jones, D.P. Reactive oxygen species (ros) as pleiotropic physiological signalling agents. Nat. Rev. Mol. Cell Biol. 2020, 21, 363–383. [Google Scholar] [CrossRef]
  42. Eom, H.J.; Ahn, J.M.; Kim, Y.; Choi, J. Hypoxia inducible factor-1 (hif-1)-flavin containing monooxygenase-2 (fmo-2) signaling acts in silver nanoparticles and silver ion toxicity in the nematode, caenorhabditis elegans. Toxicol. Appl. Pharmacol. 2013, 270, 106–113. [Google Scholar] [CrossRef] [PubMed]
  43. Baechler, B.L.; Bloemberg, D.; Quadrilatero, J. Mitophagy regulates mitochondrial network signaling, oxidative stress, and apoptosis during myoblast differentiation. Autophagy 2019, 15, 1606–1619. [Google Scholar] [CrossRef] [PubMed]
  44. Chen, B.; Yang, B.; Zhu, J.; Wu, J.; Sha, J.; Sun, J.; Bao, E.; Zhang, X. Hsp90 relieves heat stress-induced damage in mouse kidneys: Involvement of antiapoptotic pkm2-akt and autophagic hif-1alpha signaling. Int. J. Mol. Sci. 2020, 21, 1646. [Google Scholar] [CrossRef] [Green Version]
  45. Jiang, R.; Li, H.; Yang, J.; Shen, X.; Song, C.; Yang, Z.; Wang, X.; Huang, Y.; Lan, X.; Lei, C.; et al. Circrna profiling reveals an abundant circfut10 that promotes adipocyte proliferation and inhibits adipocyte differentiation via sponging let-7. Mol. Ther.-Nucleic Acids 2020, 20, 491–501. [Google Scholar] [CrossRef]
  46. Hardie, D.G.; Scott, J.W.; Pan, D.A.; Hudson, E.R. Management of cellular energy by the amp-activated protein kinase system. FEBS Lett. 2003, 546, 113–120. [Google Scholar] [CrossRef]
  47. Hardie, D.G. Amp-activated/snf1 protein kinases: Conserved guardians of cellular energy. Nat. Rev. Mol. Cell Biol. 2007, 8, 774–785. [Google Scholar] [CrossRef]
  48. Masroor, S.; Aalam, M.T.; Khan, O.; Tanuj, G.N.; Gandham, R.K.; Dhara, S.K.; Gupta, P.K.; Mishra, B.P.; Dutt, T.; Singh, G.; et al. Effect of acute heat shock on stress gene expression and DNA methylation in zebu (Bos indicus) and crossbred (Bos indicus × Bos taurus) dairy cattle. Int. J. Biometeorol. 2022, 66, 1797–1809. [Google Scholar] [CrossRef] [PubMed]
  49. Lee, J.; Lee, S.; Son, J.; Lim, H.; Kim, E.; Kim, D.; Ha, S.; Hur, T.; Lee, S.; Choi, I. Analysis of circulating-microrna expression in lactating holstein cows under summer heat stress. PLoS ONE 2020, 15, e0231125. [Google Scholar] [CrossRef]
  50. Raza, S.H.A.; Wijayanti, D.; Pant, S.D.; Abdelnour, S.A.; Hashem, N.M.; Amin, A.; Wani, A.K.; Prakash, A.; Dawood, M.A.O.; Zan, L. Exploring the physiological roles of circular rnas in livestock animals. Res. Vet. Sci. 2022, 152, 726–735. [Google Scholar] [CrossRef]
  51. Chen, X.; Dong, J.N.; Rong, J.Y.; Xiao, J.; Zhao, W.; Aschalew, N.D.; Zhang, X.F.; Wang, T.; Qin, G.X.; Sun, Z.; et al. Impact of heat stress on milk yield, antioxidative levels, and serum metabolites in primiparous and multiparous holstein cows. Trop. Anim. Health Prod. 2022, 54, 159. [Google Scholar] [CrossRef]
  52. Zhang, M.Q.; Jiang, M.; Bi, Y.; Zhu, H.; Zhou, Z.M.; Sha, J.H. Autophagy and apoptosis act as partners to induce germ cell death after heat stress in mice. PLoS ONE 2012, 7, e41412. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Maltz, E.; Kroll, O.; Barash, H.; Shamy, A.; Silanikove, N. Lactation and body weight of dairy cows: Interrelationships among heat stress, calving season and milk yield. J. Anim. Feed Sci. 2000, 9, 33–45. [Google Scholar] [CrossRef] [Green Version]
  54. McDowell, R.E.; Hooven, N.W.; Camoens, J.K. Effect of climate on preformance of holsteins in 1st laction. J. Dairy Sci. 1976, 59, 965–973. [Google Scholar] [CrossRef]
  55. Martin, M.J.E.j. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 2011, 17, 10–12. [Google Scholar] [CrossRef]
  56. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef] [Green Version]
  57. Kim, D.; Landmead, B.; Salzberg, S.L. Hisat: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef] [Green Version]
  58. Kim, D.; Salzberg, S.L. Tophat-fusion: An algorithm for discovery of novel fusion transcripts. Genome Biol. 2011, 12, R72. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Zhang, X.O.; Dong, R.; Zhang, Y.; Zhang, J.L.; Luo, Z.; Zhang, J.; Chen, L.L.; Yang, L. Diverse alternative back-splicing and alternative splicing landscape of circular rnas. Genome Res. 2016, 26, 1277–1287. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  60. Zhang, X.O.; Wang, H.B.; Zhang, Y.; Lu, X.H.; Chen, L.L.; Yang, L. Complementary sequence-mediated exon circularization. Cell 2014, 159, 134–147. [Google Scholar] [CrossRef] [Green Version]
  61. Gao, Y.; Wang, J.F.; Zhao, F.Q. Ciri: An efficient and unbiased algorithm for de novo circular rna identification. Genome Biol. 2015, 16, 4. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  62. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. Edger: A bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Szklarczyk, D.; Franceschini, A.; Kuhn, M.; Simonovic, M.; Roth, A.; Minguez, P.; Doerks, T.; Stark, M.; Muller, J.; Bork, P.; et al. The string database in 2011: Functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 2011, 39, D561–D568. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [Google Scholar] [CrossRef]
  65. Tang, Y.; Li, M.; Wang, J.; Pan, Y.; Wu, F.-X. Cytonca: A cytoscape plugin for centrality analysis and evaluation of protein interaction networks. Biosystems 2015, 127, 67–72. [Google Scholar] [CrossRef]
  66. Lewis, B.P.; Burge, C.B.; Bartel, D.P. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microrna targets. Cell 2005, 120, 15–20. [Google Scholar] [CrossRef] [Green Version]
  67. Enright, A.J.; John, B.; Gaul, U.; Tuschl, T.; Sander, C.; Marks, D.S. Microrna targets in drosophila. Genome Biol. 2004, 5, P8. [Google Scholar]
Ijms 24 10150 g001 550
Figure 1. Genomic characteristics of circRNA identified in Holstein cattle. (a) CircRNA numbers predicted in each group. (b) CircRNAs from exons, introns, or intergenic regions. (c) Distribution of circRNAs on chromosomes. (d) Number of exons contained within exonic circRNAs.
Figure 1. Genomic characteristics of circRNA identified in Holstein cattle. (a) CircRNA numbers predicted in each group. (b) CircRNAs from exons, introns, or intergenic regions. (c) Distribution of circRNAs on chromosomes. (d) Number of exons contained within exonic circRNAs.
Ijms 24 10150 g001
Ijms 24 10150 g002 550
Figure 2. DE circRNAs in response to heat stress in Holstein cattle. (a) Distribution of circRNA expression in 15 samples. (b,c) The volcano plot of DE circRNAs detected in two comparisons. Red dots represent the up-regulated circRNAs, blue dots represent the down-regulated circRNAs, and gray dots represent the unchanged circRNAs. (d) Venn diagram of DE circRNAs in comparison1 (Sum1 vs. Spr1) and comparison2 (Sum1 vs. Spr2) data sets. (e,f) The heatmaps of DE circRNAs in two comparisons.
Figure 2. DE circRNAs in response to heat stress in Holstein cattle. (a) Distribution of circRNA expression in 15 samples. (b,c) The volcano plot of DE circRNAs detected in two comparisons. Red dots represent the up-regulated circRNAs, blue dots represent the down-regulated circRNAs, and gray dots represent the unchanged circRNAs. (d) Venn diagram of DE circRNAs in comparison1 (Sum1 vs. Spr1) and comparison2 (Sum1 vs. Spr2) data sets. (e,f) The heatmaps of DE circRNAs in two comparisons.
Ijms 24 10150 g002
Ijms 24 10150 g003 550
Figure 3. Annotations and enrichment analysis of the DE circRNAs in Sum1 and Spr1 groups. (a) GO annotations of DE circRNA host genes. Diagram shows the significantly enriched GO terms (p < 0.05). The abscissa represents the terms in GO categories of biological process (BP), cellular component (CC), and molecular function (MF). (b) KEGG enrichment analysis of DE circRNAs. The abscissa indicates the rich factor, and the ordinate indicates the name of the enriched pathway. The size of the dot represents the number of enriched genes in the pathway, and the color corresponds to the p value as indicated in the figure legend.
Figure 3. Annotations and enrichment analysis of the DE circRNAs in Sum1 and Spr1 groups. (a) GO annotations of DE circRNA host genes. Diagram shows the significantly enriched GO terms (p < 0.05). The abscissa represents the terms in GO categories of biological process (BP), cellular component (CC), and molecular function (MF). (b) KEGG enrichment analysis of DE circRNAs. The abscissa indicates the rich factor, and the ordinate indicates the name of the enriched pathway. The size of the dot represents the number of enriched genes in the pathway, and the color corresponds to the p value as indicated in the figure legend.
Ijms 24 10150 g003
Ijms 24 10150 g004 550
Figure 4. Annotations and enrichment analysis of the DE circRNAs in Sum1 vs. Spr2 groups. (a) GO annotations of DE circRNA host genes. Diagram shows the significantly enriched GO terms (p < 0.05). The abscissa represents the terms in GO categories of biological process (BP), cellular component (CC), and molecular function (MF). (b) KEGG enrichment analysis of DE circRNAs. The abscissa indicates the rich factor, and the ordinate indicates the name of the enriched pathway. The size of the dot represents the number of enriched genes in the pathway, and the color corresponds to the p value as indicated in the figure legend.
Figure 4. Annotations and enrichment analysis of the DE circRNAs in Sum1 vs. Spr2 groups. (a) GO annotations of DE circRNA host genes. Diagram shows the significantly enriched GO terms (p < 0.05). The abscissa represents the terms in GO categories of biological process (BP), cellular component (CC), and molecular function (MF). (b) KEGG enrichment analysis of DE circRNAs. The abscissa indicates the rich factor, and the ordinate indicates the name of the enriched pathway. The size of the dot represents the number of enriched genes in the pathway, and the color corresponds to the p value as indicated in the figure legend.
Ijms 24 10150 g004
Ijms 24 10150 g005 550
Figure 5. Relationship between significant pathways and enriched DE circRNAs in two comparisons. The red square nodes represent comparison1 (Sum1 vs. Spr1) and comparison2 (Sum1 vs. Spr2), the green circular nodes represent pathways, the purple triangle nodes represent DE circRNAs, the solid lines represent significantly enriched pathways in comparisons, and the dotted lines represent DE circRNAs involved in enriched pathways.
Figure 5. Relationship between significant pathways and enriched DE circRNAs in two comparisons. The red square nodes represent comparison1 (Sum1 vs. Spr1) and comparison2 (Sum1 vs. Spr2), the green circular nodes represent pathways, the purple triangle nodes represent DE circRNAs, the solid lines represent significantly enriched pathways in comparisons, and the dotted lines represent DE circRNAs involved in enriched pathways.
Ijms 24 10150 g005
Ijms 24 10150 g006 550
Figure 6. Identification of hub source genes of circRNA from the PPI network. (a) PPI network of 74 source genes of circRNA. (b) PPI network of 26 source hub genes of circRNA that were extracted from the PPI network. The larger the circle, the higher the Betweenness Centrality (BC); conversely, the smaller the circle, the lower the BC. The red circles note the top 10 hub source genes of circRNA.
Figure 6. Identification of hub source genes of circRNA from the PPI network. (a) PPI network of 74 source genes of circRNA. (b) PPI network of 26 source hub genes of circRNA that were extracted from the PPI network. The larger the circle, the higher the Betweenness Centrality (BC); conversely, the smaller the circle, the lower the BC. The red circles note the top 10 hub source genes of circRNA.
Ijms 24 10150 g006
Ijms 24 10150 g007 550
Figure 7. Interaction between circRNAs and miRNAs. The red squares represent circRNA, the purple circles represent miRNA, and the lines represent the relationships between circRNAs and miRNAs.
Figure 7. Interaction between circRNAs and miRNAs. The red squares represent circRNA, the purple circles represent miRNA, and the lines represent the relationships between circRNAs and miRNAs.
Ijms 24 10150 g007
Ijms 24 10150 g008 550
Figure 8. Validation of DE circRNAs. (a) The gel electrophoresis results of circRNAs. (b) Sanger sequencing confirmed the back-splice junction sequence of the indicated circRNAs. (c) Comparison of results from RNA-seq and qRT-PCR. Six DE circRNAs were randomly selected, and their expression level was validated by qRT-PCR.
Figure 8. Validation of DE circRNAs. (a) The gel electrophoresis results of circRNAs. (b) Sanger sequencing confirmed the back-splice junction sequence of the indicated circRNAs. (c) Comparison of results from RNA-seq and qRT-PCR. Six DE circRNAs were randomly selected, and their expression level was validated by qRT-PCR.
Ijms 24 10150 g008
Table
Table 1. The performances and physiological indicators of Holstein cattle in Spr1, Spr2, and Sum1 groups.
Table 1. The performances and physiological indicators of Holstein cattle in Spr1, Spr2, and Sum1 groups.
VariableSprSump-Value 3
Spr1
(MEAN ± SD) 1		Spr2
(MEAN ± SD) 1		Sum1
(MEAN ± SD) 1
Individual151515 2NA
DIM (d)61.6 ± 10.62137.53 ± 7.63137.60 ± 10.62NA
7AMY (Kg/day)40.91 ± 4.40 ab42.93 ± 3.25 a38.99 ± 3.64 b<0.05
MY (Kg)43.95 ± 3.89 a40.88 ± 2.05 b37.51 ± 5.36 c<0.05
RT (°C)38.63 ± 0.18 b38.58 ± 0.13 b38.96 ± 0.24 a<0.05
RR (breaths/min)38.80 ± 3.63 b37.60 ± 5.14 b88.4 ± 12.39 a<0.05
THI54.96 ± 3.2781.08 ± 4.57<0.05
Note: DIM, days in milk; 7AMY, average milk yield across seven days; MY, milk yield on the blood sampling day; RT, rectal temperature; RR, respiration rate; THI, temperature–humidity index; 1 Mean ± standard deviation; 2 they are the same individuals as in Spr1; a,b,c Different letters in the same row indicate significant differences between the two groups (p < 0.05); 3 p < 0.05 indicates a significant difference among groups; NA: not applicable.
Table
Table 2. The GO term of overlapped circRNAs in two comparisons.
Table 2. The GO term of overlapped circRNAs in two comparisons.
GO TermCircRNA (Source Gene)p-Value 1
Guanyl-nucleotide exchange factor activityCircRNA27860 (ANKRD27); circRNA27332
(RALGPS2); circRNA3628 (MCF2L2); circRNA27183 (DOCK7)		0.001
Cellular response to transforming growth factor beta stimulusCiRNA3301 (WWOX); circRNA4205 (NR3C1)0.045
Small GTPase-mediated signal transductionCircRNA27332 (RALGPS2); circRNA13493 (RHEB); circRNA27183 (DOCK7)0.001
Repressing transcription factor bindingCiRNA292 (RUNX2); circRNA4205 (NR3C1)0.001
CytosolcircRNA27437 (PPHLN1); circRNA4205 (NR3C1); circRNA11103 (STAM2); circRNA13493 (RHEB); ciRNA3301 (WWOX)0.001
Note: Words in bold indicate the DE circRNAs were enriched in multiple GO and KEGG entries. 1 p-values were calculated using Student’s t-test.
Table
Table 3. The KEGG term of overlapped circRNAs in two comparisons.
Table 3. The KEGG term of overlapped circRNAs in two comparisons.
KEGG TermCircRNA (Source Gene)p-Value 1
Choline metabolism in cancerCircRNA13493 (RHEB); circRNA5677 (GPCPD1)0.005
EndocytosisCircRNA11103 (STAM2); circRNA27501 (ZFYVE16)0.025
Longevity-regulating pathwayCircRNA13493 (RHEB)0.0057
Note: Words in bold indicate the DE circRNAs were enriched in multiple GO and KEGG entries. 1 p-values were calculated using Student’s t-test.
Table
Table 4. The top 10 hub source genes of DE circRNA in 2 comparisons.
Table 4. The top 10 hub source genes of DE circRNA in 2 comparisons.
AccessionSource GeneSum1 vs. Spr1Sum1 vs. Spr2
Regulationp-Value 1Regulationp-Value 1
CircRNA1209PTENdown0.01down0.15
CiRNA1282HIF1Adown0.45down0.01
CircRNA14021PTPRCup0.04up0.47
CircRNA13493RHEBdown0.01down0.04
CircRNA22544EZH2up0.03up0.40
CircRNA292RUNX2down0.02down0.00
CircRNA2467NR3C1down0.44down0.05
CircRNA4205NR3C1down0.01down0.00
CircRNA3691PIK3CBdown0.02down0.09
CircRNA12923ROCK1down0.01down0.39
CircRNA5505PLD1down0.01down0.39
Note: 1 p-values were calculated using Student’s t-test.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, C.; Wang, S.; Hu, L.; Fang, H.; Chen, G.; Ma, X.; Yu, Y.; Wang, Y.; Xu, Q. Analysis of CircRNA Expression in Peripheral Blood of Holstein Cows in Response to Heat Stress. Int. J. Mol. Sci. 2023, 24, 10150. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms241210150

AMA Style

Zhang C, Wang S, Hu L, Fang H, Chen G, Ma X, Yu Y, Wang Y, Xu Q. Analysis of CircRNA Expression in Peripheral Blood of Holstein Cows in Response to Heat Stress. International Journal of Molecular Sciences. 2023; 24(12):10150. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms241210150

Chicago/Turabian Style

Zhang, Congcong, Shuhui Wang, Lirong Hu, Hao Fang, Gong Chen, Xiaojuan Ma, Ying Yu, Yachun Wang, and Qing Xu. 2023. "Analysis of CircRNA Expression in Peripheral Blood of Holstein Cows in Response to Heat Stress" International Journal of Molecular Sciences 24, no. 12: 10150. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms241210150

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop