Small RNA-Seq Analysis Reveals miRNA Expression of Short Distance Transportation Stress in Beef Cattle Blood

Simple Summary In this study, three miRNA modules were identified in a cattle short-distance transportation stress model, and the turquoise module showed key miRNA sets according to their correlation with hub genes. Further, hub miRNAs were identified based on their targeting relationship with hub genes in our previous study. This finding provides the potential utility for predicting and treatment of short-distance transportation stress in beef cattle. Abstract Transportation is a crucial phase in the beef cattle industry, and the annual losses caused by beef cattle transport stress are substantial. Because of its huge economic losses, such as lower growth rate and even death, long-distance transportation stress has attracted more attention from beef production practitioners because of its huge economic losses. Compared with the long-distance transportation stress, the short-distance transportation stress was ignored for the reason of no obvious symptoms in cattle. Our previous study showed that the disorder of B cell function could be a potential health risk after short-distance transportation. However, the transcriptome details of the changes in the cattle blood after short-distance transportation and the molecular mechanisms for the regulation of the developmental process are not clearly known. In this study, a total of 10 Qinchuan cattle were used to compare the molecular characteristics of blood before and after short-distance transportation. The miRNA-seq showed that 114 differentially expressed miRNAs (DEMs) were found (40 upregulated and 74 downregulated) between two groups before and after transportation. Furthermore, more than 90% of the miRNAs with counts of more than 10 were used to construct a co-expression network by weighted correlation network analysis (WGCNA), and four independent modules were identified. According to their relationship with 30 hub genes, the turquoise module was the key module in this study. The regulator network of hub genes and miRNAs in the turquoise module was constructed by miRNAs targeting genes predicting, and the miRNAs had targeting sites within hub genes that could be identified as hub-miRNAs. Further, it showed that CD40 and ITPKB had the same targeting miRNAs (miR-339a/b), and the newly discovered hub miRNAs filled the gaps in our previous study about the relationship between hub genes in short-distance transportation stress and provided the potential utility for predicting and treatment of short-distance transportation stress in beef cattle.


Introduction
With the improvement of social productivity, animal welfare issues have become increasingly prominent. Animal welfare not only affects the production efficiency of animal husbandry, but also affects human health seriously. Transportation stress is an important link in animal production, involving the economic benefit and welfare of animal production [1]. Transportation stress can be defined as the change of environment significantly that disturbs body homeostasis. Compared with the transportation stress process of other livestock, beef cattle are paid more attention because of the particularity of their digestive system, especially the rumen digestive function. Transportation stress has brought huge economic losses to the beef cattle industry, which could be divided into long-distance transportation stress and short-distance transportation stress [2,3]. On the one hand, long-distance transportation stress has attracted more attention from beef production practitioners, because it can slow down the growth rate and even cause death of beef cattle [4,5]. The mechanism of long-distance transportation was interpreted, including hormone, amino acids, glucose, lipid, cholesterol, immunocyte, gastrointestinal microbiota, nasopharyngeal microbiota, and carcass quality [6][7][8][9][10]. On the other hand, most cattle were normal apparently after short-distance transportation, but could have potential health risks. The mechanism of short-distance transportation was still poorly understood compared with long-distance transportation. Our previous study showed that the disorder of B cell differentiation, proliferation, survival, and apoptosis were the potential molecular mechanisms in short-distance transportation stress [11]. It suggested that beef cattle are susceptible to the attack of pathogenic bacteria after short-distance transportation. Considering the production efficiency and animal welfare of beef cattle, attention should be paid to the transportation stress of short-distance transport of beef cattle. According to the results of transcriptome analysis, three potential molecular markers genes were identified in 30 hub genes, but the regulatory network of hub genes is still unknown. Noncoding RNAs participate in the restoration of cellular homeostasis or adaptation to environmental conditions through changes in the gene expression programs. As an important factor of noncoding RNAs, microRNAs (miRNAs) play key roles in the regulation of cellular homeostasis in eukaryotic organisms. While the uses of miRNA in the diagnosis of transportation stress and its mitigations are in their early phases, successful miRNA-based therapeutics have been established in essential exploration. The transcriptome details of what changes occur in the cattle blood after short-distance transportation and the molecular mechanisms for the regulation of the developmental process are not clearly known.
This study aimed to explore the change in blood transcriptome of miRNAs after short-distance transportation, thereby updating the theoretical basis for the diagnosis of Beef cattle transport stress syndrome (TSSBC) and providing targeting sites of potential miRNA-based therapeutics for the management of beef cattle production and safeguarding of animal welfare.

