Identification of Metabolic Pathways and Hub Genes Associated with Ultrasound Subcutaneous Fat and Muscle Depth of the Longissimus Muscle in Cull Beef Cows Using Gene Co-Expression Analysis

1. Introduction

Beef is an excellent source of protein, and its demand is increasing globally [1]. Beef production efficiency is directly linked to profitability and sustainability of the beef cattle industry [2]. Beef production is a complex system and influenced by different factors such as age, breed, diet, genetic factors (genes), weather, genotype-by-environment interaction and other factors [3,4]. Most of the genetic or molecular factors involved in beef production are not completely understood. Few genes related to beef production in cows have been reported. For example, ACACA (Acetyl-CoA carboxylase alpha) gene is reported as a major player in intramuscular deposition in ruminants. In beef cows, a positive correlation between ACACA expression level and total lipids and trans-fatty acids were reported [5]. Pleomorphic adenoma gene 1 (PLAG1) and Pleomorphic adenoma gene 2 (PLAG2) were reported to be associated with growth, stature, skeletal muscle growth, fat thickness, and fatty acid composition in cattle [6]. A haplotype analysis found evidence that the mutation in PLAG1 mutation affects body size and weight in modern cattle [7].

Next-generation sequencing allows detailed analysis of complex traits at the level of their genomes (genomics), transcriptomes (transcriptomics), or epigenome (epigenomics) [8]. Among these, transcriptomics is widely utilized to identify regulators of complex traits at gene/transcript level. Ribonucleic acid (RNA) sequencing (RNA-Seq) is a popular approach for transcriptome profiling that uses next-generation sequencing analysis. RNA-Seq gives detailed quantification of transcripts, and their isoforms compared to other techniques [9]. RNA-Seq has been utilized for studying gene expression profiles of different animals, e.g., cow, horse, sheep, buffalo and others [10,11,12,13]. Most of these studies have identified differentially expressed genes and their associated pathways. However, additional identification of genes and their relationship to important traits can shed more light on how genes regulate complex traits. Weighted Gene Co-Expression Network Analysis (WGCNA) is a widely used technique for co-expression network analysis to screen gene expression data and relationships with phenotypic traits [14]. WGCNA has been widely used to study complex traits and their correlation with gene expression data [15]. In cows, WGCNA has been used to investigate feed efficiency. One study identified ATP7B (ATPase Copper Transporting Beta) as a hub gene (central regulator) significantly associated with residual feed intake [16] In Nellore cattle, relationships of between gene expression and carcass traits were revealed using WGCNA [17]. Authors reported that energy and lipid metabolism pathways, highlighting fatty acid metabolism, were the central pathways associated with rib-eye area.

In an effort to further understand the complex relationships of production traits and gene expression in beef cows, we employed RNA-Seq based transcriptome data and WGCNA. In our study, body weight, body condition score, back fat and muscle depth were traits (variables) associated with beef production which were monitored. Total RNA was extracted from the longissimus dorsi with RNA-Seq analyses conducted. Weighted gene co-expression network analysis was applied to associate changes in beef production related traits with gene expression. Relationships between measured traits, metabolic pathways, and hub genes were then identified.

2. Materials and Methods

2.1. Experimental Site

