BACKGROUND
Milk fat is a major determinant of milk quality and pricing, and breed-related differences in milk composition are well recognized. The authors note that indigenous cattle in India vary widely in milk quantity and quality, making them a useful model for studying the genetics of milk traits. Prior work has identified important bovine milk QTL genes such as `DGAT1`, `GHR`, `ABCG2`, and `SPP1`, but the authors emphasize that no previous study had specifically evaluated variation in fat QTLs across indigenous breeds by integrating gene expression and whole-genome variant analysis. The goal of this study was therefore to identify genomic variation within milk-fat-associated QTL genes that were also differentially expressed, comparing high-milk-yield and low-milk-yield cattle breeds.
METHODS
The study first mined the Animal QTLdb for genes associated with milk fat yield and milk fat percentage. The authors identified `286` genes for milk fat yield and `256` genes for milk fat percentage, yielding `417` unique genes overall and `125` genes common to both traits after duplicate removal. Protein interaction network analysis of the `417` genes was performed using STRING 11.0 and Cytoscape 3.8.0, and `74` genes were designated as hub and/or bottleneck genes. In parallel, publicly available milk transcriptome data from bioproject `PRJNA419906` were analyzed. This dataset included mammary epithelial cell RNA-seq from `six` lactating cows, with `three Jersey` cows and `three Kashmiri` cows sampled on Day `15`, Day `90`, and Day `250`, representing early, mid-, and late lactation. Jersey was used as a higher-fat breed with fat content ranging from `4.10–4.86%`, and Kashmiri as a lower-fat breed with fat content ranging from `3.20–3.94%`. Reads were quality controlled and mapped to the `Bos indicus` genome, and differential expression was assessed using DESeq2.
For genomic analysis, the authors studied `14` indigenous cattle divided into high- and low-milk-yield groups. The high-milk-yield group had average milk yield per day of `8 kg` and included Sahiwal (`n = 4`) and Gir (`n = 4`). The low-milk-yield group had average milk yield per day of `2.5 kg` and included `6` animals representing Gaolao, Deoni, Hallikar, Dangi, Pulikulam, and Amritmahal. Whole-genome sequencing libraries were generated on an Illumina platform for paired-end sequencing (`2150 bp`, as reported). Variant calling was performed with freebayes and GATK, SNPs were filtered, and common SNPs across callers were extracted. Non-synonymous SNPs were annotated using SnpEff. Three nsSNPs in `GHR`, `LPIN1`, and `TLR4` were selected for validation by pyrosequencing in `14` animals.
KEY RESULTS
Functional enrichment of the `125` genes common to milk fat yield and percentage showed involvement in biosynthetic, catabolic, regulatory, transportation, and cellular-response metabolic processes. The protein interaction network contained `403` nodes and `671` edges. Based on network ranking, `50` hub genes and `50` bottleneck genes were selected, with `74` total unique genes classified as hub and/or bottleneck. Of these, `25` genes were differentially expressed between Jersey and Kashmiri, accounting for `18` hubs and `17` bottleneck genes; `10` genes had both hub and bottleneck characteristics.
RNA-seq analysis showed that `70` genes were upregulated and `52` genes were downregulated in Jersey compared with Kashmiri. Among the top network-ranked genes, `SRC` had the highest degree of association at `30` with a `log2` fold change of `1.480587`, while `DGAT1` had a degree of `25` with a `log2FC` of `0.921104`. The discussion also highlights `DGKG` as highly upregulated with `log2FC = 4.03`, `SLC6A9` with `log2FC = 4.49`, `UGDH` with `log2FC = 0.75`, and `IDH1` as downregulated with `log2FC = −1.632`.
Whole-genome sequencing of the `14` samples generated `12.77 billion` paired-end reads. After preprocessing, `11.02 billion` clean reads remained, corresponding to approximately `1516 Gb` of data. Each dataset had a minimum sequencing depth of `≥30x`, average GC content of `45.26%`, mean `97.91%` Q20 bases, and mean `93.97%` Q30 bases. Alignment to the Brahman reference genome exceeded `>95%`. Initial variant calling produced `63,357,363` variants. After filtering on `Q > 20`, `33,976,892` SNPs were identified. GATK analysis yielded `39,625,917` variants that passed recalibration, and `25,956,231` SNPs were common to both variant callers.
From the `25` differentially expressed hub/bottleneck milk fat QTL genes, `20` genes had non-synonymous substitutions in coding regions. A fixed SNP pattern in high-milk-yield breeds relative to variable patterns in low-milk-yield breeds was seen in `GHR`, `TLR4`, `LPIN1`, `CACNA1C`, `ZBTB16`, `ITGA1`, `ANK1`, and `NTG5E`. The reverse pattern was observed in `MFGE8`, `FGF2`, `TLR4`, `LPIN1`, `NUP98`, `PTK2`, `ZTB16`, `DDIT3`, and `NT5E`. Pyrosequencing confirmed SNPs `C/G`, `C/A`, and `G/A` in `GHR`, `TLR4`, and `LPIN1` in the low-milk-yield breeds Amritmahal, Pulikulam, and Dangi, versus `C/C`, `C/C`, and `G/G` in the high-milk-yield breeds Gir and Sahiwal. Specific examples included `TLR4` variant `g.107083326A>C` in the low-milk-yield group but fixed in the high-milk-yield group, `LPIN1` variant `g.85211528C>G` in the low-milk-yield group but fixed in the high-milk-yield group, and `NUP98` variant `g.32707374G>A` plus `LPIN1` variant `g.85205642T>G` in the high-milk-yield group while fixed as `g.32707374G>G` and `g.85205642T>T` in the low-milk-yield group. `LPIN1` and `ITGA1` had the highest SNP count at `10`, whereas `PTK2`, `IGF1R`, `DDIT3`, `CXCL8`, and `LPL` had the lowest SNP counts.
CLINICAL IMPLICATIONS
This is not a clinical treatment study, but it has practical importance for animal breeding and dairy production. The data suggest that milk-fat-associated biology differs between high- and low-milk-yield cattle not only at the transcriptome level but also through potentially functional coding variation. The identified nsSNP patterns in genes such as `GHR`, `TLR4`, and `LPIN1` provide plausible marker candidates for future validation in larger breed panels. However, the authors appropriately conclude that the functional impact of these SNPs on fat yield, fat percentage, or overall milk yield still needs further study before they can be used reliably in selection programs.