Transcriptome profiling of longissimus lumborum in Holstein bulls and steers with different beef qualities

Previous research regarding Holstein cows has mainly focused on increasing milk yield. However, in order to maximize the economical profits of Holstein cattle farming, it is necessary to fully take advantage of Holstein bulls to produce high-grade beef. The present study aims to investigate different transcriptomic profiling of Holstein bulls and steers, via high-throughput RNA-sequencing (RNA-seq). The growth and beef quality traits of Holstein steers and bulls were characterized via assessment of weight, rib eye area, marbling score, shear force and intramuscular fat percentage of the longissimus lumborum (LL) muscle. The results indicated that castration improved the meat quality, yet reduced the meat yield. Subsequently, RNA-seq of the LL muscle from Holstein steers and bulls revealed a total of 56 differentially expressed genes (DEGs). We performed the functional enrichment analysis in Gene Ontology (GO) annotations of the DEGs using GOseq R package software and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using KOBAS tool. Through the integrated analysis of DEGs with reported QTLs and SNPs, seven promising candidate genes potentially affecting the beef quality of LL muscle following castration were discovered, including muscle structural protein coding genes (MYH1, MYH4, MYH10) and functional protein coding genes (GADL1, CYP2R1, EEPD1, SHISA3). Among them, MYH10, GADL1, CYP2R1, EEPD1 and SHISA3 were novel candidate genes associated with beef quality traits. Notably, EEPD1 was associated with both meat quality and reproduction traits, thus indicating its overlapping role in responding to hormone change, and subsequently inducing beef quality improvement. Our findings provide a complete dataset of gene expression profile of LL in Holstein bulls and steers, and will aid in understanding how castration influence meat yield and quality.


Introduction
Holstein cattle, a common cattle breed, has been mainly raised for dairy purpose. However, to meet the growing demand for beef, more and more Holstein cattle are also used to produce Food Co., Ltd. Feeding management of the two groups was kept consistent before selection and throughout the experimental period. Referring to Nutrition Requirement of Beef Cattle (2000) and Japanese Feeding Standard for Beef Cattle (2008), the basic diet formulation was adjusted along with the growth of body weight (S1 Table). The dietary nutrition levels of the two groups were identical. Following a 7 d adaption period, the experiment lasted for a total of 448 d.

Slaughtering, tissue and blood sampling
After 448 d of feeding, all of the cattle were driven to the abattoir (Fucheng Wufeng Food Co., Ltd., Langfang, China), located next to the feedlot. The cattle were feed fasted for 24 h and water fasted for 3 h before slaughter. Following electric shock, all cattle were slaughtered by exsanguination. The slaughtering process was performed according to the operating procedure of cattle slaughtering outlined in the China National Standard. Meanwhile, the hot carcass weight was recorded. Next, the carcasses were divided into two parts, weighed and chilled at 0~4˚C for 48 h. After chilling, 1 kg of LL muscle on the right carcass side between the 12th and 13th ribs was removed. Half of the muscle samples were stored at -20˚C for beef quality characterization, while the other half was cut into several pieces and placed quickly in liquid nitrogen, followed by storage at -80˚C until analysis.

Measurements of beef quality parameters
All of the measurements below were performed on the -20˚C stored muscle samples. The perimeter of the rib eye area was traced on a sheet of sulfuric paper, followed by calculation with a planimeter (Jilin University, China). The marbling grade was determined according to the Japanese Marbling scores (5-Excellent; 4-Good; 3-Average; 2-Below average; 1-Poor). Immediately after this, the meat sample was boiled and cooled, and the shear force was measured using a Warner-Bratzler shear force machine (Bodine Electric Co., Chicago, IL, USA) according to the manufacturer's instructions. According to the AOAC official methods, the meat samples were dried in a freeze dryer (LaboGene, Allerod, Denmark) to determine the contents of H 2 O, crude protein and IMF. The Kjeldahl method was then applied to analyze the crude protein content using a Kjeltec 8400 Analyzer Unit (Foss Analytical, Höganäs, Sweden). Finally, the total IMF percentage was determined by the Soxhlet extraction method.

