Transcriptome Characterization of Short Distance Transport Stress in Beef Cattle Blood

The transportation is a crucial phase in beef cattle industry, and the annual losses caused by beef cattle transport stress are substantial. Several studies have described the effect of long distance transportation stress on animal health, such as disorder in nervous, endocrine, immune, and metabolic system. However, molecular mechanisms underlying short distance transportation stress is still poorly understood. Present study aims to investigate the effect of short distance transportation by measuring the hematological indices and transcriptomic analysis. In this study, a total 10 Qinchuan cattle were used to compare the molecular characteristics of blood before and after transportation. We have found that a stress-related marker “white blood cell count (WBC)” increased significantly after transportation. The decrease in triglyceride (TG), cholestenone (CHO), high-density lipoprotein (HDL), and low-density lipoprotein (LDL) showed that energy expenditure was increased after transportation, but not enough to activate fatty decomposition. Intriguingly, the decrease of malondialdehyde (MDA) showed that cattle were more resilience to oxidative stress. The RNA-seq showed that 1,092 differentially expressed genes (DEGs) were found (329 up-regulated and 763 down-regulated) between group before and group after. The GO and KEGG enrichment showed that the metabolic pathway and B cell function related pathways were enriched. Furthermore, median absolute deviation (MAD) top 5,000 genes were used to construct a co-expression network by weighted correlation network analysis (WGCNA), and 11 independent modules were identified. Combing with protein-protein interaction (PPI) analysis, the verification of quantitative real-time PCR (qPCR) and the correlation of B cell function, structural maintenance of chromosomes 3 (SMC3), jun proto-oncogene (JUN), and C-X-C motif chemokine ligand 10 (CXCL10) were suggested as potential molecular markers in identification of short distance transportation. Collectively, the blood RNA-seq analysis and WGCNA indicated that the disorder of B cell differentiation, proliferation, survival, and apoptosis were the potential molecular mechanism in short distance transportation stress. In conclusion, our results provide the novel insight about potential biomarkers for short distance transportation stress, which may serve as for diagnosing and preventing this condition in beef industry.


INTRODUCTION
Freshly cut meat has become a major product in the meat market worldwide. Livestock is bought and sold more frequently in the market economy and transportation is an integral part of the livestock industry. The key factors that can affect the welfare of animals and the economic benefits of this industry during transportation include a higher susceptibility to disease and increased disease spread (Schwartzkopf-Genswein et al., 2012;Deng et al., 2017). Transportation stress syndrome (TSS) often appears in beef cattle after ground transportation, and can affect many tissues, including the nervous, endocrine, immune, and energy supply systems (Van Engen and Coetzee, 2018).
The stressors responsible for inducing transport stress syndrome are complex, including temperature, wind, rain, hunger, thirst, crowding, shock, turbulence, social interaction, feed change, physical exertion, environmental change, and invasion of pathogenic bacteria (Delic et al., 2014). Corticosterone secretion and lipid mobilization are normal phenomena under stress conditions (Chirase et al., 2004;Kinlein et al., 2019). Although TSS gradually subside through vitamin in feed, however, in some cases, a few animals do not recover by vitamins supplement, eventually decrease the production efficiency or die (Kegley et al., 2016;Takemoto et al., 2017;Takemoto et al., 2018). Greater attention should be paid to transport stress syndrome of beef cattle (TSSBC), as their rumen digestive system needs time to recover from the dysbiosis caused by transportation stress (Deng et al., 2017;Li et al., 2019;Alfaro et al., 2020). To avoid the losses caused by TSS and better safeguard animal welfare, it is essential to provide a better condition and immune enhancing supplements for cattle during transportation. Beside extra-care it is also important to gain an in-depth understanding of the molecular mechanisms underlying TSS.
Comparing to long distance transportation stress, the molecular mechanism of short distance transportation stress is poorly understood. RNA sequencing (RNA-seq) is the method of choice for the high-throughput analysis of gene expression due to its low cost and ease-of-use, and offers a better guide for exploring the unknown information in transportation stress (Poplawski et al., 2016;Hrdlickova et al., 2017). Considering the potential application in transportation stress, blood is a suitable medium for samples collection and has low damage to cattle. This study aimed to examine the change in hematological indices and blood transcriptome after transportation, thereby providing the corresponding theoretical basis for the diagnose of TSSBC and a reference for the cost-control of beef cattle production and the safeguarding of animal welfare.

