We evaluated two groups of divergent animals, with high and low genetic value based on the respective sire’s EBV estimated over 150 bulls and phenotype records on the standard milk yield (SMY) of multi-parous (6 to 7 parities) animals at the government farm of ICAR-Central Institute for Research on Buffaloes, Hisar, India. The EBV for the SMY for each bull comprised the first lactation records of 100 daughters of the respective bull, thus comprising 4.575 million milk records compared with available records for contemporary daughters at the farm. Parental inheritance contribution from bulls is more than 35% toward their daughter buffaloes as per progeny testing, thus expressing high and/or low genetic merits of production as per the inherited potential from their respective sires. All the animals were in mid-lactation physiologically and were kept on a uniform feeding system. Conclusive inference is based on the functional expression of milk yield trait with respect to high and low genetic merit (sire’s EBV and phenotypic records for milk production) in buffaloes. The average SMY over parities was 1869.57 ± 189.36 and 1893.00 ± 275.45 in low milk yielding animals, i.e. MY1 and MY2, respectively. The SMY for high milk yielding individuals was 2909.50 ± 492.63 and 2419.83 ± 361.49 in MY3 and MY4, respectively, as replicates in the two study groups were compared for their RNA-seq data to identify differentially expressed genes(DEGs). RNA-seq data were obtained selectively, considering high and low performance extremes with respect to four economically important traits: residual feed intake, age at first calving, inter-calving period between the first two calvings, and milk volume for growth to develop the “reference” assembly.

RNA extraction was performed using whole blood from B. bubalis. RNA quality was assessed using an Agilent 2100 Bioanalyzer. mRNA was purified using oligo-dT beads (TruSeq RNA Kit, Illumina) with 1 μg of RNA (RIN = 8.0). The pure RNA was fragmented at 90°C and reverse transcribed using random hexamers and Superscript II Reverse Transcriptase (Life Technologies). Complementary DNA (cDNA) was synthesized using RNaseH and DNA polymerase I and was cleaned for SPRI cleanup using Beckman Coulter Agencourt Ampure XP SPRI beads. The cDNA molecules were ligated with Illumina adapters after end repair and were subsequently polyadenylated. The cDNA library was further amplified using PCR to enrich the adapter-ligated fragments. Nucleic acid quantification was performed for individual libraries using a NanoDrop spectrophotometer, and the quality was validated using a Bioanalyzer (Agilent Technologies). Sequencing was performed on an IlluminaHiSeq 2500 platform. Paired-end FASTQ files with Phred scores >20 were used to obtain high-quality (HQ)-filtered reads.

Considering the vast natural diversity with respect to major physiological criteria of genetic merit indicating the economic significance of buffaloes, the “reference” assembly was generated for animals exhibiting extremely low and high performances with respect to early maturity, those achieving first calving at a younger age with the ability to further produce a calf per year, and those taking optimum time for establishing subsequent pregnancy during the reproductive life span as a result of efficient feed utilization selectively to achieve better growth and milk production status.

Therefore, the reference assembly was prepared using paired-end RNA-seq library, as per the protocol recommended by Illumina, (GenBank NCBI databases under SRA ID SUB5692555 and the Bio Project ID: PRJNA546485) including 16 test animals of B.bubalis selected from both low- and high-performance extremes with respect to the following criteria of physiological merit to achieve true assembly covering vast diversity: milk production, age at first calving, inter-calving period, and feed efficiency. Replicates in each of the high- and low-performing groups in each physiological category were sequenced using 100 bp paired end module sequencing and re-assembled into 40–60 million reads for each library. The reads were subjected to quality control using the NGSQC tool kit (Patel and Jain, 2012). Paired end sequence reads with Phred score >Q30 were selected for further analysis.

