Differential gene expression in the Longissimus dorsi of Nguni and Bonsmara bulls finished on low and high energy diets

Objectives of this research were to examine differential gene expression profiles of Nguni and Bonsmara cattle fed diets differing in their energy density. The ultimate goal was to improve understanding of the mechanisms that underlie differences between these breeds and the potential interactions of the differences between breeds with the nutritive environment. The experiment was designed as a 2 × 2 factorial arrangement of breed and diet (12.5 MJ/kg DM vs. 10.9 MJ/kg DM). The initial feeding trial had 10 bull calves per treatment. However, financial constraints limited RNA sequencing to six animals per treatment and the RNA generated from one animal was of insufficient quality to be useful. Transcripts with false discovery rate P -values <0.05 and fold-changes >2.0 were considered significant. Bonsmara had a faster growth rate, heavier live and carcass weights, and better feed conversion compared to Nguni. However, lower levels of fat were observed in Nguni. Twenty different genes were differentially expressed, with three exhibiting interaction effects and all 20 having differences in transcript abundance between the breeds. A dietary effect was only observed for the one gene and that gene was also subject to an interaction effect with breed. Observed differences in gene expression between Bonsmara and Nguni by several genes affecting the structure or function of the mitochondria imply differences in energy metabolism between the breeds. Interaction effects on the abundance of some gene transcripts indicate the need to consider the diet when evaluating breed differences and conversely, consider breed when evaluating diets.


Introduction
The goal of a differential gene expression study is to determine which genes are expressed at different levels between the experimentally imposed conditions. These genes can offer biological insight into the processes affected by the condition(s) of interest. Identification of differentially expressed genes between specific conditions is a key in the understanding phenotypic variation (Costa-Silva et al., 2017;McDermaid et al., 2019). Examples of the successful use of this approach include providing biological insight into 1) parental genome effects on bovine foetal development between taurine and indicine breeds; and 2) differences in skeletal muscle differentiation between Wagyu and Piedmontese (Lehnert et al., 2007). Higgins et al. (2019) used differential gene expression to investigate the hepatic transcriptome of Charolais and Holstein-Friesian steers divergent in their residual feed intake when fed concentrate and forage-based diets. Thus, the use of differential gene expression methodology was implemented to study differences between Nguni and Bonsmara cattle when fed diets that differed in their energy density.
Nguni cattle are a tropically adapted and admixed Sanga-type breed that are indigenous to southern Africa (Scholtz et al., 2011a;Makina et al., 2016). These cattle are also widely used in communal and small holder farming systems on natural grazing where they contribute to the food security of resource-poor farmers in South Africa (Bester et al., 2003;Mapiye et al., 2009;Morokong, 2016). Females of this small-framed breed are often used in crossbreeding systems due to their low maintenance requirements, high fertility, and good mothering ability (Scholtz et al., 2011b). However, because Nguni cattle have a propensity to deposit excess fat at less-than-optimal carcass weights when fed high-energy diets, they are not preferred for finishing under feedlot systems (Scholtz et al., 2008;Strydom, 2008).
The majority of beef in South Africa is produced in feedlots using commercially-formulated, high energy diets, with preference given to medium-and large-framed and later maturing cattle (Scholtz et al., 2008;Strydom, 2008). Due to the Nguni's early maturity and lower meat yield, it is not preferred as a feedlot animal (Scholtz et al., 2008). At age-constant endpoints, du Plessis and Hoffman (2007) found Nguni steers were clearly fatter than Bonsmara-cross steers, even when finished on natural pastures. Nguni and Bonsmara cattle have been shown to produce similarly high quality meat (Strydom et al., 2000;Muchenje et al., 2008). One time-honoured strategy for feeding early maturing breeds to a heavier weight without depositing excess fat is to provide them a less energy-dense diet for a longer period of time (Owens et al., 1995).
The objective of this research was to examine differences in the gene expression profiles in the longissimus thoracis et lumborum of Nguni and Bonsmara cattle fed diets of differing energy density. The ultimate goal of this research was to obtain an improved understanding of the molecular mechanisms that underlie the differences between these breeds and potential interactions with their nutritive environment. At this time, this is the only study of differential gene expression in livestock that are indigenous to South Africa other than a study of Boar and village goats under intensive and extensive management (Ncube et al., 2022).