Animal Model of Short Distance Transport Stress
Ten healthy Qinchuan beef cattle (3∼4 years old) were used in the study. The animals were kept in loose housing conditions before transportation and fed on a total mixed ration (TMR). The beef cattle were deprived of food and water until after transportation. Blood samples were collected from each animal's jugular vein before transportation, and marked B1∼B7. Blood was collected about 17 h after last feeding. Thereafter, transportation has begun 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 a head which was similar to production process. The transportation distance was 70 km, and the vehicle was a single-layer Foton cart with a maximum speed of 60 km/h and an average speed of 30 km/h. The temperature ranged from 26 to 32 • C, and the humidity was 70%. There was no other animal transported with them. After transportation, blood samples were collected again from each animal's jugular vein immediately, and marked A1∼A7. The interval time between the two blood samples collection was 6-8 h.

RNA Isolation and Library Preparation
Blood for RNA-seq and quantitative real-time PCR (qPCR) was collected in sodium heparin anticoagulant containing tubes and stored at −80 • C until further analysis. To ensure the quality of RNA in samples, the following steps were performed as quickly as possible. ACK lysis buffer lysed erythrocytes, karyocytes were harvested following centrifugation at 4,000 rpm for 6 min at 4 • C, and then transfer into eppendorf tubes containing trizol, and finally stored in liquid nitrogen (Braga et al., 2017;Huang et al., 2017). Total RNA from the blood samples (n = 7) were isolated from leukocytes using the RNAiso Plus Kit (Takara, Tokyo, Japan). RNA concentration and quality were evaluated using Qubit2.0 fluorometer (Invitrogen, California, United States) and NanoDrop 1000 (Thermo Fisher Scientific, Massachusetts, United States). cDNA libraries were prepared using Hieff NGS TM MaxUp Dual-mode mRNA Library Prep Kit for Illumina R (YEASEN Biotech Co., Ltd, Shanghai, China). Libraries were quantified using Qubit2.0 Fluorometer DNA Assay Kit (YEASEN Biotech Co., Ltd, Shanghai, China) and submitted for sequencing (Illumina Xten, San Diego, United States).

Primer Design and qPCR
qPCR was performed to verify the RNA-seq expression pattern of blood before and after transportation. β-actin (ACTB) was used as the housekeeping gene ( Table 1). Reverse transcription of RNA to cDNA was performed before qPCR, carried out in the 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. All expression levels were normalized to that of ACTB and quantified using the 2 − Ct method (Livak and Schmittgen, 2001). ACTB was not different between animal groups before and after transportation. Twelve differentially expressed genes (DEGs) (fold-change > 2 and False discover rate (FDR) < 0.05) were randomly chosen for verification by qPCR. Hub gene (degree top 30) belonged to DEGs were also chosen for verification by qPCR.

The Effects of Transportation Stress on Hematological and Biochemical Indices
The analysis of hematological indices showed significant changes in 11 of the 18 indices evaluated. The WBC, MON#, Gran#, Gran%, MCHC, PLT and PCT% were higher after transportation compare to before transportation. The RBC, HGB, LYMPH% and HCT% were lower after transportation compare to before transportation (Figure 1).
The analysis of biochemical parameters showed significant differences in 5 of the 7 indices assessed. The TG, CHO, HDL, LDL, and MDA contents were lower after transportation compare before transportation, whereas GLU and SOD levels did not differ significantly between the two stages (Figure 2).