The HQ-filtered reads obtained from each group were provided as input to a Trinity assembly with an auto-kmer setting with default parameters. The output of the Trinity assembly was subjected to an improved amelioration pipeline to filter putative false positive transcripts using coverage and depth criteria in each library. This resulted in true positive transcripts in each library, which were subjected to CD-HIT EST clustering pipeline for eliminating redundant transcripts, thereby creating a reference unigene transcriptome for each library of B.bubalis. This reference unigene transcriptome was further clustered into a single non-redundant transcriptome assembly used as a “reference” assembly. HQ reads obtained from individual libraries generated from low milk yield performers (MY1 & MY2) and high performers (MY3 & MY4) were aligned with the reference buffalo genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_003121395.1/UOA_WB_1) using the Bowtie tool for quantification of expressed transcripts. Transcripts with RPKM ≥1 were considered as expressed in each of the sequenced library.

The resultant “reference” transcriptome assembly was used as reference, and reads from each individual libraries were aligned to the reference using the Bowtie tool for quantification of expressed transcripts to validate sequencing. Replicate reproducibility was assessed using a correlation matrix (R²). Distribution of inter and intra stages of physiology specifically expressing transcripts (in replicates) were checked via the Venn diagram approach. Read count for each transcript was normalized to obtain RPKM≥1, which qualified it as an expressed transcript in each sequenced library. The final transcriptome was subjected to SNP calling using SAM tools and VarScan and filtered (Depth≥30). SNPs were divided based on the four physiological-stage-specific groups. Milk production specific transcripts and SNPs in the genes were annotated with Gene Ontology (GO) and Pathways to perform biological enrichment analysis.

