Integrative Analysis of Transcriptome and GWAS Data to Identify the Hub Genes Associated With Milk Yield Trait in Buffalo

Abstract

The mammary gland is the production organ in mammals that is of great importance for milk production and quality. However, characterization of the buffalo mammary gland transcriptome and identification of the valuable candidate genes that affect milk production is limited. Here, we performed the differential expressed genes (DEGs) analysis of mammary gland tissue on day 7, 50, 140, and 280 after calving and conducted gene-based genome-wide association studies (GWAS) of milk yield in 935 Mediterranean buffaloes. We then employed weighted gene co-expression network analysis (WGCNA) to identify specific modules and hub genes related to milk yield based on gene expression profiles and GWAS data. The results of the DEGs analysis showed that a total of 1,420 DEGs were detected across different lactation points. In the gene-based analysis, 976 genes were found to have genome-wide association (P ≤ 0.05) that could be defined as the nominally significant GWAS geneset (NSGG), 9 of which were suggestively associated with milk yield (P < 10^-4). Using the WGCNA analysis, 544 and 225 genes associated with milk yield in the turquoise module were identified from DEGs and NSGG datasets, respectively. Several genes (including BNIPL, TUBA1C, C2CD4B, DCP1B, MAP3K5, PDCD11, SRGAP1, GDPD5, BARX2, SCARA3, CTU2, and RPL27A) were identified and considered as the hub genes because they were involved in multiple pathways related to milk production. Our findings provide an insight into the dynamic characterization of the buffalo mammary gland transcriptome, and these potential candidate genes may be valuable for future functional characterization of the buffalo mammary gland.

Keywords: RNA-seq; WGCNA; buffalo; genome-wide association studies; hub genes; milk yield.

FIGURE 1

Expression profiles of buffalo mammary gland tissues during four lactation point. (A) The cluster heat map of all mRNAs expression at different stages during lactation. (B) Bar plots are showing top10 genes with highest TPM values across the four points. (C) Bar plots showing differentially expressed genes by the pairwise comparison. (D) Hierarchical clustering is showing often up and down-regulated mRNAs across the four points.

FIGURE 2

Manhattan plot (A) of -log10 (P-values) and Quantile-Quantile plot (B) of P-values for milk yield from the gene-based method. The blue horizontal line indicates the suggestive significance level [-log10(1e-4)].

FIGURE 3

Identification of modules and functional annotation analysis for the module genes. (A) Module detection for the DEGs dataset and GO analysis for the module genes. (B) Module detection for the NSGG dataset and GO analysis for the module genes. (C) KEGG enrichment analysis for module genes from the DEGs dataset. (D) KEGG enrichment analysis for module genes from the NSGG dataset.

FIGURE 4

Heatmap of the correlation between module eigengenes and milk yield. (A) module-traits analysis for the DEGs dataset. (B) module-traits analysis for the NSGG dataset.

FIGURE 5

Hub genes detection and network construction analysis. (A) Scatter plot of module eigengenes in the turquoise module from DEGs dataset. (B) Scatter plot of module eigengenes in the turquoise module from NSGG dataset. (C) The Venn diagram of the DEGs and NSGG hub genes. (D) Hub gene interaction network of in the turquoise module from the DEGs and NSGG dataset. The color intensity in each node was proportional to the TOM values calculated by WGCNA (the higher TOM values were in a circle with red, whereas the lower TOM values were in a circle with white).