RNA-Seq Data and Verification by qPCR
Fourteen libraries representing seven animals in each group (before and after transportation) were prepared from total leukocyte mRNA, and group B was the reference group. The range of reads count for the 14 samples was 34888216∼43650486. For each sample, the number of total mapped reads was >90%, the Q20 base ratio was >97%, and the Q30 base ratio was >92% ( Table 2). Correlation analysis showed that 14 samples were FIGURE 2 | Comparison of serum biochemical parameters between before and after transport. TG, triglyceride; CHO, cholestenone; HDL, high-density lipoprotein; LDL, low-density lipoprotein; GLU, glucose; SOD, superoxide dismutase; MDA, malondialdehyde.  classified into two groups (group B and group A; Figure 3A). The result of principal component analysis (PCA) was same as correlation analysis, 14 samples also were divided into two groups (group A and group B; Figures 3B,C). There were 1,092 DEGs with FDR <0.05 and fold change >2 as the cut-off criterion, including 329 up-regulated and 763 down-regulated genes (Supplementary Material 1). Twelve DEGs (fold-change >2 and FDR <0.05) were randomly chosen for verification by qPCR. Nine genes showed the same trend as the RNA-seq data (S100A9, MS4A1, RGS1, CD19, CEBPD, ATP6, CD79B, CEBPB, and BIRC3), two genes showed the opposite trend (CD14 and BPI) and another gene was not significantly different between before and after transportation by qPCR (BANK1). The consistency of verification was 75% (Figure 4).

Enrichment Analysis of DEGs Between Before Transport and After Transport
To identify key differences between before and after transportation, GO and KEGG enrichment were performed to determine DEGs' function. The biological process (BP) of GO terms showed significant genes those are involved in cell activation, immune effector process, activation of the immune response, immune system process, regulation of immune system process, etc (Figure 5). The molecular function (MF) of GO terms shows the protease binding, regulatory region nucleic acid binding, RNA polymerase II transcription factor complex, RNA polymerase II regulatory region sequence-specific DNA binding, C-5 sterol desaturase activity, etc (Figure 5).

The most enriched GO terms in cellular component (CC)
were chromosome/telomeric region, nucleus, intracellular, nucleoplasm, centrosome, etc. (Figure 5). The KEGG analysis showed 200 pathways were significantly enriched. The pathways were mainly cell cycle-yeast, platelet activation, fatty acid metabolism, steroid biosynthesis, and intestinal immune for IgA. The KEGG analysis results showed that the DEGs were mostly associated with the immune response, lipid metabolism, hormone regulation, etc (Figure 5). The other pathways have been presented in Supplementary Materials 2, 3.

Co-expression Analysis, Key Module, and Hub Genes Identification
The median absolute deviation (MAD) top 5,000 genes were used to construct a co-expression network (Supplementary Material  4). 14 samples were clustered in two branches ( Figure 6A).
According to the results of scale independence and mean connectivity, we took power as 8 when correlation index Fold-change of DEG levels for both Log (TPM) and qPCR. S100A9, S100 calcium binding protein A9; MS4A1, membrane spanning 4-domains A1; RGS1, regulator of G protein signaling 1; CD19, CD19 molecule; CD14, CD14 molecule; CEBPD, CCAAT enhancer binding protein delta; ATP6, ATP synthase F0 subunit 6; CD79B, CD79b molecule; BANK1, B cell scaffold protein with ankyrin repeats 1; CEBPB, CCAAT enhancer binding protein beta; BIRC3, baculoviral IAP repeat containing 3; BPI, bactericidal permeability increasing protein; ACTB, actin beta. reached 0.85 (Figures 6B,C). All the genes were clustered into eleven modules ( Figure 6D). The modules genes see Supplementary Material 5. The results of network heatmap indicated a significant difference among modules ( Figure 6E). The most extensive module was the turquoise module, containing 2,546 genes. The smallest module was purple module, containing 68 genes ( Figure 7A). After KEGG enrichment in different modules (Supplementary Materials 6, 7), the turquoise module was the most essential module associated with short distance transportation stress, including immune reaction and metabolic pathway ( Figure 7C). Then PPI protein network of turquoise was constructed in STRING (Supplementary Material 8). The Cytohubba plugin based on Cytoscape was used to perform the PPI network analysis (Supplementary Material 9) and the top 30 genes based on degree algorithm were identified as potential hub genes (Figure 7D), which contains 9 DEGs ( Figure 7B). Considering fold change of qPCR was more comfortable to detect in practical production and the gene function of B cell receptor signal pathway, structural maintenance of chromosomes 3 (SMC3), jun proto-oncogene (JUN) and C-X-C motif chemokine ligand 10 (CXCL10) could be considered as molecular markers in short-distance transportation stress (Figure 8).