All the procedures were conducted according to a protocol approved by the Institutional Animal Care and Use Committee at the University of Kentucky (protocol #2022-4039). The facility was inspected and approved by the Institutional Animal Care and Use Committee. Animals were individually penned (2.4 m x 14.6 m pens) in a partially covered concrete floored pen. The covered area was bedded with sawdust. Each pen contained an automatic waterer shared between two pens. Feed was delivered in concrete bunks each morning.

2.2. Beef Production Related Trait Analysis

Twenty-two mature cows were sourced from local livestock markets between May 23 and June 6, 2022. Purchased cows were of unknown age and genetic make-up. The majority of cows had black haircoats likely being of Angus influence while some were noted to have Charolais and Hereford influence determined by their coat color and pattern. Arrival body condition scores by three trained individuals and body weights were collected on June 10, 2022. After a few days, three of the purchased cows were removed due to lameness or other pre-existing conditions which were not improving. A replacement open cow, predominately Angus breeding, from the university herd was enrolled in the study to allow for 20 cows to be sampled at the start of the study. Cows were provided hay ad libitum for the first few days after arrival. Cows were then transitioned to a diet consisting primarily of corn silage and a protein supplement. The diet was limit-fed based on arrival weights to 1.5x maintenance energy requirement levels.

2.3. Longissimus Dorsi Muscle Biopsy

Muscle biopsies of the Longissimus dorsi were performed between 12th and 13th ribs on July 29, 2022. Ultrasonography was used to assess subcutaneous fat thickness over the Longissimus dorsi and muscle depth. The ultrasound was used to target the mid-point of the muscle to provide a consistent sample region between animals. Muscle samples were obtained from the right side of the animal. The target region was clipped to remove hair, scrubbed with a betadine solution and then alcohol. A local anesthesia injection of 5-ml of a lidocaine solution (2%) was given before biopsy incision. The incision through the skin was made using a 2.54 cm long 14-gauge needle. Muscle samples were collected using semi-automatic biopsy needles (Argon medical device, catalogue number: 701114090) to collect approximately 200 mg of muscle tissue. Excessive fat was removed from the biopsies. Samples collected for RNA extraction were suspended in RNAlater (RNAlater™ stabilization solution, Invitrogen, AM7020) to prevent RNA degradation during sampling. All the samples were flash frozen in liquid nitrogen and stored at -80^oC until RNA extraction.

2.4. RNA Extraction, Library Preparation and Sequencing

Muscle samples were subjected to total RNA extraction using the Zymo Research quick RNA micro prep kit (R01050). The RNA’s purity was assessed through gel electrophoresis and the Agilent 2100 Bioanalyzer. Isolation of pure messenger RNA (mRNA) was achieved using oligo (dT) beads. Following this, cDNA libraries were generated utilizing the Kapa Stranded mRNA Library Preparation Kit (Roche, KR0960) in accordance with the manufacturer’s guidelines. The cDNA library underwent size selection (150 bp) and was subsequently enriched through PCR. The quality assessment of the libraries was performed using a bioanalyzer (Agilent 2100 Bioanalyzer system). The sequencing of cDNA libraries was carried out on the Illumina NovaSeq S4 platform, producing paired-end (2 × 150 bp) data in FASTQ format. The sequencing library preparation and sequencing procedures were conducted by Quick Biology Inc.

2.5. RNA Sequencing Data Analysis

Raw sequencing reads (FASTQ files) were subjected to quality control by FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). In brief, FastQC detected contaminants including overrepresented sequences, sequencing adapters, poor quality reads and PCR duplicates. All these contaminants were filtered using Cutadapt 1.4 (Command line : cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq) [18,19]. The cow reference transcriptome files (Bos taurus ARS_UCD 1.2) for mapping and quantification were obtained from Ensembl genome browser (https://ftp.ensembl.org/pub/release-108/fasta/bos_taurus/cdna/). The initial step involved indexing the cow reference transcript file using the Salmon index option. (command line: salmon_index -t cow_transcripts.fasta -d decoys.txt -p 12 -i cow_salmon_index). The transcriptome-wide quantification in the form of TPM (transcript per million) on clean reads were performed using salmon v 1.6.0 (command line: salmon quant -i salmon_index -l A -1 Clean_Read_1.fq -2 Clean_Read_2.fq --validateMappings -o salmon_quant -p 12 --seqBias) [19]. The workflow of RNA sequencing and data analyses are shown in Figure 1.

2.6. Weighted Gene Co-Expression Network Analysis (WGCNA)

The R package, Weighted Gene Co-expression Network Analysis (WGCNA) was used to construct co-expression networks and find co-expressed genes [20]. Gene expression data (TPM) of 18 cows with four different traits (WT, BF, MD and BCS) were selected for WGCNA analysis. Two cows were removed due to being outliers. At first, the low count genes and outliers were filtered. To choose modules associated with traits (WT, BF, MD and BCS) of our interest, Pearson correlation analysis between each module’s “eigengene” and the traits were determined. Modules are groups of genes with similar expression profiles and tend to be functionally related and co-regulated. Modules with a module-trait correlation >0.4 or <-0.4 for at least one trait (P ≤ 0.10) were considered significant. Gene significance (GS) was calculated for each gene as the correlation between gene expression counts and traits (WT, BF, MD and BCS). Hub genes were identified by choosing genes with high gene significance and module membership in the modules of interest.

2.7. Functional Enrichment Analysis

Associated modules were selected (R > 0.4 or <-0.4, p≤ 0.10) for functional enrichment analysis. KEGG pathways and gene ontology terms enriched in these modules were identified using gprofiler interface (https://biit.cs.ut.ee/gprofiler/) with threshold of a Benjamini-Hochberg FDR (false discovery rate) of 0.05 [21]. We chose the “g:GOSt – Functional profiling of gene lists” option for conducting the metabolic pathway and gene ontology analysis within the selected Bos taurus reference database. The graphical representation of most significant pathways was done using R package ggplot2 (https://ggplot2.tidyverse.org/index.html). [22]

3. Results

3.1. Trait Relationship Analysis

We recorded four traits in these cows (Table 1). Covariation between these four traits showed that all four traits are positively correlated with each other as expected. Among these traits BCS and WT had correlation of 0.581, whereas MD and BCS also were positively correlated being 0.463. (Figure 2)

3.2. RNA Sequencing

RNA sequencing yielded an average of 33 million raw reads per sample. On average, 98.80% of the raw reads successfully passed the quality check. The average percentage of mappable clean reads per sample were 70.13%, ranging from 64.3 to 74.7% (Supplementary File S1). The transcript per million (TPM) for all 37,926 cow genes were estimated, the average TPM per sample was 26 (Supplementary File S2). Out of 37,926 genes, 7,655 genes were removed from further analysis due to excessive missing samples or zero variance.

3.3. Weighted Gene Co-Expression Network Analysis (WGCNA)

Sample dendrogram produced by hierarchical clustering is shown in Figure 3. Sample dendrogram based on gene expression and heatmap using trait data clearly divided the 18 cows into two groups: group 1- with mostly light weight cows (8 cows) and group 2- mostly cows having heavier live weights (10 cows). The WGCNA analysis constructed 21 modules (i.e., cluster of co-expressed genes) identified by a different color. The WGCNA also yielded correlations between the genes of each module and measured traits (WT, BF, MD and BCS). Figure 4 illustrates the association between modules and traits, which represents Pearson’s correlation coefficients measured between each single module and trait. Given the vast amount of data, we focused on two modules i.e., dark-green and magenta which showed significant positive and negative correlation with all four traits. Data for all genes, their respective modules (21 eigengene modules) and correlation values are given in Supplemental File S3.

3.4. Functional Enrichment Analysis

KEGG pathway analysis of the dark-green module showed ‘olfactory transduction’ as the most significant pathway. Gene ontology molecular function and biology process analysis of dark-green showed RNA binding and spliceosome snRNA activity as significantly enriched GO terms (Figure 5 and Supplementary file S4). Additionally, KEGG pathway analysis of the magenta module showed hematopoietic cell lineage as the most significant pathway. Gene ontology molecular function and biology process analysis of magenta showed protein binding and immune system process as significantly enriched GO terms (Figure 6 and Supplementary File S4).

3.5. Hub Gene Identification

For the dark-green module, CCDC88C (Coiled-coil domain containing 88C) and CACNA2D3 (Calcium voltage-gated channel 3) were identified as hub genes. These genes have gene significance >0.5 with back fat and module membership >0.9 with dark-green module eigengene. In the magenta module, SEPTIN9 and NONO (Non-POU domain containing octamer binding) genes were defined as hub genes based on their gene significance >0.8 with muscle depth and module membership >0.9 with magenta module eigengene. (Figure 7)

4. Discussion

The objective of the study was to identify candidate gene sets and metabolic pathways that could provide a genetic insight into tissue accretion of cull beef cows. Using gene expression data from the longissimus dorsi muscle tissue, gene co-expression modules were built and correlated to traits associated with beef production. WGCNA identified metabolic pathways and candidate genes that are expected to have a biological role in beef production. Out of 21 gene modules, the modules having the greatest significant correlations to the traits measured were further analyzed and discussed. Dark-green module (356 genes) had a significant positive correlation with back fat and muscle depth. The functional analysis highlighted olfactory transduction as a significant pathway and RNA activity (processing, splicing, and binding) as a significant GO term. In the olfactory transduction pathway, odorant or food is detected with signals sent to the brain or other body parts [23]. These signals are in general influenced by calcium, cyclic AMP and G protein molecules. In animals, olfactory transduction is directed linked with appetite, feed preference and weight gain [24,25]. In the pig, olfactory transduction pathway exhibited direct association with residual feed intake [26]. In beef cattle (SimAngus), quantitative trait loci detected for residual feed intake was flanked by eight olfactory receptor genes [27]. Thus, biologically animals with an upregulation of feed intake would be expected to have increased tissue accretion or at least for the traits measured in this study.

Another noteworthy pathway found in the dark green module is the hippo signaling pathway. Previous reports have shown that the hippo signaling pathway regulates the expression of Yes-associated protein (Yap) and transcriptional co-activator with PDZ-binding motif (Taz). Both Yap and Taz genes are critical regulators of skeletal muscle mass. Additionally, these genes are involved in controlling tissue growth in various cell types [28]. The significant gene ontology terms of dark-green modules were RNA binding, snRNA activity, RNA processing and RNA splicing. All these terms come under RNA processing and modifications. According to previous literature, muscle development relies on RNA processing. It is mainly driven by coordination among RNA binding proteins. Additionally, misregulation of RNA processing causes muscle diseases [29,30]. Alternative splicing is a mechanism that enables single genes to produce multiple mature mRNAs. Among various tissues, skeletal muscle cells, along with the brain and heart, display the most prominent levels of tissue-specific and evolutionarily conserved alternative splicing events. These molecular alterations are accompanied by a multitude of transcriptional and posttranscriptional changes, many of which are governed by alternative splicing mechanisms [31,32,33,34]. In mice, RNA splicing of the Fragile-X mental retardation autosomal homologue-1 (FXR1) gene directly contributes to muscle development. The FXR1 gene is hypothesized to play a regulatory role in mRNA translation, localization, and stability within muscle cells [35]. Thus, cull beef cows which are emaciated due to restricted nutrient intakes and then provided an increased plane of nutrition may elicit responses for protein accretion.

The magenta and turquoise modules had negative correlation with back fat, muscle depth and body condition score. The functional analysis of magenta highlighted hematopoietic cell lineage pathway as the most significant pathway. In animals, tissue self-renew, tissue growth and muscle tissue regeneration are mainly controlled by Hematopoietic cell lineage [36]. During muscle growth, Hematopoietic cell lineage pathway is activated. However, in the current study the opposite was observed showing less expression for increasing weight, BCS, MD and BF. Detailed characterization of this pathway is necessary to study the exact role of hematopoietic cell lineage in beef production. It is possible this response observed in this group of cull beef cows is an artifact of previous management in which nutrient intake was less than the animal’s requirements. In regard to the turquoise module, it was noted that the oxidative phosphorylation pathway was the most highly significant pathway. Oxidative phosphorylation provides energy in the form of ATP to support muscle growth and maintenance [37]. In buffalo, the WGCNA analysis identified pathways related to energy metabolism, such as oxidative phosphorylation, that play a role in fat deposition [38]. Thus, if animals were in a negative energy state previously, it would be plausible the expression of genes in this pathway may suppressed.

Hub genes are candidate genes which are important in the regulation of the expression of several other genes in a module and are potential biomarkers for the selection of a trait [39]. Such an approach has already been employed in ruminants to identify genes governing traits [40,41,42,43,44]. The WGCNA analysis in buffalo identified six hub genes [FH (Fumarate Hydratase), MECR (Mitochondrial Trans-2-Enoyl-CoA Reductase), GPI (Glucose-6-Phosphate Isomerase), PANK3 (Pantothenate Kinase 3), ATP6V1A (ATPase H+ Transporting V1 Subunit A), PHYH (phytanoyl-CoA Hydroxylase)]. These hub genes were postulated to have associations with a range of aspects, encompassing growth and development, fat deposition, and the levels of amino acids within muscle tissue in buffalo [38]. In the dark-green module, CCDC (Coiled-coil domain containing 88C) and CACNA2D3 (Calcium voltage-gated channel 3) were reported as central hub genes. Calcium voltage-gated channel is regulated by calcium and responsible for muscle contraction [45]. The calcium channels also supply calcium ions to muscles, and these ions serve as secondary messengers, governing cellular functions within muscle tissue [46].

In magenta module, SEPTIN9 and NONO (Non-POU domain containing octamer binding) were identified as the most significant hub genes. Septin 9 is a cell cycle-related gene and indispensable for coordinating myosin motor proteins during cytokinesis [47,48]. In animals, the NONO gene is a multifunctional gene involved in transcription regulation and alternative splicing. Notably, in mice, the NONO gene plays a direct and critical role in B cell development through the ERK (Extracellular Signal-Regulated Kinase), AKT (PI3K/AKT/mTOR), and NF-κB (Nuclear Factor-κB) pathways [49]. This underscores the importance of NONO in the intricate process of cellular and muscle development in animals[49]. Further experiments are necessary to validate the modules, metabolic pathways, and hub genes uncovered in the present study.

5. Conclusions

Identifying genes and metabolic pathways offers tremendous help for marker-assisted selection and plays a key role toward the goal of increasing beef production efficiency. Genes and metabolic pathways highlighted in this study revealed both positive and negative correlations with the beef production related traits measured. Olfactory transduction, hematogenic cell lineage and oxidative phosphorylation were central pathways associated with beef production. Hub genes, including CCDC88C, CACNA2D3, SEPTIN9 and NONO were considered putative candidates related to beef production. The molecular mechanisms underlying beef production warrants further explorations to enhance production efficiency.