Materials & Methods
Ethical approval was received from both the Animal Ethical Committee of the Agricultural Research Council (ARC) and the AEC of the University of Pretoria (ec090-15).
At Irene, the calves were subjected to a 28-d adaptation period, followed by an 84-d feeding trial. The trial was designed as a 2 × 2 factorial arrangement of breed and feed level (high vs low energy diets). For each breed-diet treatment, there were 10 animals. The low energy and the high energy diets are shown in Table 1. The bulls were individually fed using Calan gates (American Calan, USA) and water was freely available. The calves were weighed weekly and images of the longissimus muscle and rump were collected at the beginning and end of the feeding period with real-time ultrasound (RTU) (MyLabOne, Esaote, USA). Traits measured on the RTU images were longissimus muscle area (EMA), subcutaneous fat depths over the rib (RF) and rump (P8), and intramuscular fat content of the longissimus (IMF). After being fed 102 days, the animals were slaughtered at the local abattoir and the carcasses were weighed. The abattoir was 3 km from where the animals were fed and they were transported the morning of slaughter.
Data from the feeding trial were analysed using the PROC MIXED of SAS v 9.4 (SAS Institute, Inc., Cary, North Carolina, USA) as a 2 × 2 × 2 factorial of farm of origin, breed, and diet. All effects were tested for significance against the residual variance. Interaction effects that included farm of origin were generally not significant (P >0.05) and were not considered further.
At slaughter, after each calf was stunned, hoisted and bled, a cut was made in its dorsal right side at the 13th rib. An 8 mm Miltex® biopsy punch (Integra LifeSciences, Princeton, New Jersey, USA) was used to sample the longissimus muscle. Forceps were used to place the muscle sample into a cryotube, which was then flash frozen in liquid nitrogen and the samples were stored in a -80 ºC freezer until RNA extraction. Due to funding limitations, samples from 24 animals (six per treatment) were selected for RNA extraction. Three replicate tissue samples (120 mg each) from each animal were processed. TRIzol ® Reagent (1 ml per sample) (Ambion, USA) was used to extract the RNA followed by homogenization. The sample was placed in a GenoGrinder 2010 (SPEX sample prep, USA) for 15-20 min at 1730 rpm. Chloroform was added to separate the DNA and the RNA. The homogenized and separated tissue was subsequently precipitated with isopropanol. In the final step, the precipitated pellet was washed with 1 ml 70% ethanol. Each sample was treated with the RNase-free DNase set (Qiagen, Hilden, Germany) to decrease the risk of genomic DNA contamination, and it was purified with the RNeasy mini kit according to the manufacturer's guidelines (Qiagen, Hilden, Germany). Extracted RNA from the samples was stored in a -80 ºC freezer. The extracted RNA was quantified using a Qubit Fluorescent meter (Invitrogen, USA) according to manufacturer's protocols and with a 1% gel using ethidium bromide and Tris-borate-ethylenediaminetetraacetic acid buffer. Following quality control, the RNA from one Bonsmara calf that was fed the high energy diet was discarded. Samples from the remaining 23 animals had RNA with concentrations between 221 μg/ml and 570 μg/ml and were of sufficient quality for sequencing.
Sample preparation was done using the TruSeq stranded mRNA protocol (Illumina, USA) with 3-5 µg RNA for each sample. Sequencing was performed using a HiSeq 2500 (Illumina, USA) for generation of pair-ended 2 × 125 bp reads. Of the 69 sequences, 13 were removed due to poor quality reads. Thus, there were 15, 3, and 5 animals with 3, 2, and 1 technical replicates, respectively. Technical replicate sequences were combined into a single file for each animal. Sequence reads were imported into CLC Genomics Workbench, version 22.0 (Giagen, Aarhus, Denmark) for bioinformatics analysis. Reads were trimmed using Q = 20 as a quality threshold and adapters were removed prior to genome mapping. Reads were mapped to the ARS-UCD1.2 bovine genome assembly (Rosen et al., 2020) using a length fraction of 0.8 and minimum similarity fraction of 0.8. Mapped read counts were normalized with the trimmed mean of M samples (TMM) method (Robinson and Oshlack, 2010). Pairwise comparisons between 1) breeds, 2) diets, and 3) diet-by-treatment interactions were completed with a generalized linear model assuming a negative binomial probability distribution of read counts (McCarthy et al., 2012), with treatment as the independent variable. Transcripts with Benjamini-Hochberg False Discovery Rate (FDR) P-values <0.05 (Benjamini and Hochberg, 1995) and foldchanges >2.0 were considered statistically significant.