Animal Model of Short-Distance Transport Stress
Ten healthy, unrelated female beef cattle (3~4 years old) were used to construct the short-distance transportation stress model. Before transportation, the cattle were kept in loose housing conditions and fed on a total mixed ration (TMR). The beef cattle were deprived of food and water that day until after transportation. Blood samples were collected from each cattle before transportation, marked group B. After transportation, blood samples were collected immediately and marked group A. Group B was control group. The transportation route was: from Qinbao Cattle Industry Co., Ltd in Yangling county, Shaanxi province to Qinbao Cattle Industry Co., Ltd in Qishan county, Shaanxi province. The average transportation density was 1.28 m 2 /head, which was similar to the production process. The transportation distance was 70 km, and the vehicle was a singlelayer Foton cart with a maximum speed of 60 km/h and an average speed of 30 km/h. The interval time between the two blood samples collection was 6-8 h. The temperature ranged from 26~32 • C, and the humidity was 70% [11].

RNA Extraction, sRNA Library Construction and Sequencing
Twenty blood samples were collected in sodium heparin anticoagulant containing tubes. ACK lysis buffer lysed erythrocytes, karyocytes were harvested following centrifugation at 4000 rpm for 6 min at 4 • C, and then transferred into Eppendorf tubes containing TRIzol, and finally stored in liquid nitrogen. Total RNA was isolated using the RNAiso Plus kit (Takara, Tokyo, Japan). RNA concentration and quality were evaluated using Qubit2.0 RNA detection kit (Life, California, USA) and nanodrop 1000 (Thermo, Massachusetts, USA). Considering the quality of RNA, fourteen miRNA cDNA libraries were built in this study. The steps of miRNA cDNA library construction were the following: Ligated with 3 and 5 adapters (Illumina, San Diego, CA, USA) using T4 ligase (New England Biolabs), purified RNA was reverse-transcribed into the first strand cDNA and amplified by PCR using primers complementary to the adaptor sequences. The final small RNA sequencing library was prepared by purified the nucleotide fractions at 140~150 bp length. After then, each library was loaded into a single Illumina Hiseq (Illumina, San Diego, USA) lane with 75 bp single-end sequencing.

Primer Design and qPCR
The primers of miRNA were designed by tailing reaction on Sangon online software (https://www.sangon.com/newPrimerDesign, accessed on 20 August 2020) ( Table 1). miRNA First Strand cDNA Synthesis (Tailing Reaction) (Sangon, Shanghai, China) was used for reverse transcription of RNA to cDNA. qPCR was carried out in Y480 Real-Time PCR Detection System (Roche, Basel, Switzerland) utilizing SYBR green detection (Takara Bio, Tokyo, Japan). The amplification protocol was as follows: 95 • C for 30 s, followed by 50 cycles of 95 • C 10 s, and 60 • C for 30 s. Melt curve analysis was performed between 55 and 95 • C, with a 0.5 • C increment every 5 s. Samples were run in triplicate. Forward primers were listed in Table 1 and reverse primers were provided in the cDNA Synthesis kit. U6 was used as the reference gene, its primers were kept secret and also provided in cDNA Synthesis kit. All expression levels were normalized to that of U6 and quantified using the 2-∆∆Ct method [12]. Twelve differentially expressed miRNA (DEMs) (fold-change > 2 and false discovery rate (FDR) < 0.05) were randomly chosen for verification by qPCR. Table 1. List of the primers used in this study.

Preliminary Analysis of RNA-Seq Data and Verification by qPCR
Fourteen libraries representing seven animals in each group (before and after transportation) were prepared from total leukocyte miRNA, and group B was the reference group. The range of total reads counts for the 14 samples were 11,003,312~18,920,754 ( Table 2). The range of clean reads counts for the 14 samples were 9,269,728~16,096,861 ( Table 2). The range of mapped ratio for the 14 samples were 79.41~89.20 ( Table 2). The range of average reads length were 22.50~24.30, most of reads distributed from 19 to 25 bp ( Figure 1A). For each sample, the Q20 base ratio was >97.98%, and the Q30 base ratio was >96.82% ( Table 2). The percentage of miRNA in total reads was between 21.84~41.12% ( Figure 1B). Known miRNA were more than 90% of all miRNAs, the rest were novel miRNA ( Figure 1C,D). Principal component analysis (PCA) showed that 14 samples were classified into two groups (group B and group A; Figure 2A). There were 114 DEMs (fold-change > 2 and FDR < 0.05) between group B and group A (Table S1), including 40 upregulated miR-NAs and 74 downregulated miRNAs ( Figure 2B). All the DEMs were shown in the heatmap ( Figure 2C). Twelve DEMs (fold-change > 2 and FDR < 0.05) were randomly chosen for verification by qPCR. All the DEMs showed the same trend as the RNA-seq data ( Figure 3).

Enrichment of DEMs between before and after Transportation
To identify key differences between before and after transportation, GO and K  Figure 4C). The KEGG analys showed 104 pathways were significantly enriched (FDR < 0.05). The pathways we mainly focal adhesion, axon guidance, ECM-receptor interaction, etc. ( Figure 4D).