DISCUSSION
Beef cattle transportation is an inseparable part of the modern beef cattle industry, and since 2 decades, both transnational and cross-regional transportation were increased significantly. The neglect of beef cattle transportation stress by producers renders it difficult to accurately estimate the beef cattle industry (Deters and Hansen, 2020). It has been estimated that the death rate due to beef cattle transportation is in the range of 1-14%. Transportation stress is multifaceted but can be divided into long distance and short distance transportation stress (Marti et al., 2017). Long distance transportation stress can lead to a decline in the body's resistance, as well as a comprehensive, pathological, whole-body reaction syndrome (such as fever, depression, loss of appetite, thick coat, cough, asthma, sticky or purulent nasal fluid, diarrhea, bloody stools, arthritis, conjunctivitis, and even death) (Padalino et al., 2017;Miranda-de la Lama et al., 2018).
Comparing with long distance transportation stressors, short distance transportation stressors were less, mainly include fear, cold or heat, and lack of food and water (Allano et al., 2016). Additionally, less attention has been paid to short distance transportation stress as it does not result in obvious economic losses, and most of the symptoms disappear automatically (Simova et al., 2017). In rare cases, post-stress symptoms can occur, and drug intervention is needed (Van Engen et al., 2016;Capik et al., 2017;Wu et al., 2020). Mitigating beef cattle transportation stress mostly relies on experience and lacks a theoretical scientific basis. Consequently, the relevant underlying molecular mechanisms must be determined to alleviate transportation stress. In this study, we constructed a short distance transportation stress model to obtain the biological blood information and identify the molecular mechanism related to the TSSBC, thereby providing a reference for alleviating this condition. We collected the venous blood of beef cattle before and after transportation, and subsequently measured the blood hematological and biochemical indices, and analyzed the blood transcriptome.
Numerous indices were recorded in hematological parameters. The increases of MON#, Gran#, Gran%, MCHC, PLT, and PCT% were related to enhanced immune response following external stimulation. The increases of WBC was an essential phenotype in transportation stress (Dalin et al., 1993;Endo et al., 2017). Meanwhile, the decreases of RBC, HGB, and HCT% may simply have resulted from irregular drinking and eating before and after transportation. Finally, the high WBC and PLT were related to handling, loading, traveling and intense physical activity during loading (Pascual-Alonso et al., 2017). We found that the TG, CHO, HDL and LDL contents were significantly decreased for blood biochemistry parameters. This might be related to the fasting and lack of drinking water during transportation. Short term fasting could promote energy uptake by tissues and organs, but the time of fasting in this study could not enough for glycogenolysis and the following lipidolysis. Besides, the level of MDA showed that cattle have increased resistance to oxidative stress (Huang et al., 2000). SOD level was insignificant between group B and group A, but had an upward trend. It showed that most cattle improved ability to scavenge oxygen free radicals (Weismann et al., 2011).
WBC were the only nucleated cell in the blood, and were used to compare the gene expression between two groups by RNAseq. We first evaluated the quality of the sequencing data, and found that the number of total mapped reads was >95%, the Q20 base ratio was >97%, and the Q30 base ratio was >92% for each sample, indicating that the quality of RNA sequencing data was reliable and could be used for further investigation. Next, we identified DEGs based on the sequencing data by DESeq2. To verify the accuracy of the sequencing data, we randomly selected 12 DEGs for qPCR and found that the results had 75% accuracy, which implied that the data were suitable for further analysis. FIGURE 9 | Workflow and molecular mechanism in this study. In this study, we constructed a short distance transportation stress model of beef cattle. Blood samples were collected and performed to RNA-seq, GO, KEGG, and WGCNA. Comparing with hematological and biochemical indices, we concluded that B cell function, lipid metabolism and steroid function were the mean molecular mechanism in short distance transportation stress of beef cattle. CD19 and CD40, the markers of B cell, were down-regulated in transcriptome data. The results of WGCNA showed that B cell function could be the core of transportation stress. After hub gene identification and qPCR verification, SMC3, JUN, and CXCL10 could be suggested as potential biomarkers for the diagnosis and recover of this condition.
We performed GO and KEGG enrichment analysis of DEGs. Combining with the BP enrichment results top 20 terms showed that the level of humoral immunity with B cells as the core decreased after transportation (Honda et al., 2015). The decrease of humoral immunity explained that cattle were more vulnerable to disease after transportation (Quinteiro-Filho et al., 2010Calefi et al., 2014;Gomes et al., 2014). It was suggested that good management should be provided after transportation to help cattle recover from transportation stress. Like BP enrichment, KEGG enrichment results showed that B cell receptor function had a significant change after transportation. Besides, fatty acid metabolism and steroid biosynthesis were enriched between two groups, which was same as the results of the changing of TG, CHO, HDL, LDL in blood after transportation, and showed that hypothalamus-endocrine organ axis could be the core of transportation stress regulation. We hypothesized that B cell dysfunction might be the hematological molecular mechanism of short-distance transport stress, and the decrease of humoral immunity could be the potential factor behind the economic loss.
The main objective of this study was utilize a global approach to construct a gene co-expression network that could predict clusters of candidate genes involved in transportation stress. We hypothesized that tightly co-expressed gene modules, enriched in shared functional annotation, would provide the most effective predictions of candidate gene sets that might underlie a given biological process. Here, WGCNA was applied to investigate MAD Top 5,000 genes among the 14 samples in this study, and constructed 11 independent modules. We applied KEGG to evaluate the function of genes in different modules. Finally, turquoise module was enriched in metabolic pathway and immune reaction related to transportation stress. The pathway enriched in other modules related to transportation stress had no direct connection with transportation stress. It is reminding that housing disinfection and more attention to cattle's health were effective methods for reducing losses even short distance transportation. All the genes in turquoise module were used to construct PPI network, and degree top 30 genes were identified as hub genes by Cytohubba plugin based on Cytoscape. CD19 and CD40 were two specific B cell surface markers and have an important essential function of B cell maturation. The other hub genes were associated with the function of immune, differentiation, proliferation, survival, and apoptosis. Combing with the results of GO and KEGG enrichment, the function of hub genes was centered around B cell function (Hu et al., 2018;Gazon et al., 2018;Karin and Razon, 2018;Tokunaga et al., 2018;Xiao et al., 2019;Huang et al., 2020;Yang et al., 2020). Considering fold change of qPCR was easier to detect in practical production and the gene function of B cell receptor signal pathway, SMC3, JUN and CXCL10 could be suggested as molecular markers in transportation stress. SMC3, important member in cell cycle pathway, could be the regulator in the differentiation and activation of B cell function (Hu et al., 2018). JUN was belonged to many pathways related to B cell function, and had other functions containing differentiation, proliferation, survival, and apoptosis (Gazon et al., 2018 ; Figure 9). As a critical factor in immune activation through paracrine signaling, the decreasing of CXCL10 showed that immune cell differentiation was weakened, and the cattle susceptible to infection by germs (Karin and Razon, 2018;Tokunaga et al., 2018).
For the first time of this study, WGCNA was used to identify co-expression gene sets associated with short distance transportation stress. Turquoise module was found to be highly enriched for genes involved in the immune response and metabolism pathway. Collectively, the blood RNA-seq analysis and WGCNA indicated that the changing of B cell differentiation, proliferation, survival, and apoptosis were the potential threats in short distance transportation stress. Besides, SMC3, JUN, and CXCL10 were three potential molecular markers for short distance transportation diagnosis and recover.

DATA AVAILABILITY STATEMENT
The high-throughput sequencing data of the original RNAseq has been saved in the NCBI Sequence Reading Archive (https://www.ncbi.nlm.nih.gov/sra/), with the accession number PRJNA658391 (SRR12800407-SRR12800420).

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee of Northwest A&F University.

AUTHOR CONTRIBUTIONS
HZ, MW, and SW analyzed the data. HZ, XT, XY, and SL wrote the manuscript. HZ, XT, MW, QL, XY, SL, and JJ collected the samples. QL performed the qPCR. XY and SL reviewed and edited the manuscript. XS designed the experiment. All authors contributed to the interpretation of the results and writing of the manuscript.

ACKNOWLEDGMENTS
We thank Qinbao Cattle Industry Co., Ltd. for sample collection.