Identification of single nucleotide polymorphisms (SNPs) potentially associated with residual feed intake in Qinchuan beef cattle by hypothalamus and duodenum RNA-Seq data

Experimental animals and data collection

Based on our previous study (Yang et al., 2021), 30 Qinchuan bulls with similar age (15 ± 1 months) and weight (280.6 ± 30.9 kg) were selected from a farm in Ningxia, China. The study subjects were given a standardized feeding regimen throughout the experimental period, and free access to water and food was ensured. Body weight measurements were taken monthly throughout the 81-day experimental period, then daily feed intake, average daily gain (ADG), and the midpoint metabolic body weight (MMBW0.75) was calculated based on feed intake (FI). To classify the cattle into high and low RFI groups, we used multiple linear regression of FI on the midpoint MMBW^0.75 and ADG to estimate individual RFI (Yang et al., 2021).

RNA extraction and sequencing

Based on the results of the RFI calculation, five individuals with extremely low RFI (LRFI, high efficiency) and high RFI (HRFI, low efficiency) phenotypes were selected for slaughter after a 16-h fasting period. All experimental procedures involving animals were conducted by the Guidelines for Ethical Review of Laboratory Animal Welfare of Ningxia University (NXUC20211015). The hypothalamus (including the arcuate strong nucleus, parabrachial nucleus, supraoptic nucleus, dorsal/ventral medial nucleus, and other brain tissues) and the descending duodenum (mucosa, submucosa, and external muscular propria) were collected post-slaughter. Our tissue samples include hypothalamus and duodenum tissues from five high RFI and five low RFI cattle, totaling 20 samples. These tissue samples were washed with PBS to remove blood and other impurities, then cut into small pieces to increase surface area, and placed into sterile tubes for further processing. Total RNA was extracted from 500 mg of tissue samples using TRIzol method (TaKaRa Bio, Beijing, China), following the manufacturer’s instructions. After extraction, the RNA was further purified using column purification to enhance purity and remove residual contaminants. DNase treatment was performed to eliminate genomic DNA contamination. The quality and integrity of the extracted RNA were assessed using 1% agarose gel electrophoresis, Nanodrop, and Agilent 2100 to ensure a sample concentration of ≥500 ng/µL, 28S:18S > 1.0, and RIN ≥ 7. For library construction, the RNA was first reverse transcribed into cDNA. Adapter ligation was then performed to facilitate the sequencing process, followed by amplification to enrich the library. The library’s initial quantification was carried out using Qubit 2.0, and the insert size was verified using an Agilent 2100. The effective concentration of the library (effective concentration > 2 nM) was accurately determined using qRT-PCR. The hypothalamus tissue samples underwent whole transcriptome sequencing to capture both coding and non-coding RNA species, providing a comprehensive view of the transcriptome. The duodenum tissue samples were subjected to standard transcriptome sequencing. Finally, pair-end sequencing data (raw data) with 150 bp read length were generated using the Illumina HiSeq 4000 platform. High RFI sequencing data of hypothalamus and duodenum were named Q_H1~Q_H5 and S_H1~S_H5, respectively, while low RFI sequencing data were named Q_L1~Q_L5 and S_L1~S_L5, respectively.

Quality control, mapping and transcript assembly