RNA extraction and quality analysis
Three cattle were randomly selected from both the bull and steer groups for RNA-sequencing analysis. The total RNA was extracted from the tissues using TRIzol Reagent (Life Technologies, CA, USA) according to the manufacturer's instructions. A Qubit 1 RNA Assay Kit (Invitrogen, CA, USA) was used to determine the RNA concentration with a Qubit 1 2.0 Flurometer (Life Technologies, CA, USA). The RNA degradation and contamination were monitored by 1% agarose gel. The RNA purity was measured by a NanoPhotometer 1 spectrophotometer (IMPLEN, CA, USA). Finally, the RNA integrity was checked using a Bioanalyzer 2100 RNA 6000 Nano Kit (Agilent Technologies, CA, USA).

Library preparation and RNA sequencing
The RNA sequencing libraries were constructed according to the manufacturer's instructions of the NEBNext 1 Ultra™ RNA Library Prep Kit for Illumina 1 (NEB, MA, USA). Index codes were added, so as to attribute sequences to each sample. Next, mRNA was isolated using poly-T oligo-attached magnetic beads and fragmented by divalent cations in a NEB Next First Strand Synthesis Reaction Buffer (5X) under increased temperature. Subsequently, the first and second strand cDNA were synthesized using random hexamer primer. The 3' ends of DNA fragments were adenylated, followed by ligation with a NEBNext Adaptor. cDNA fragments of 150~200 bp were selected with an AMPure XP system (Beckman Coulter, CA, USA). PCR was performed with Universal PCR primers and an Index (X) Primer. The library quality was checked by an Agilent Bioanalyzer 2100 system. TruSeq PE Cluster Kit v3-cBot-HS (Illumina, CA, USA) was used for cluster generation. The final RNA-seq libraries were constructed using an Illumina Hiseq 2500 platform (Illumina, CA, USA), which generated 125 bp/150 bp paired-end reads.

Sequencing data analysis
In order to obtain clean data, raw data (fastq format) were processed through in-house perl scripts. In this step, clean data were obtained by removing reads containing adapters, low quality reads (the proportion of the read bases with Phred quality score�20 is over 50% of the reads), and reads with proportion of N greater than 10%. Subsequently, the Q20 (the percentage of read bases with Phred quality score >20), Q30 (the percentage of read bases with Phred quality score >30), and GC contents of the clean data were calculated. The downstream study analyzed the clean data with high quality.

Gene expression quantification and differential expression analysis
The gene expression level was estimated as fragments per kilobase of transcript per millions fragments mapped reads (FPKM), which was calculated based on the length of the gene and reads count mapped to this gene. HTSeq v0.6.1 was utilized for counting the read numbers mapped to each gene. Differential expression analysis of Holstein bulls and steers was performed using the DESeq R package (1.18.0). Benjamini and Hochberg's approach was applied to correct the resulting P-values to control the false discovery rate. Genes with a corrected Pvalue <0.05 were considered as differentially expressed genes (DEGs).