Results and Discussion
Average weights of the Nguni and Bonsmara bulls at the start of the trial were 193 ± 4.29 kg and 206.49 ± 4.53 kg, respectively. Live weight at slaughter, ADG, FCR, P8, and EMA differed markedly between the breeds (Table 2). Diet had an effect on IMF and carcass weight (P <0.05). A breed by diet interaction (P <0.05) also affected IMF. Farm of origin had no effect on the various traits (P >0.05). Bonsmara had a faster growth rate, heavier live and carcass weights, and better feed conversion than Nguni cattle. However, lower levels of fat were observed in Nguni. Previous research has indicated that Nguni will be slower growing and weigh less at slaughter than Bonsmara, but the breeds will be similar in feed conversion and carcass fat content (Strydom et al., 2001). Thus, the samples of the breeds used in the current study may be considered more or less consistent with their previous characterizations. Samples from the animals produced a median of 99.5 million RNA transcript reads in pairs. On average, 93.6% of the reads mapped to the ARS-UCD1.2 bovine genome assembly (Rosen et al., 2020). Broken read pairs (those whose two reads didn't map to the same gene) were discarded, resulting in an average of 78.5% of reads subsequently being used. Of these counted reads, an average of 73.4% were mapped onto known genes (62.3% on exons and 11.3% on introns), while 26.6% mapped onto intergenic regions of the reference genome, suggesting the possibility of additional genes beyond those that are currently annotated.
Twenty genes were differentially expressed between Nguni and Bonsmara breeds, with 12 of these being more abundant in Nguni (Table 3). Differences in transcript abundance may be correlated with phenotypic differences in tissue growth and function. However, transcriptome differences should not be over-interpreted. Differences in transcript abundance do not always result in different protein abundance nor do they necessarily cause changes in phenotype. In this work, the transcripts were mapped to the genome and differences in gene expression were then quantified at the gene level. These gene transcripts corresponded to functions associated with ATP production and mitochondrial function, growth, appetite stimulation, and ribosome and cytoskeletal structure (Clark et al., 2016). Ubiquinol-Cytochrome C Reductase Core Protein 2 -7.61 4.9 x 10 -2 PECAM1 Platelet and Endothelial Cell Adhesion Molecule 1 1.96 4.9 x 10 -2 FDR: False discovery rate 1 Positive values indicate greater abundance in Bonsmara and negative values indicate greater abundance in Nguni 2 Identified by searching Homo sapiens and Mus musculus gene lists using Panther v17 (pantherdb.org) Three gene transcripts were differentially abundant in the contrast that measured the breed by diet interaction. These transcripts illustrate that there may be instances when the joint consideration of both breed and diet are needed for the proper interpretation of treatment effects. The ubiquinolcytochrome c reductase core protein 2 (Uqcrc2; fold-change = -19.1; FDR P = 6.0 x 10 -7 ) is located in the mitochondrion, where it is part of the ubiquinol-cytochrome c reductase complex. This complex constitutes a part of the mitochondrial respiratory chain responsible for electron transport (Guo et al., 2017). Previously, Kelly et al. (2011) failed to detect differences in the expression level of Uqcrc2 in the longissimus dorsi muscle between groups of Limousin × Friesian heifers fed high and low forage diets. Expression of transcripts of ENSBTAG00000031205 (fold-change = 31.5; FDR P = 2.5 x 10 -5 ) that encode 40S ribosomal protein S4 were likewise affected by the breed by diet interaction. Waerp et al. (2019) previously observed transcripts of this gene to be differentially abundant in adipose tissue from dairy heifers at 12 months of age versus when they were seven months pregnant when the heifers were fed a high energy diet but not when a lower energy diet was fed. The ENSBTAG00000049315 (foldchange = -178.5; FDR P = 6.6 x 10 -3 ) was the final transcript affected by the breed by diet interaction. These interaction effects provide further evidence for the necessity of considering the diet when evaluating breed differences and conversely, the breed when evaluating different diets (da Costa et al., 2013).
Only one uncharacterized gene transcript (ENSBTAG00000031205; fold-change = 19.8; FDR P = 4.0 x 10 -3 ) was observed more frequently in animals that were fed the low energy diet. A plausible explanation for this was that the differences between the diets were not sufficient to lead to changes in gene expression. As discussed above, the number of transcripts of this gene was also affected by the breed by diet interaction, indicating the effect of diet on its level of expression to be breed-specific.
Four genes associated with mitochondria or the electron transport chain were differentially expressed between Nguni and Bonsmara. Previous reports have documented relationships of feed efficiency with mitochondrial function in cattle (Kolath et al., 2006;Alam et al., 2012). The cytochrome c (Cycs) gene encodes a protein that is directly involved in the electron transport system, which produces ATP (Dayhoff, 1972). Similarly, the ubiquinol-cytochrome c reductase core protein 2 (Uqcrc2) encodes a protein that is part of complex III of the electron transport chain (Guo et al., 2017). Defects in electron transport also may be responsible for greater oxidized protein levels in less feed-efficient phenotypes (Bottje et al., 2006;Bottje and Carstens, 2009). The mitochondrial ribosomal RNA methyltransferase 2 (Mrm2) gene encodes a protein responsible for methylation of uridine in mitochondrial 21S rRNA (Pintard et al., 2002); deletion of this gene in yeast leads to a defect in cell respiration (Garone et al., 2017). The solute carrier family 25 member 6 (Slc25a6) gene encodes an ADP/ATP carrier that transports ADP/ATP across the mitochondrial membrane (Clémençon et al., 2013). Three of these four gene transcripts were more highly abundant in Nguni skeletal muscle, the exception being Slc25a6. Higher abundance of gene transcripts involved in mitochondria and oxidative phosphorylation suggest that ATP (and energy) production might be higher in muscle of the less feedefficient Nguni at time of slaughter. Higher ATP production might result from increased expression of oxidative phosphorylation genes, increased number of mitochondria, or both. Bottje and Kong (2013) found it logical that key genes associated with mitochondrial function, activities, or development would be differentially expressed between broilers that differed in feed efficiency.
Intriguingly, six genes associated with cytoskeleton structure or muscle fibres were differentially expressed between Nguni and Bonsmara skeletal muscle. Four of these gene transcripts (Dennd2a, Pwwp3a, Zwint, and Chrnd) were more highly abundant in Nguni. In chickens and turkeys, low feed efficiency has been associated with decreased abundance of transcripts associated with cytoskeleton structure and muscle fibres (Bottje and Kong, 2013;Pezeshkian et al., 2022). The Nguni calves in this study had higher FCR ( Table 2), suggesting that the greater expression of these genes in Nguni may contribute to their lower feed efficiency.
Two transcripts normally associated with the hypothalamus-anterior pituitary axis were more abundant in skeletal muscle tissue from Nguni calves. The growth hormone releasing hormone (Ghrh) gene encodes a peptide secreted from the hypothalamus that acts on the anterior pituitary gland to stimulate release of growth hormone (Grings et al., 1988;Alba and Salvatori, 2004). Melanin concentrating hormone receptor-1 (Mchr-1) gene encodes a receptor that binds melanin-concentrating hormone, which regulates appetite (Ingvartsen and Boisclair, 2001;Lord et al., 2021). Increased abundance of these transcripts in Nguni skeletal muscle may suggest increased feed intake at this age compared to Bonsmara calves. However, the action of both genes in skeletal muscle tissue has not been characterized. Further, increased transcript abundance of Ghrh and Mchr-1 may not lead to phenotypic differences in feed intake and growth rate.

Conclusions
Despite uncovering relatively few differentially expressed genes, fundamental differences in gene expression that underpin the metabolic processes of Bonsmara and Nguni were observed. In particular, several significant effects were noted for genes affecting the structure or function of the mitochondria. The presence of interaction effects on the abundance of some gene transcripts indicates the importance of considering the diet when evaluating breed differences and conversely, considering breed when evaluating diets.