The statistical power of this experimental design, calculated in RNASeqPower (https://bioconductor.org/packages/release/bioc/html/RNASeqPower.html) is 0.86 (The corresponding code can be found in Script S1). The quality of raw data was assessed using the fastQC v.0.11.9. Subsequently, the Trimmomatic v.0.39 was used to perform quality control on the data. This included removing adapter sequences, trimming bases with Phred scores below 30 at the beginning and end of reads, applying a sliding window approach with a window size of 5 bp to remove bases with an average Phred score below 20, and discarding reads shorter than 75 bp. The cleaned data were then reevaluated using the fastQC software to ensure they met the requirements for subsequent analysis. The clean data were aligned to the bovine reference genome (ARS-UCD1.2; INSDC Assembly) using the STAR v.2.7.3a with the following parameters: --outSAMtype BAM Unsorted SortedByCoordinate, -outFilterMismatchNmax 999, --outFilterMismatchNoverReadLmax 0.04, --outFilterMultimapNmax 1. The resulting alignment files were further processed using the AddOrReplaceReadGroups tool of the PICARD v.2.27.4. This added the sample ReadGroups (RG) information to each alignment file. Additionally, the MarkDuplicates tool was applied to remove duplicate amplifications resulting from the PCR process during library construction.

Merging of sample data

For increasing the number of reads per variant locus, enhancing the depth coverage of reads across the entire transcriptome, as well as the depth coverage and quality of variant calls (Lam et al., 2020), in this study, the data from hypothalamic and duodenal tissues were merged into two BAM files based on phenotype (high RFI group and low RFI group). This merging process, performed using the “merge” command of the samtools v.1.16.1, aimed to balance sequencing depth between samples and minimize the impact on SNP analysis results. Both the high RFI and low RFI groups in the subsequent analysis referred to the combined group data (Fig. 1).

SNPs recognition, filtering and annotation

BCFtools v.1.16 was utilized to execute the variant calling on the combined data of the high RFI group and low RFI group respectively, enabling identification of SNP sites and generating BCF files containing variant information. The “norm” parameter of BCFtools was then employed to normalize the variant information, thereby eliminating ambiguity caused by varying methods. Subsequently, the low-quality SNPs data underwent further filtering to reduce the likelihood of false positives and alleviate computational resource requirements for subsequent analysis. The software BCFtools and VCFtools v.0.1.16) were employed for variant filtering, employing the following criteria: (1) Removal of SNPs within a 5 bp range near indels; (2) setting a minimum coverage (DP) of 10; (3) enforcing a minimum allele frequency not less than 0.2 and a secondary allele depth not less than 2; (4) filtering loci with quality scores below 30. Finally, the functional annotation of SNPs was performed using the SnpEff v.5.1d with the built-in ARS-UCD1.2.105 database. The thresholds of above software are referenced from previous study (Lam et al., 2020).

Identifying and annotating high and low RFI group-specific SNPs

Using the SnpEff software, the VCF files underwent annotation, allowing for the identification of SNPs specific to the high and low RFI groups. SnpSift v.5.1d was then employed to screen for SNP loci with significant functional and modifier-type impacts. This enabled the selection of candidate genes associated with these SNP loci.

Gene function enrichment analysis and protein interaction network analysis