Gene Ontology (GO) and pathway enrichment analysis of DEGs
Gene Ontology (GO) of DEGs was implemented using the GOseq R package. KOBAS software was applied to calculate the statistical enrichment of DEGs in the KEGG pathways (http:// www.genome.jp/kegg/).

QTL-SNP screening analysis of DEGs
To associate the DEGs with beef quality traits, a QTL-SNP screening analysis was performed to determine the candidate genes [26,27]. Firstly, the DEGs with the physical position located within the region of reported QTLs related with beef quality traits (https://www. animalgenome.org/cgi-bin/QTLdb/BT/index) were selected. The genetic distance between the DEGs and QTL peak was calculated. Secondly, we created a pool of SNPs related with beef quality traits reported by previous GWAS. Physical position of the DEGs and the reported SNPs were compared, and the DEGs with distance to significant SNPs less than 5 Mb were selected. Finally, combining the QTL and SNP screening results, the overlapped DEGs those were both located inside of QTLs and near significant SNPs were selected as candidate genes associated with beef quality traits.

qRT-PCR
To confirm the sequencing results, nine DEGs were randomly selected for qRT-PCR. The total RNA (1 μg) was reverse transcribed into first strand cDNA by a Quantitect 1 reverse transcription kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Primers for qPCR were designed using the primer-BLAST tool on the NCBI website (S2 Table). The transcript levels of the tested genes were normalized to β-actin and GAPDH in qRT-PCR, and calculated using the 2 -44ct method. The qPCR reactions were performed in triplicates using iQ SYBR Green Supermix (BioRad, CA, USA). The amplified conditions were as follows: 95˚C for 3 min; 40 cycles of 95˚C for 15 s, 60˚C for 15 s and 72˚C for 20 s; 72˚C for 10 min.

Statistical analysis
The data of the growth and beef quality characteristics from the steers' group and the bulls' group obeyed normal distribution with the Kolmogorov-Smirnov test, followed by analysis with the Independent-Sample T test using SPSS 19.0 software. The results were shown as mean value ± standard error. A probability of P �0.05 was considered as significant, while P�0.01 was highly significant, and 0.05<P<0.1 was trendy. Pearson's correlations between transcript abundance from RNA-seq and qRT-PCR were analyzed using Minitab18 statistical software.

Growth and beef quality characteristics of Holstein bulls and steers
The growth and beef qualities characteristics of the bulls and steers are listed in Table 1. Castration negatively affected the carcass weight of the steers (P<0.05). In comparison with the bulls, the steers showed significantly decreased rib eye area (p<0.01), thus it can be concluded that castration resulted in slower growth and lower meat yield. However, the steers showed superior marbling (P<0.05) and higher IMF content (P<0.01) of LL muscle than the bulls. At the same time, the LL of the steers showed slightly lower shear force than the bulls, indicating that the LL of the steers possessed superior tenderness. Therefore, it was concluded that the castration significantly improved steer beef quality, and the LL samples could be used for RNA-seq to detect genes associated with beef quality traits.

RNA sequencing of bovine LL muscle
Next, a total of 305.57 million clean reads, with an average of 50.92 million for each sample (ranging from 47.80 to 54.43 million) were generated (S3 Table). The respective quality values of Q20 and Q30 were 95.29% and 88.84%. In this study, 88.32% of the total reads which  Table). Reads which could not be mapped (8.75%), or mapped to multiple positions (2.93%), were excluded from the following analysis (S4 Table). Through FPKM calculation, a total of 11,521 and 11,511 expressed genes were respectively detected in the bulls and steers. Among the expressed genes, a total of 10437 genes were common for the bulls and steers. In order to better characterize the difference of gene expression intensities, the FPKM values of gene expression were divided into five expression levels, from lowest to highest (S5 Table).

Differentially expressed genes (DEGs)
In order to characterize how the gene expression profiles of the LL muscle were impacted by castration, we compared the gene expression levels of the steers versus the bulls using DESeq method. A total of 56 significant DEGs were identified, among which 37 were upregulated and 19 were downregulated in the steers. Among these, 44 DEGs were known transcripts, of which 40 were transcripts of annotated genes and 4 were transcripts of pseudogenes. Twelve novel transcripts were detected, with 10 upregulated and 2 downregulated. The transcription profiles of the bulls and steers are shown in Fig 1. All of the annotated DEGs identified from comparison between the LL muscle of steers and bulls are listed in Table 2.

Validation of RNA-seq by qRT-PCR
In order to verify the RNA-seq results, we randomly selected 9 genes from the 41 annotated DEGs for qRT-PCR, namely IGFBP5, MYH1, PLCL1, SLC6A1, BDH1, MYH4, SLC30A3,  PETREG1 and ME2. The comparisons of transcript abundance detected by qRT-PCR and RNA-Seq are illustrated in Fig 2, showing the correlated gene expression levels using these two approaches. Consequently, it was shown that the RNA-seq data of this study are reproducible and convincible.

GO and pathway analysis of the DEGs
In order to further investigate the physiological characteristics that were affected by castration in the LL muscle, 56 DEGs were analyzed for their functions via Gene Ontology (GO). Most of the GO categories were associated with muscle development, cell division, various enzymatic activities and nucleotide metabolism (Fig 3). Specifically, the most enriched cellular components were myosin filament, cleavage furrow, cell surface furrow and myosin complex, most of which relate to muscle growth and development. The most highly enriched molecular functions include deaminase activities, oxidoreductase activities, vitamin D 3 25-hydroxylase activity and hydrolase activities, all of which either directly or indirectly participated in fatty acid metabolism. An additional GO term clearly associated with fatty acid synthesis in the top 30 enriched categories is a response to ketones. Meanwhile, the top biological processes consist of IMP salvage, purine nucleotide salvage, purine-containing compound salvage and nucleotide salvage, which were all related with the purine nucleotide cycle, which is the ultimate function in muscle energy production.

PLOS ONE
To discover the metabolic pathways associated with beef quality, we performed metabolic pathway analysis on the DEGs. The details of the top 12 enriched pathways in the LL muscle by comparing bulls and steers are listed in Table 3. The pathways which were related with fatty acids metabolism were synthesis and degradation of ketone bodies, butanoate metabolism, purine metabolism, galactose metabolism, linoleic acid metabolism, and pyruvate metabolism. Other pathways were involved in muscle cell growth, such as tight junction, oxytocin signaling pathway and proteoglycans in cancer.

Candidate genes associated with beef quality traits by QTL-SNP screening
In order to further screen the DEGs for the candidate genes related with beef quality traits, we analyzed the DEGs in the animal QTL database (https://www.animalgenome.org/cgi-bin/ QTLdb/index). By comparing the gene position of the DEGs on the chromosome with the QTLs region, 9 out of 40 annotated DEGs were discovered (S6 Table). In addition, we created a genome-wide association studies (GWAS) database by pooling the published SNPs related with meat quality traits for further analysis. The differentially expressed genes are considered as related to the specific SNP associated traits only if the position of the gene on the chromosome is at a distance of less than 5 Mb from that of the SNP (S7 Table). A total of 39 out of 40 annotated DEGs were found as potential genes related to beef quality. Taken together, seven overlapping DEGs, namely GADL1, CYP2R1, MYH1, EEPD1, MYH10, SHISA3 and MYH4, are considered as promising candidate genes responsible for the differences in LL muscle qualities between Holstein bulls and steers (Table 4).

Fig 2. Correlations of mRNA expression levels of 9 random DEGs in LL muscle of steers versus bulls by qRT-PCR and
RNA-seq. The x-axis indicates the log 2 (ratio of mRNA levels) using RNA-seq, and the y-axis displays the log 2 (ratio of mRNA level) measured by qRT-PCR. The blue dots represent the tested genes. The red line indicates the scatterplot of qPCR vs RNA seq. https://doi.org/10.1371/journal.pone.0235218.g002

Fig 3. Most highly enriched GO terms of DEGs between Holstein steers and bulls.
The x-axis displays the number of DEGs, and the y-axis represents the GO terms. The bar colors correspond to different GO categories, with yellow for molecular function, red for biological process and blue for cellular component. 1 Hydrolase activity, acting on carbonnitrogen (but not peptide) bonds, in cyclic amidines. 2 Oxidoreductase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor. https://doi.org/10.1371/journal.pone.0235218.g003

Discussion
In this study, the LL muscle of Holstein steers and bulls were biochemically characterized, and the transcriptomes of the LL muscle were profiled. The results showed that castration improved beef quality yet reduced the meat yield of Holstein bulls. We obtained a total of 56 DEGs, with 37 being upregulated and 19 being downregulated. Among them, seven promising candidate genes which were associated with beef quality upon castration were identified upon our integrated analysis, which included GO analysis, KEGG pathways enrichment analysis and screening with reported QTLs and SNPs. The seven genes included both muscle structural proteins coding genes (MYH1, MYH4, MYH10) and functional proteins coding genes (GADL1, CYP2R1, EEPD1, SHISA3). With the exceptions of MYH1 and MYH4, the remaining five genes are novel reported genes related to beef quality traits [34,35]. Our findings generate a complete data set of genes in the LL muscle that are affected by castration. MYH1, MYH4 and MYH10 are different isoforms of the myosin heavy chain gene family (MYH), and are responsible for the diversities of skeletal muscle fiber types. Myosin fibers are mainly composed of four types, namely class I, class IIA, class IIX and class IIB, which correlate to the muscle morphologies and biochemical characteristics, ultimately impacting meat quality traits [36]. MYH1 encodes the motor domain of class 2X isoform, while MYH4 and MYH10 both encode the motor domain of class 2B isoform. A previous study regarding the meat qualities of pork revealed that class 2B and 2X are the most abundant fibers representing the muscle mass, and affect meat quality based on a positive correlation with pH decrease ratio and marbling score [37]. Interestingly, our QTL-SNP screening for beef qualify related genes revealed that MYH1, MYH4 and MYH10 were all located within QTL regions associated with IMF, and close to SNPs associated with fatty acids composition and carcass conformation, thus indicating that these myosin heavy chain genes mainly affected beef quality through meat Table 4. Candidate DEGs associated with meat quality traits. composition (Table 4). Proteomic study revealed that MYH1 was associated with superior beef tenderness [38]. A microarray analysis at the transcriptomic level showed that MYH1 was related to the shear force of beef muscle [34]. Similarly, our RNA-seq analysis showed that MYH1 expression was significantly increased in the LL of steers with better beef quality, coordinating with its identified role in improving beef quality. Whole genome resequencing demonstrated that MYH4, as a positive marker of meat quality improvement, was related with the skeletal muscle development and muscle fiber type for red Angus cattle [35]. Our results revealed that the MYH4 mRNA level significantly increased in the steers in comparison to the bulls, indicating MYH4 may play a similar role in positively affecting beef quality traits. Although beef quality of the steers was significantly higher than that of the bulls, but the MYH10 level was relatively lower in the steers' muscle. Hence, MYH10 showed negative correlations with beef quality traits. This is also the first report that myosin heave chain gene MYH10 is associated with beef quality traits. Coincidently, a transcriptome analysis of porcine longissimus dorsi muscle suggested that, compared to purebred Iberian pigs, the MYH10 level was higher in crossbred Iberian pigs which contained lower IMF [39]. Therefore, both studies indicated that MYH10 is negatively correlated with meat quality traits.

Gene Information Reported QTLs Reported SNPs
Glutamate decarboxylase, such as protein 1 (GADL1), encodes a protein that is specifically expressed in cattle muscle and functions in β-alanine and taurine production, and is ultimately involved in carnosine and pantothenate biosynthesis [40]. GADL1 transcript level was upregulated in the steers, thus indicating its role in improving beef quality. Our beef quality related candidate genes screening revealed that GADL1 is located within a QTL region related to lean meat yield, at a position close to significant SNPs with carcass conformation and LL muscle area (Table 4). It has been previously observed that carnosine and taurine can promote the stability of IMF lipids and beef color [41]. In addition, nutrition studies regarding pigs and geese revealed that supplementing with pantothenic acid improves carcass characteristics, lipid metabolism and meat quality, thus indicating that pantothenic acid may play a similar nutrient role in cattle [42,43]. Our data suggested that the promising functions of carnosine, taurine and pantothenic acid in improving meat quality may be related to the gene GADL1. Overall, we proposed that GADL1 is likely associated with lipid metabolism and carcass, but this requires further investigation.
Cytochrome P450 family 2 subfamily R member 1 (CYP2R1) encodes an enzyme of the cytochrome P450 superfamily. CYP2R1 is highly conserved from zebrafish, chicken to mammals. Its human homolog catalyzes steroid biosynthesis, thereby implicating a similar role in cattle [44]. Since steroid is a precursor of androgen, the decrease of CYP2R1 mRNA level in steers may be straightly explained by lower androgen level after castration. At the same time, it was previously reported that vitamin D 3 supplementation improved beef tenderness in steers [45]. Coincidingly, the results of our QTL-SNP screening suggested that the CYP2R1 gene loci are close to an SNP associated with IMF percentage which is responsible for tenderness ( Table 4). The KEGG analysis suggested that CYP2R1 catalyzed vitamin D 3 , so as to produce calcidiol. The decrease of the CYP2R1 expression level in the LL of steers may cause the accumulation of vitamin D 3 , which further improved the beef tenderness. Therefore, we propose that CYP2R1 may be an important candidate gene that is related to beef tenderness, mainly through its involvement in vitamin D 3 accumulation. Shisa family member 3 (SHISA3) encodes a single-transmembrane protein of the Shisa family. Our QTL-SNP screening showed that SHISA3 is associated with beef quality traits such as shear force and IMF percentage, indicating the function of SHISA3 might be related with adipose deposition. Previous study revealed that SHISA3 plays a role as a tumor suppressor in lung cancer, by inhibiting the Wnt/β-catenin signaling pathway in lung adenocarcinoma cells [46]. Coincidently, Wnt/β-catenin signaling pathway was reported as potent inhibitor of adipogenesis [47]. A study on Korean cattle indicated that Wnt/β-catenin signaling pathway was downregulated following castration, which negatively corelates with increased IMF in longissimus dorsi muscle [48]. If SHISA3 regulates Wnt/β-catenin signaling pathway in muscle cells in the same way as lung adenocarcinoma cells, expression of SHISA3 should have shown positive correlations with beef quality traits. However, our study revealed that the transcription level of SHISA3 was downregulated in the steers compared to the bulls, showing a negative correlation with higher IMF content in steers. Since the function of SHISA3 in cattle has not been studied yet, it is worth to investigate the mechanism how SHISA3 impacts adipogenesis or adipose deposition in future studies. Regulation of SHISA3 on Wnt/β-catenin signaling pathway in LL muscle will be of particular interest in order to better explain our findings.
The endonuclease/exonuclease/phosphatase family domain containing 1 (EEPD1) encoded gene that mainly plays a role in repairing stressed replication forks during DNA damage in humans [49]. EEPD1 in cattle has not been studied extensively, but a transcriptome analysis of chicken adipose tissue revealed that a decrease of EEPD1 correlated with both fasting and insulin-neutralization, thus indicating its negative regulation in fatty acids synthesis [50]. Our QTL-SNP screening of DEG-related beef quality suggested that EEPD1 was located within a QTL region associated with carcass weight and fat thickness at the 12th rib, and near SNPs impact carcass conformation and marbling score (Table 4). Therefore, EEPD1 may affect lipid metabolism or fat deposition. Opposite to the results of the chicken study, in our study EEPD1 was upregulated in the LL of steers, which contains significantly higher IMF, indicating that it may play a positive role in fat deposition in cattle muscle. In addition, a study by Nelson et al. on human and murine macrophages proved that EEPD1 is regulated by LXR, which is a responsive receptor of high cellular sterol load, while EEPD1 contributes to the efflux of excess cellular cholesterol [51]. These findings indicate that EEPD1 may play a potential role as an androgen sensor in muscle. Subsequently, in order to detect whether EEPD1 could be a potential candidate gene associated with reproduction traits, we compared the physical location of the EEPD1 gene with the reported QTL regions and SNPs associated with cattle reproduction traits. It was clear that the physical location of the EEPD1 dropped into the QTL regions related with scrotal circumference and gestation length, and was very close to the SNPs associated with poor sperm motility and cow conception rate (S8 Table). Although the indicated QTLs and SNPs were discovered in reproductive organs, traits such as scrotal circumference, sperm motility and gestation length are all corelated with sex hormone levels directly or indirectly. From this prospective, and linking the research by the Nelson Group, we cannot exclude the possibility that the upregulation of EEPD1 in LL muscle of steers is a response to higher cellular steroid level induced by castration. Therefore, EEPD1 became a strong candidate gene with overlapping functions of sensing castration effects, subsequently improving beef quality.
Overall, our study reported 56 DEGs in the LL of Holstein bulls and steers. Our QTL-SNP screening analysis revealed seven candidate genes potentially associated with beef quality improvement caused by castration, including muscle structural proteins coding genes (MYH1, MYH4, MYH10) and functional proteins coding genes (GADL1, CYP2R1, EEPD1, SHISA3). Among them, GADL1, CYP2R1, EEPD1, SHISA3 and MYH10 were novel candidate genes associated with beef quality traits. EEPD1 is associated with both beef quality and reproduction traits, providing important insight as to how castration affected the beef quality at the gene level. Our study confirmed that the reported positive markers of meat quality traits MYH4 and MYH1 are associated with the beef quality improvement following castration. We also identified another myosin heavy chain family gene MYH10, which is negatively associated with beef quality traits. Hence, these data demonstrated that the myosin heavy chain family plays an important role in regulating beef quality. Due to the fact that our QTL-SNP screening was biased on screening genes associated with meat quality traits, it is possible that the DEGs related only to the sex hormone responses, but not meat quality traits, were omitted. The overall analysis indicated castration affected gene expression in LL muscle, resulting in beef quality improvement in Holstein steers.
Supporting information S1