The reference genome of the Mediterranean water buffalo (B. bubalis) (https://www.ncbi.nlm.nih.gov/assembly/GCF_003121395.1/) was used for read alignment of transcripts, identification of coding regions, and annotation of transcripts using Kallisto (Nicolas et al., 2016). Replicate transcripts obtained from the low- and high-merit subgroups were used for differential expression analysis using DESEQ2 (Love et al., 2014). Transcript log2ratio ≥2.0 was the limit to select differentially expressed transcripts. Volcano plots were generated using the R package “Enhanced Volcano” to identify the differential gene clusters. Gene ontology (GO) and pathways involved in high or low levels of milk production were elucidated using the DAVID Functional Annotation Tool (DAVID Bioinformatics Resources 6.8, NIAID/NIH) (Da et al., 2007). Differentially expressed transcripts obtained from individuals of low and high merit for milk production were identified.

GO and KEGG pathways were identified using the DAVID Functional Annotation Tool. Integrated networks were deduced using the Cytoscape software (Shannon, 2003).

Deep sequencing of RNA obtained from HiSeq platform resulted in an average of approximately 55 million raw reads of 210bp average insert size per sample library ( Table 1 ). The raw FASTQ sequences were filtered using NGSQC tool kit to obtain HQ reads, an average of approximately 50 million HQ reads per sample (~90%). Each library datapoint had an SD of <30% within and across indicating uniformity in data generation, which is critical for estimating transcript coverage and normalization, while the minimum transcript length cutoff of the assembled transcriptome was 500 bases and the maximum transcript length observed was 27 kbp. N50 of the assembled transcriptome was approximately 1.97 kbp, indicating good integrity of the transcriptome available for expression profiling, resulting in identification of an average of approximately 57,000 transcripts expressed in one sample.

Transcript distribution provided the expression dynamics specific to high and low milk yield. Furthermore, the transcripts obtained with respect to differentials between and within high and low milk yield subgroups were annotated using GO and KEGG pathways for biological enrichment analysis. Transcripts obtained from RNA-seq were subjected to differential expression analysis. Data retrieval from the HiSeq platform with respect to the replicates between the high and low milk production subgroups yielded 0.2044 billion reads. An average of 51.1 million raw paired-end reads ( Table 1 ) with an average insert size of 210bp was used to determine DEGs.

The transcriptome assembly comprised 83,783 sequences with an average sequence length of 1518.46 and G + C% of 48.93, based on which replicates were compared. HQ reads obtained from the individual library were mapped over the B. bubalis reference genome, and around 88.7% of the HQ reads were mapped to the human reference genome using KALLISTO pipeline, suggesting an acceptable alignment ( Figure 1 ). Approximately 20% difference in alignment of reads with reference genome indicated underlying diversity in the genetic merit of individual animals under study, as the ratio of HQ reads to total reads obtained in each compared sample library was uniform (0.91–0.93). A high degree of correlation between biological replicates within the performance subgroup, with an average R² value of 0.98 among all the profiled replicates, indicated acceptable phenotype selection of replicates and reliable quality of experimentation and sequencing.

Notable differences were observed among the transcripts of animals of diverse genetic merit referred through the de-regulated genes ( Figure 2 ). Distinct gene patterns obtained from the RNAseq data of low and high milk yield subgroups’ animals depicted through heat maps ( Figure 3 ) covered 30,551 transcripts, representing over 25,000 genes, among which 887 were differentially expressed transcripts (fold-change value <-2 and padj value <0.05). Significantly differentiated transcripts were ascribed to the functional expression of low and high merit in lactating buffaloes.

The downregulated expression of genes was apparently less, including for 373 genes, than the upregulated expression of 514 genes ( Figure 4 ); 83 and 147 genes were significantly(p<0.05) downregulated ( Table 2 ) and upregulated ( Table 3 ), respectively.

Gene enrichment analysis was performed using standard procedures (https://david.ncifcrf.gov/). Significant (p ≤ 0.05) biological processes, cellular locations, and molecular functions that are related to the production of milk volume ( Table 4 ) were selected to plot the networks using System Biology of Milk Yield in Buffaloes. The gene clusters (p<0.05) identified in this study represented 57 categories of annotated GO terms and KEGG pathways govern cell functions, including the cytoplasm, nucleus, nucleolus, and membrane, to propagate milk production ( Figure 5 ). GO:0005654~nucleoplasm (FDR=0.000257993148089497), bta04668: TNF α signaling pathway (FDR=0.00685236886627249), and GO:0016020~membrane (FDR= 0.0814138896939487) were among the most significant (p <0.05) functional categories. The number of genes annotated under each functional category along with the enlisted significant genes are listed in Table 4 . The main biological processes and pathways were zinc ion binding (GO:0008270), negative regulation of extrinsic apoptotic signaling pathway(GO:2001237), negative regulation of intracellular estrogen receptor signaling pathway(GO:0033147), ATP binding (GO:0005524), ubiquitin-mediated proteolysis, UB4B(bta04120), TNF αsignaling (bta:04668), sphingolipid signaling (bta04071), amino acid transmembrane transporter activity (GO:0015171), and RNA splicing regulation (GO:0043484), which are involved in regulating the milk volume ( Figure 5 ). The gene network and protein–protein interaction determined via Cytoscape analysis are depicted in Figure 4 ; Table 4 .

GO:0005654~nucleoplasm emerged as a significant function location ( Table 4 ) while inferring the differences related to high and low milk yield in buffaloes. Other significant (padj<0.05) processes included dystroglycan 1, heterogeneous nuclear ribonucleoprotein F, cleavage and polyadenylation specific factor 6, protein phosphatase 2 phosphatase activator, DExH-box helicase 9, nuclear transcription factor Y subunit alpha, pterin-4 alpha-carbinolamine dehydratase 1, transmembrane and ubiquitin-like domain containing 1, CXXC finger protein 5, and host cell factor C1 regulator 1.TNF signaling pathway (bta04668) emerged as a significant (p< 0.05, FDR<0.05) pathway; the up-regulated genes involved were TRAF1, TRAF2 (TNF ligand), MMP9 (SCC in milk), CCL5 (recruiting leukocytes,T-cells, and macrophages to inflammatory sites), CREB3L1 (homeostasis of the secretary pathway), PIK3R5 (cell growth, proliferation, differentiation, and motility), and LTA (innate immunity). Similarly, down-regulated (p< 0.05) genes were CXCL2 (chemokine antimicrobial protein), BIRC3, ATF2 (inhibits cellular proliferation and induces cell death), RPS6KA5, AKT3, CFLAR (apoptosis regulation), and GRO ( Table 4 ).

TNF α holds significance because it regulates apoptosis, cell survival, inflammation, and immunity as a cytokine, thus propagating the onset of lactation and enhancing the prolactin promoter in the pituitary. The NF-κB signaling pathway (bta04064) corroborates with TNF α signaling (Soünke et al., 2006) in these animals. Estrogen interacts with TNF-α (GO:0033147), regulating PRL expression. Active TNF emerged as a marker protein for milk production in this study. Binding of TRAF2 to TNFR protein molecules may trigger signal transduction to activate genes of the NF-kappa B pathway and the MAPK PI3K-AKT signaling cascade, determining MEC survival and its capacity for milk production. PIK3R5was differentially up-regulated (p<0.05), supporting PI3K-AKT signaling cascades regulating the cytoskeleton and cell organization (GO:0015629)( Figure 4 ). ActivatedPI3K phosphorylates AKT and regulates apoptosis, protein synthesis, metabolism, and the cell cycle. PI3K executes actin reorganization and cellular polarization through G-protein-coupled receptors and the secretion of chemokines and signaling proteins (GO:0008009 and bta04062).The PI3K-Akt-mTOR-signaling pathway thus creates a structural link between membrane receptors and the actin cytoskeleton, harboring signaling molecules (kinases and phosphatases) to regulate gene functions and cell survival. Associated pathways initiate the innate immune response against pathogens via Toll-like receptors (TLRs), producing pro-inflammatory cytokines using activated NF-κB and MAPK pathways.TLR ligand binding triggers responses, such as protein modifiers, helper cofactors, and post-transcriptional/epigenetic regulators.

Baculoviral IAP repeat-containing gene (BIRC3) emerged as a significantly downregulated gene in this study. BIRC3 regulates caspases, apoptosis, inflammatory signaling, immunity, mitogenic kinase signaling, cell proliferation, and cell invasion. BIRC3 acts as an E3 ubiquitin-protein ligase, regulating the NF-kappa-B signaling. The present study confirmed the collaborative role of metabolic and immune regulatory pathways to establish the genetic merit of milk production, thus corroborating with a previous hypothesis (Bu et al., 2017) regarding mutual crosstalk between the liver and mammary tissue for optimum lactation. Nevertheless, the whole blood transcriptome, a non-invasive resource for study, enabled projections of lactation physiology corroborative with the genes and functions that emerged as lactation enhancers with potential as markers for selecting high-yielding (high EBV) buffaloes.

In the present study, whole blood transcriptome reflected differential gene expression with respect to buffaloes with high and low genetic merit for milk production. Identification of DEGs using non-invasive methods is useful in large dairy animals. Most transcripts were associated with transcriptional regulation, primarily RNA synthesis, mRNA translation, cell signaling, and proliferative capacity of mammary cells either regulated by hormones and/or growth factors, whereas cell metabolic activity is controlled by the liver and is thus regulated by diet (Bu et al., 2017). These findings indicate the major role of transcription factors in the nucleus as regulators for the repression or activation of biological processes of cell/ribosomes, and epigenetic regulation by amino acids and amino acid transporters as inducers and availability of energy (NADH) across the membranes which are the key factors to achieve adequate milk production, consistent with previous reports (Dai et al., 2018). Other than the structural variation depicted as genomic variants reported earlier (Mishra et al., 2021). Pre-mRNA modifications introduce expression variants at the level of individual animals. Differential gene expression reflects the physiological differences in mammary tissue functioning between high and low milk producing buffaloes. A decrease in methyl-deoxycytidine levels in buffaloes is associated with increased milk production due to de-condensation of chromatin, consistent with previous reports (Choi et al., 1998). Casein production is enhanced by reduced cytidine methylation, indicating differential epigenetic regulation in high and low milk producing animals.

The present study revealed the impact of cell growth and differentiation factors, cytokines/chemokines, and transcription repressors as signaling molecules on epithelial cell differentiation and proliferation, migration, vasculogenesis, apoptosis, innate immune cells, inflammation (i.e., the NFκB pathway), and translocation of signal carriers from external membranes to nuclear regulatory factors through receptors/signal sensors in the mammary tissue. Enrichment of these functions indicate that lactation results from the combined reciprocal control of signals induced as metabolites/molecules originating from the pituitary, liver, and mammary tissue in synergy with immune function regulation in buffaloes.