Co-Expression Analysis and Their Relationship with Hub Genes
All the miRNAs with their counts greater than 10 in more than 90% samples we used to construct a co-expression network. After the samples cluster, all the samples we used to construct a co-expression network ( Figure 5A). The first power number when co rection index more than 0.85 were selected for the next analysis (power = 6, Figure 5B When power = 6, the mean connectivity was less than 10 and suited for the next analys ( Figure 5C). According to the parameters list in 2.4 RNA-seq analysis and statistics, thr modules (blue, brown, and turquoise) were identified in 14 samples, which respective have 98, 42, and 235 miRNAs within (Table S4 and Figure 6A). The results of netwo heatmap indicated a significant difference among modules ( Figure 5D-F). The correlatio of hub genes and miRNA modules showed that the turquoise module was the key modu in the short-distance transportation stress model ( Figure 6B). According to the relatio ship between miRNAs in turquoise module and hub genes in our previous study, w found CD40 and ITPKB have the same targeting miRNAs ( Figure 6C).

Enrichment of DEMs between before and after Transportation
To identify key differences between before and after transportation, GO and KEGG enrichment were performed to determine DEMs' function (Tables S2 and S3). The biological process (BP) of GO terms showed DEMs those are involved in intracellular signal transduction, protein phosphorylation, cell adhesion, etc. ( Figure 4A). The cellular component (CC) of GO terms showed DEMs those are involved in receptor complex, basement membrane, focal adhesion, etc. ( Figure 4B). The molecular function (MF) of GO terms showed DEMs those are involved in protein kinase activity, ATP binding, transcription factor activity, sequence-specific DNA binding, etc. ( Figure 4C). The KEGG analysis showed 104 pathways were significantly enriched (FDR < 0.05). The pathways were mainly focal adhesion, axon guidance, ECM-receptor interaction, etc. ( Figure 4D).

Co-Expression Analysis and Their Relationship with Hub Genes
All the miRNAs with their counts greater than 10 in more than 90% samples were used to construct a co-expression network. After the samples cluster, all the samples were used to construct a co-expression network ( Figure 5A). The first power number when correction index more than 0.85 were selected for the next analysis (power = 6, Figure 5B). When power = 6, the mean connectivity was less than 10 and suited for the next analysis ( Figure 5C). According to the parameters list in 2.4 RNA-seq analysis and statistics, three modules (blue, brown, and turquoise) were identified in 14 samples, which respectively have 98, 42, and 235 miRNAs within (Table S4 and Figure 6A). The results of network heatmap indicated a significant difference among modules ( Figure 5D-F). The correlation of hub genes and miRNA modules showed that the turquoise module was the key module in the short-distance transportation stress model ( Figure 6B). According to the relationship between miRNAs in turquoise module and hub genes in our previous study, we found CD40 and ITPKB have the same targeting miRNAs ( Figure 6C).

Discussion
Transportation stress is closely related to animal welfare and animal productivity. Transportation stress is essentially a phenomenon of shock to the homeostasis of the organism. It was found that cortisol levels were higher after the transportation than before in cattle [1]. In order to reduce the damage to animals caused by transport stress, it is of great significance to construct transcriptional regulatory networks based on transport stress. miRNA may be involved in post-transcriptional regulation and can explain part of the gene expression that cannot be explained by transcriptional level. Our previous study did not reveal the hematological molecular mechanism of transportation stress fully, revealing the expression pattern of miRNA after transportation could elucidate the internal relationships between hub genes and enrich the regulation mechanism of transcriptional transport stress was enriched at the epigenetic level. Because there was no obvious symptom after short-distance transportation, it had not been sufficiently touched upon as they are supposed to. Based on the model of short-distance transportation stress in beef cattle, the purpose of this study was to reveal the role of miRNA in short-distance transportation stress in cattle blood, in order to provide a theoretical basis and potential therapeutic target for the prevention and treatment of beef cattle short-distance transportation stress.

Discussion
Transportation stress is closely related to animal welfare and animal productivity. Transportation stress is essentially a phenomenon of shock to the homeostasis of the organism. It was found that cortisol levels were higher after the transportation than before in cattle [1]. In order to reduce the damage to animals caused by transport stress, it is of great significance to construct transcriptional regulatory networks based on transport stress. miRNA may be involved in post-transcriptional regulation and can explain part of the gene expression that cannot be explained by transcriptional level. Our previous study did not reveal the hematological molecular mechanism of transportation stress fully, revealing the expression pattern of miRNA after transportation could elucidate the internal relationships between hub genes and enrich the regulation mechanism of transcriptional transport stress was enriched at the epigenetic level. Because there was no obvious symptom after shortdistance transportation, it had not been sufficiently touched upon as they are supposed to. Based on the model of short-distance transportation stress in beef cattle, the purpose of this study was to reveal the role of miRNA in short-distance transportation stress in cattle blood, in order to provide a theoretical basis and potential therapeutic target for the prevention and treatment of beef cattle short-distance transportation stress.
The results of Figure 1 showed that the reads length was enriched from 19~25 bp and 21%~41% of the reads were miRNAs. Most of reads were known miRNAs and part of novel miRNAs were predicted in this study. Through qPCR verification, 12 randomly selected miRNAs have the same tendency as miRNA sequencing, which showed sequencing data has certain credibility in this study. GO and KEGG enrichment were used to identify the function of miRNA between before and after transportation, and the pathways related to energy supply were enriched, such as the insulin signaling pathway, type II diabetes mellitus, and sphingolipid signaling pathway. It indicated that the energy supply was a key factor in short-distance transportation stress, and confirmed the results of our previous study about mRNA expression pattern in short-distance transportation stress. The reason for transportation stress was complex, including temperature, wind, rain, hunger, thirst, crowding, shock, turbulence, social interaction, feed change, physical exertion, environmental change, and invasion of pathogenic bacteria, and part of miRNA were confirmed in recent studies [18][19][20][21]. miRNA sequencing in chronic heat stress showed that miR-199b was highly expressed and had the same expression tendency in this study [22]. miR-2284 families were also verified in circulating-miRNA expression in lactating Holstein cows under summer heat stress [23]. It suggested that miR-199b and miR-2284 families were related to heat stress response. The conditions of transportation stress were unique and unrepeatable, and the molecular mechanism of transportation stress was incomplete, especially in short-distance transportation stress. Numerous studies have begun to address the details of transportation stress, such as fasting and management [24,25]. Their focus was on hormone and mRNA expression, and the part of miRNA expression was lacking. The transcription factor NF-κB is involved in the inflammatory response and the process of immunity [26]. It has been reported that there is a direct link between stress-induced NF-κB activation and the upregulation of some miRNA [27]. It suggested that miRNA plays an important role in immune response during stress. In addition, studies have shown that stress conditions can change the function of miRNA and the expression of mRNA targets. In addition to the expression levels of miRNA and mRNA targets, the function of miRNA depends on the miRNA that has a similar target near the target [28]. It was with regret that no hub miRNA was identified through GO and KEGG enrichment analysis, and then WGCNA was used to identify hub miRNAs in short-distance transportation stress. Three modules were identified in the short-distance transportation stress model. Based on their correlation of hub genes set in our previous study, the turquoise module was considered as a key module. Combined with the targeting relationship between hub genes and miRNAs in turquoise module, nine miRNAs were considered as hub miRNAs in short-distance transportation stress. miR-339a, miR-339b, miR-538-mature, and novel-654 mature all have target sites of CD40 and ITPKB gene, and this revealed the regulatory network of hub genes at the epigenetic level. The CD40 was a specific B cell surface marker and was associated with the differentiation of B cells to memory B cells, and ITPKB have an important essential function of survival of naïve mature B cells, which indicated that the miRNA could affect the immune response of cattle by CD40 and ITPKB in transport stress [29][30][31][32][33]. These newly discovered hub miRNAs fill the gaps in our previous study about the relationship between hub genes in short-distance transportation stress. Despite that limitation, we believe our findings have the potential utility for predicting and treatment of short-distance transportation stress. Interestingly, previous studies have also identified the miRNAs that we recognized as influential on short-distance transportation stress. Further investigations with larger numbers and different conditions of transportation stress models are appropriate to validate the application of using these biomarkers and the proposed scoring mechanism.

Conclusions
In this study, 10 Qinchuan cattle were used to compare the changes of blood molecular characteristics before and after in short-distance transportation. One hundred and fourteen differentially expressed miRNAs were identified by miRNA-seq and four key modules were identified by WGCNA. Combined with the targeting relationship between miRNA and hub gene in previous research, we further screened four significant miRNAs, which provided a certain basis for the diagnosis of short-distance transport stress.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ani11102850/s1, Table S1. The differential expression miRNAs (DEMs) list. Table S2. The GO analysis of the DEMs. Table S3. The KEGG analysis of the DEMs. Table S4

Data Availability Statement:
The high-throughput sequencing data of the small RNA-seq have been saved in the NCBI Sequence Reading Archive (https://www.ncbi.nlm.nih.gov/sra/), with the accession number PRJNA764198 (SRR15944129-SRR15944142).