We employ clusterProfiler to conduct GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis. The filter parameter is set as p value < 0.05. The GO enrichment can further clarify the main biological functions of the genes where the specific SNPs are located. The KEGG pathway enrichment can be used to understand the signal pathway regulated by the genes. Using the STRING database (https://cn.string-db.org/), we perform protein interaction analysis on the relevant genes to select core genes that have interaction effects.

Statistical analysis

The data from the experiment were analyzed and visualized using R v.4.3.0 (R Core Team, 2023) and Prism v.10.1.1 software. Statistical significance between the treatment and control groups was assessed using non-parametric tests or t-tests. A p-value of less than 0.05 was considered indicative of a significant difference. The corresponding code scripts can be found in Script S1.

RNA-seq sequencing data quality and comparisons

From the RNA-seq data of the hypothalamus and duodenum, we obtained a total of 1,019 million and 275 million paired sequencing reads, respectively. After quality control, all 20 samples in this study had a Q30 score (error rate p < 0.001) above 96%. The GC% content was approximately 50% (Table S1). The clean data was aligned to the bovine reference genome, where the percentage of reads aligning to the reference genome was above 91%, and the percentage of reads with a unique alignment position ranged from 81.07% to 93.81% (Table S2). Analysis of the transcript expression for each sample indicated relatively consistent transcript abundances (Fig. 2). These results demonstrate the high quality of the obtained data, reducing the impact of sequencing errors on subsequent analysis.

SNPs screening and analysis

SNPs in the high and low RFI groups were identified by combining RFI phenotype with RNA-Seq data from hypothalamic and duodenal tissues. The numbers of homozygous and heterozygous mutations in the combined samples were counted (Table 1), revealing that there were 270,410 and 255,120 specific SNPs in the high and low RFI groups, respectively. Among these SNPs, 11,991 (4.43%) and 14,007 (5.49%) were homozygous mutations, with heterozygous mutations far outnumbering homozygous mutations. Statistical analysis based on the different types of mutations in the SNPs (Fig. 3) showed that the total number of transition types (A–G, C–T) was higher than the total number of transversion types (A–T, C–G, A–C, G–T). Among the six types of single nucleotide variations, A–G and C–T had the highest occurrence rates, while the occurrence rates of the other four types of transitions were relatively lower. In the high RFI group, A–G accounted for 36.5% and C–T accounted for 35.4%. In the low RFI group, A–G accounted for 35.7% and C-T accounted for 35.8%. The Ts/Tv ratios in the high and low RFI groups were 2.55 and 2.50, respectively. The differences in occurrence rates of the six types of single nucleotide variations between the high and low RFI groups were small, with a mean occurrence rate of transitions being 71.72% and transversions being 28.29%.

Distribution statistics of SNPs

The study investigated the distribution and variation of SNPs across different chromosomes, providing insights into the genetic diversity among genes. Analysis of the combined high and low RFI groups revealed no significant difference (p > 0.05) in distribution between the two groups. Chromosome 1 exhibited the highest number of SNPs, while the mitochondria (MT) showed the lowest distribution (Fig. 4A). To account for differences in chromosome length, the ratio of SNPs number to chromosome length was calculated, revealing that the MT had the highest variant rate among both groups, indicating a higher density of SNPs per unit length (Fig. 4B). The chromosomal distribution of SNPs within the high and low RFI groups is delineated in Figs. 4C and 4D.

Comprehensive statistical analysis of SNP loci distribution across the genome was conducted for both high and low RFI datasets, emphasizing different genomic functional regions, such as downstream, exon, intergenic, intron, untranslated region, etc. (Table 2). The statistical analysis of the SNP locations in the genome for the high RFI and low RFI groups shows ten different distributions (DOWNSTREAM, EXON, INTERGENIC, INTRON, SPLICE_SITE_ACCEPTOR, SPLICE_SITE_DONOR, SPLICE_REGION, UPSTREAM, UTR_3_PRIME, UTR_5_PRIME) (Table 2). A single SNP may be located in multiple transcript regions. The analysis found that the INTRON region had the most SNP locations in both groups, with 429,995 areas annotated in the HRFI group and 435,881 areas in the LRFI group, significantly more than other functional regions. The next most abundant functional elements are INTERGENIC and DOWNSTREAM, while the remaining functional regions are less common. SPLICE_SITE_ACCEPTOR has the fewest functional regions, with 37 in the high RFI group and 73 in the low RFI group. We hope to find most of the SNPs in the exonic regions, but coding regions generally experience higher selective pressure compared to non-coding regions (Zhao et al., 2003). The annotation of SNPs in the high and low RFI groups accounts for 1.11% and 1.38% in the exonic regions, respectively. At the same time, this also explains our detection results: the higher distribution of SNPs in intron regions is partly due to the fact that unspliced transcripts are also detected during sequencing, and partly because intron regions constitute 47.51% of the whole genome, which is significantly higher than the length of exonic regions (Fig. 4E). SNPs located in intergenic regions may be found in new genes or gene portions that have not been annotated yet.

Influence prediction and amino acid change

By using the SnpEff software to evaluate the potential effects of SNP mutations on codons, it was found that over 98% of the SNPs in both groups were classified as modifiers, having minimal effect on genes and proteins. However, 159 SNPs in the high RFI group and 293 SNPs in the low RFI group were predicted to have a high effect (Figs. 5A, 5B; Tables S3, S4), warranting further investigation. This situation is as we expected, most SNPs are located in intron regions and intergenic areas, making it difficult to directly affect protein coding. Therefore, high-impact SNPs will be relatively fewer. Additionally, analyzing the overall levels of each amino acid can provide important insights into evolutionary pressures and adaptation mechanisms. This approach helps identify patterns and frequencies of amino acid substitutions in different biological contexts, enhancing our understanding of how these changes affect protein stability, function, and interactions. To analyze the potential impact of SNPs on genes and proteins, the effect of intergroup-specific SNPs (Tables S5, S6) on codons and subsequent amino acids was assessed. The analysis revealed that the amino acids most affected in both the high and low RFI groups were alanine-threonine, alanine-valine, and isoleucine-valine (Figs. 5C, 5D). By identifying SNP loci that significantly impacted both the high and low RFI groups and mapping them to the corresponding genes using the SNPs annotation files, a total of 83 genes were identified in the high RFI group and 97 genes in the low RFI group. Interestingly, one gene, JSP.1, belonging to the MHC class I family, was common in both groups, and played a key role in regulating animal health within the immune system (Hewitt, 2003).

Genes function annotation of high-impact SNP loci

GO functional annotation and enrichment analysis were conducted for the aforementioned genes (Figs. 6A, 6B). The results revealed that the enriched genes in the high and low RFI groups were primarily associated with protein binding and enzyme binding processes. Notably, a significant number of genes related to NADH activity were found in the low RFI (high feed efficiency) group. These genes were associated with oxidoreductase activity, acting on NADH or NADPH; NADH dehydrogenase activity; NADH dehydrogenase (ubiquinone) activity; NADH dehydrogenase (quinone) activity; and oxidoreductase activity, acting on NADH or NADPH, quinone, or similar compounds as acceptors. NADH and its phosphorylation product NADPH play pivotal roles as coenzymes in various metabolic activities, including cell signaling, protein modification, energy metabolism, mitochondrial function, calcium homeostasis, antioxidative stress, biosynthesis, and cell death (Berger, Ramirez-Hernandez & Ziegler, 2004; Patterson et al., 2005; Xiao et al., 2018; Ying, 2006, 2007, 2008). In particular, the enrichment genes ND4, ND5, and ND6 are core subunits of the mitochondrial respiratory chain NADH dehydrogenase (complex I). They facilitate the transfer of electrons from NADH through the respiratory chain, utilizing ubiquinone as an electron acceptor, and are crucial for the catalysis and assembly of complex I (UniProt, 2023).

The KEGG enrichment analysis outcomes indicated enrichment of pathways related to thyroid hormone synthesis, pancreatic secretion, gastric acid secretion, cAMP signaling pathway, thermogenesis, parathyroid hormone synthesis and secretion, glycerolipid metabolism, TNF signaling pathway and beta-alanine metabolism in the high and low RFI groups (Figs. 6C, 6D). Notably, the thermogenesis pathway exhibited enrichment of ND4, ND5, and ND6 genes. Additionally, ATP1A2, SLC9A4, and PLA2G5 were identified as genes associated with energy metabolism (Lingrel, 1992; Sakuta et al., 2020; Sun et al., 2004).

Protein-protein interaction analysis of high-impact SNP loci

Protein-protein interaction analysis is a method used to study the interactions between proteins, which can be employed to uncover the relationships and networks among proteins, consequently explaining functional interactions and illustrating the intricate interconnections between proteins. Our results revealed distinct patterns of core genes and interaction relationships between the two groups. In the high RFI group, we identified 29 core genes and 23 interaction relationships, while in the low RFI group, we found 42 core genes and 41 interaction relationships (Fig. 7). Several genes, such as HSP90AA1, EIF2AK3, PAK1, MAP3K7, PGM2L1, DNM1L and CYB5R3, were found to be related to energy metabolism, fat deposition and muscle development (Badri et al., 2018; Charoensook et al., 2012; Chen et al., 2019; Chiang & Jin, 2014; Hogarth et al., 2018; Liu et al., 2024; Lopez-Bellon et al., 2022; Zhang, O’Keefe & Jonason, 2017).

Analysis of candidate gene SNP loci

Based on the results of GO and KEGG analysis, we focused on phenotype-related terms. In the high RFI group we screened GO terms: positive regulation of metabolic process and multicellular organismal development; KEGG terms: thyroid hormone synthesis, gastric acid secretion, cAMP signaling pathway, pancreas signaling pathway. In the low RFI group we screened GO terms: oxidoreductase activity, acting on NADH or NADPH; KEGG terms: thermogenesis, metabolic pathways, TNF signaling pathway. Finally, 18 genes were identified in the high RFI group and 21 genes in the LRFI group (Tables S3, S4, S7). Also combining protein interaction analysis and existing studies, we finally focused on 14 genes. In these genes, we combined the prediction results of SNP impact, showing high-impact SNP sites in the genes. We found that these SNPs are mostly located in the exon regions, and most are A-G mutation types. This type of mutation might change the coded amino acids, and affect the structure and function of the protein (Table 3). However, their specific function and roles would be detected in the future studies for clarifying their variant effect to phenotype.