Epigenetic Inheritance of Maternal Disease Status in Dairy Cattle and Effect on Offspring Performance

Progress Report Results:

Biological Samples and Epigenetic Data Progress (January 2023): As of January 6, 2023, 39 Holstein heifer and 5 bull calves have been successfully enrolled in the study marked by successful blood collection from both calf and dam within 24 hours of birth. Of the calves enrolled, 6 heifer calves and all bulls were sold and are no longer part of the study. One heifer calf was euthanized. Of the 32 calves that have remained in the herd, 21 were their dam’s first calf, 10 were born to multiparous dam’s with no history of transition disease in the previous lactation, and 1 was born to a dam previously treated for hypocalcemia. They have been sired by 8 different sires. 7 calves have been treated for scours and 1 was treated for a navel infection. 14 have been successfully weaned and are being transferred to the heifer raiser. The remaining heifer calves are still at the farm and are being monitored for growth and disease events until weaning.

We do have concerns about the lack of calves being born to dams with a history of transition disease. Historically, this has not been a problem, but recently breeding strategies of the VTD herd have placed more emphasis on sourcing replacement animals from primiparous animals. Most heifers are bred with sexed semen; whereas, most multiparous animals are bred with conventional dairy semen or beef semen. Furthermore, the likelihood that an animal is bred with beef semen likely increases if she is a “problem cow” that has a history of disease events. The targeted use of sexed semen and beef semen is a progressive breeding strategy allowing herds to maximize genetic progress and gain a higher sale price for calves not intended as replacements. For this study, it is a problem because there are a small proportion of Holstein heifer calves that could qualify for objectives 2-4. We are considering partnering with a larger commercial herd for sourcing the biological samples for this part of the project.

Tentatively, enrollment of all Holstein heifer calves born at the VTD will continue until June 2023.  All calves enrolled will be followed to their first calving (beginning spring of 2024). At 4 weeks into their lactations, we will choose 12 heifers based on their own and their dam’s disease status (as previously described in the proposal). An additional blood sample at 4 weeks into the lactation will be collected from these animals. Birth and lactation blood samples of the 12 selected dams and calves will be analyzed by RRBS DNA methylation sequencing. The sequences will be used for Objectives 2-4.       

Objective 1 Progress (January 2023): The dam disease status for all current cows in the herd is being extracted from the herd’s archived DairyComp records. To date, all current first and second lactation cows have been classified as either born to a primiparous dam (excluded from study), born to a dam with no record of clinical disease in the lactation coinciding with their gestation (“healthy dam” cohort), or born to a dam with recorded clinical disease in the lactation coinciding with their gestation (“diseased dam” cohort). Dam disease events have included mastitis, metritis, retained placenta, ketosis, displaced abomasum, and hypocalcemia and have not yet been filtered by timing in lactation. First lactation production and disease events have also been extracted for all current first and second lactation animals in the herd. Third and greater lactation animals currently in the herd will be classified by dam disease status and their archived first lactation records will be extracted in the coming months. Statistical analysis of this data will then be conducted to investigate the relationship between dam disease history and offspring performance.

Objective 2 – 4 Progress (Jan 2023 Report): DMR analysis will commence once study subjects have been selected and RRBS sequences have been received (estimated summer 2024).  

Biological Samples and Epigenetic Data Progress (January 2024 Report):

All sampling for the adjusted study design has been completed as of November 2023. All heifers due to calve starting in July 2023 were sampled 1x a week starting two weeks before their due date until they calved. They were then sampled again at 30 ± 4 DIM. This was done routinely until the goal of 6 matched pairs were enrolled. In total, 24 heifers were enrolled making eight case-control pairs. Six of the cases had ketosis, one had ketosis and metritis, and one had metritis. The six pairs with ketosis have been chosen for the study.

The multiparous dam case-control pairs were selected in July 2023 using the farm’s electronic records. Eight were chosen as a precaution in case there were stillborn calves and one calf was stillborn. Of the six pairs of calves successfully enrolled, 1 calf was born to a dam with a history of ketosis and the other 5 case calves had dams with a history of metritis. 8 calves were male and 4 were female. The female calves who were selected to remain in the herd had an additional blood collection between 36-48 hours after birth for serum total protein reading.

DNA has been isolated and stored for all blood samples collected for the project. The samples needed for sequencing have been selected (36 in total) and were sent to Novogene for sequencing in early January.

Final Results and Discussion (December 2024): 

Objective 1: Investigate associations between dam history of disease and offspring’s performance using a retrospective study design

Descriptive Statistics:

A total of 4,605 unique offspring records were extracted from electronic records of both farms after removing animals who did not have a dam listed in the records or who had a dam who could not be classified by her disease status (Fig. 4). These records spanned 12 years and overall CURC had more animal records than the VTD. In both farms there were more offspring from multiparous dams than primiparous dams in total, but the number of offspring from primiparous dams increased in recent years in both farms (Fig. 5). In the CURC dataset, there were more offspring from dams with a history of disease in total and that seems to be due to an increase in offspring from these dams in more recent years (Fig. 6). In the VTD dataset, the opposite was true; overall there were less offspring to dams with a history of disease and there was a decreasing trend in the number of offspring from diseased dams in recent years (Fig. 6). 

Disease outcomes were separated into two timeframes for assessing offspring performance: from birth until 1^st calving and in her 1^st lactation. Disease events were grouped into digestive, respiratory, and other diseases from birth to 1^st calving and into digestive, metabolic, uterine, mastitis, and other diseases within 1^st lactation. The number of cases of each disease type fluctuated from year to year and across farms (Fig. 7) but overall, there was a total disease prevalence of 57% from birth to 1^st lactation and 44% in 1^st lactation. Culling and milk production were also evaluated as offspring performance outcomes. Overall, 35% of animals left from birth to 1^st lactation and 13% of animals that made it to 1^st lactation left in 1^st lactation. Standardized 305 day mature equivalent milk production was numerically higher in CURC animals than VTD animals (Fig. 8).

Effect of Dam Parity on Offspring performance:

All animals were classified based on their dam’s parity. All animals who were her dam’s first calf were in the Primi-Dam group and animals who were her dam’s second or greater calf were grouped in Multi-Dam. The effect of dam parity on offspring milk production and disease risk was evaluated (Table 1). The calves from the Multi-Dam group were less likely to have disease (OR = 0.80, p-value = 0.02), specifically respiratory disease (OR = 0.86, p-value = 0.07) from birth to first lactation than the calves from the Primi-Dam group. Of the offspring that made it to first lactation, offspring of multiparous dams were more likely to have disease (OR = 1.22, p-value = 0.04) specifically metabolic disease (OR=1.34, p-value=0.08). There was no significant difference in milk production between offspring of differing dam parities. There was no significant difference in lifetime survival probability between Multi-Dam and Primi-Dam offspring.

Effect of Dam Disease History within dam parity groups:

Within dam parity groups, all animals were classified based on their dam’s disease history in the period preceding their birth. Within the Primi-Dam group, the offspring was classified based on the dam’s disease records from her birth to her 1^st calving and placed in the Primi D if there was any disease reported or in Primi NoD if there was no disease. Within the Multi-Dam group, offspring were placed in the Multi D group if the dam had any disease reported in the lactation before the offspring was born or Multi NoD if there was no disease reported in that lactation. The effect of dam disease status on offspring milk production and disease risk was evaluated (Table 2 and 3). There were no significant differences between Primi D and NoD offspring from birth to first calving. In first lactation, Primi NoD offspring tended to be less likely to have mastitis (OR=0.68, p-value = 0.06).  Across their whole lives, offspring of Primi NoD were at less risk of leaving the herd than offspring of Primi D (HR=0.88, p-value = 0.04; Fig. 9). There were no significant differences in disease, milk production, or culling outcomes between Multi D and Multi NoD offspring.  

Discussion and Limitations:

In our study we saw a bigger effect of dam parity on offspring performance than dam disease status. Similarly to Carvalho, et al. we saw offspring of multiparous dams were more likely to have disease in 1^st lactation than offspring of primiparous dams and that there was no significant effect of dam disease status on offspring milk production. Carvalho, et al. observed a lower risk for disease in offspring whose dam had a history of disease. Based on the Carvalho Study, we hypothesized that there was a protective effect of dam disease history on offspring disease susceptibility.  Overall, our study does not support this hypothesis. Other than the trend for Primi NoD offspring to have less mastitis than Primi D offspring, we did not see an effect of dam disease history on offspring disease risk. Primi D offspring also had a lower risk of leaving the herd over their whole lifetime than Primi NoD offspring. The Carvalho study did not evaluate the effect of dam disease in primiparous dams. Variation in disease diagnosis and monitoring across years, personnel, and farms was hard to completely control for in the analysis. This is likely contributing to the variation between our results and previously published studies and makes studies like this difficult. 

Objective 2: Characterize associations between DNA methylation and disease status  

Description of Animal Characteristics and Global Methylation:

The pre-calving and 30 DIM samples from six case heifers with subclinical ketosis (TD) and their matched healthy controls (NoTD) were chosen for WGBS. Cases all with the same transition disease were chosen to minimize variation as it would be expected that different transition diseases would have different epigenetic effects. The disease status, calving date, disease date, and sampling dates are displayed in table 4 and case control pairs are listed together.

Across the 24 samples (12 animals at 2 time points), we identified 3,511,854 cytosine sites. The majority of them were CpG sites (76%) and 3.1% and 3.5% were CHG and CHH sites, respectively (Fig. 10A and D). Of the CpG methylation sites, most were either highly or lowly methylated and there were not any large differences in the distribution of genome methylation levels between TD and NoTD or between pre-partum and 30 DIM samples which is as expected (Fig. 10B-C, E-F). After smoothing across CpG methylation levels in 20kb windows, there were no large global differences in CpG methylation between TD and NoTD Heifers at either sampling time point (Fig. 10G-H). This is to be expected as all samples are from the same tissue type, blood, and it is more likely that TD would affect DNAm of a smaller number of specific sites.

Possible TD Risk or Protective DNAm loci:

Across the genome, 15,004 DMRs were identified (empirical p-value <0.05) between heifers that went on to have transition disease (TD) and healthy heifers (NoTD) in their pre-partum blood sample (Fig. 11A). Of the significant DMR, 64% were hypermethylated and 36% were hypomethylated in TD heifers versus NoTD heifers. DMR were significantly enriched in CpG islands, shores, shelves, and open sea contexts (Fig. 11B). They were also enriched in promoters, exons, introns, 3’ UTR, and intergenic functional regions. Hypomethylated DMRs were also enriched in 5’UTR (Fig. 11C).

Two machine learning algorithms were used to classify TD and NoTD heifers based on DMR identified pre-partum. There were 33 DMRs that were ranked important by both algorithms (Fig. 12). Only five of these DMR were hypomethylated and 28 were hypermethylated in TD versus NoTD heifers. DMR were annotated to mostly intergenic regions but 2 were annotated to promoters, 3 to exons, and 4 to introns. These DMR could represent potential risk/susceptibility loci for ketosis as all TD heifers had ketosis specifically. These regions will be characterized further in Obj. 4.

DNAm changes during the transition period:

Across the genome, 7,624 DMRs were identified (empirical p-value <0.05) between the pre-partum and 30 DIM samples of TD heifers (Fig. 13A). These DMR were significantly enriched in CpG islands, shores, shelves, and open sea contexts (Fig. 13B). They were also enriched in promoters, exons, introns, downstream, and intergenic functional regions. Hypermethylated DMRs were also enriched in 3’UTR (Fig. 13C).

Across the genome, 8,667 DMRs were identified (empirical p-value <0.05) between the pre-partum and 30 DIM samples of NoTD heifers (Fig. 14A). These DMR were significantly enriched in CpG islands, shores, shelves, and open sea contexts (Fig 14B). They were also enriched in promoters, exons, introns, 3’UTR, downstream, and intergenic functional regions (Fig 14C).

Two machine learning algorithms were used to classify pre-partum and 30DIM samples of TD heifers. There were 19 DMRs that were ranked important by both algorithms. Nine were hypomethylated and 10 were hypermethylated in 30 DIM samples compared to Pre-Partum samples (Fig. 15). These DMR were annotated to mostly intergenic regions but 2 were annotated to promoters, and 2 to introns.

Two machine learning algorithms were used to classify pre-partum and 30DIM samples of NoTD heifers. There were 24 DMRs that were ranked important by both algorithms (Fig. 16).  Half the DMR were hypomethylated and half were hypermethylated in 30 DIM samples versus pre-partum samples. These DMR were mostly annotated to intergenic regions but two were annotated to promoters, one to an exon, and six to introns.

DMRs identified between pre-partum and 30 DIM samples represent the changes in DNA methylation patterns animals experience in the transition period. DMR identified between these points in NoTD animals likely represent the changes healthy animals go through over this time and DMR within TD animals could represent the changes sick animals go through. There were 191 regions where both groups experienced DNAm changes (Fig. 17). The 7,433 DMR only seen in TD heifers could represent changes specifically related to their case of ketosis. The 8,476 DMR only seen in NoTD heifers could represent normal changes in DNAm that are prevented in animals with transition disease. These regions will be characterized further in Obj.4.

Objective 3: Characterize DNA methylation inheritance with a focus on dam disease status and the association with offspring disease susceptibility

Description of Animal Characteristics and Global Methylation:

All 12 calves sampled were sent for WGBS. The dam disease status, sex, birthdate, and sampling date are displayed in table 5.

Across the 12 samples, we measured on average 3,263,675 cytosine sites and 77% of them were in CpG contexts, 3.1% in CHG contexts and 3.6% in CHH contexts (Fig. 18A). Of CpG sites, most were either were highly or lowly methylated (Fig. 18B-C) which is as expected.  After smoothing across CpG methylation levels in 20kb windows across the genome, there were no large global differences in CpG methylation between Dam_TD and Dam_NoTD calves as shown by the PCA plot (Fig. 18D). We did not expect dam TD to have a large global effect on DNAm of her future offspring. It is more likely there are smaller changes.   

Offspring DNAm differences associated with Dam Disease Status:

To identify smaller differences in DNAm patterns between case and control calves, differential methylation analysis was conducted to compare 5 case calves whose dam experienced metritis in the lactation before their birth to their 5 matched control calves who had historically healthy dams. For this analysis, the calf whose dam had ketosis and his matched control were removed. This allowed us to have a more targeted analysis on understanding the effects of metritis. We expect different transition diseases could have different epigenetic programming effects. Overall, we identified 12,941 DMRs between calves whose dams did or did not have metritis (empirical p-value <0.05; Fig. 19A).  The DMRs were distributed across the genome and 41% were hypomethylated and 59% were hypermethylated in Dam TD calves. DMR were enriched in CpG Islands, shores, shelves, and open sea contexts (Fig. 19B). Hypermethylated DMR were enriched in promoters, exons, 3’UTR, and intergenic regions. Hypomethylated DMRs were enriched in promoters, introns, and intergenic regions (Fig. 19C). 

Machine learning algorithms were used to identify the most informative DMR to classify Dam TD and Dam NoTD calves. Nineteen DMR were commonly picked as informative by both machine learning algorithms. Of those, 4 were hypomethylated and 15 were hypermethylated in TD calves (Fig. 20). Most of these DMR were annotated to intergenic regions but one falls within a promoter and 2 within introns.

Objective 4: Identify differentially methylated regions with putative biological functions relevant to health, growth, production, fertility, and longevity relationships observed

Enrichment and Functional Analysis of TD Risk Loci:

DMRs were annotated to the nearest gene resulting in 5,763 differentially methylated genes (DMG) between TD and NoTD heifers pre-partum. Gene ontology (GO) analysis was run with all of these DMGs (Fig. 21). Enriched biological process GO terms included hormone metabolism and cell junction and adhesion. Enriched cellular component and molecular function GO terms also point to possible differences in cell signaling and transport between TD and NoTD heifers pre-partum.

Of all the DMG identified between TD and NoTD heifers pre-partum, 594 were genes in which the DMR fell within a promoter. On chromosome 10, one DMR that was identified as an important DMR through both machine learning algorithms falls within the promoter of the TSHR gene (Fig. 12). This gene codes for thyroid stimulating hormone receptor and thyroid activity differences in cattle with ketosis have been previously reported (Pehrson, et al., 1966). This promoter is hypermethylated in TD heifers suggesting a reduction in the expression of TSHR. Another important DMR from the machine learning classification falls within the promoter of the VCAN gene on chromosome 7 (Fig. 12). This gene codes for versican which is an extracellular matrix proteoglycan important in cell adhesion. In cattle, VCAN has been previously associated with oocyte number, embryo quality, and body measurements (Rocha, et al., 2024; Yan, et al. 2024).

Lastly, we wanted to compare our list of DMG’s to other genes previously implicated in ketosis. Soares, et al. previously conducted a meta-analysis of transcriptomic and proteomic studies about ketosis and negative energy balance in cattle (2021). They published a list of 2,740 differentially expressed genes (DEG) from transcriptomic studies. Yan, et al. conducted a multi-omics study related to ketosis in cattle and proposed a list of 5 candidate genes (2020). Our DMG were significantly overlapped with the Soares study genes. This means we identified changes in DNAm in and around genes that have been previously shown to be differentially expressed during ketosis (Fig. 22B). One DMG, GRINA, was also a candidate gene from the Yan study. This gene has also been previously associated with milk protein percent in cattle (Pedrosa, et al., 2021). This supports that the differential methylation seen between our TD and NoTD heifers pre-partum could be associated with their risk for ketosis and provides more support for the association of these genes with the disease.   

Enrichment and Functional Analysis of acquired DNAm changes with TD:

Of the DMR representing DNAm changes in the transition period only seen in NoTD heifers, 4,525 DMG were identified. Gene ontology analysis was run with all these DMGs (Fig. 23A). Enriched biological processes were related to the circulatory system. The DMG identified may be associated with changes in blood cell proportions during the transition period which has been previously demonstrated (Sun, et al., 2022).  Of the DMR representing DNAm changes in the transition period only seen in TD heifers, 4,184 DMG were annotated. Gene ontology analysis was run with all these DMGs (Fig. 23B). Enriched biological processes included epithelial cell proliferation and anatomical structure morphogenesis. Enriched molecular functions included peptidase activity and voltage-gated cation channel activity.

Of all the DMG identified between pre-partum and 30 DIM samples in TD heifers, 429 were genes in which the DMR fell within a promoter and 186 were genes in which the DMR fell within a gene body (exon or intron). NPSR1 was a DMG found between TD heifers pre-partum and at 30DIM and was not seen as differentially methylated in NoTD heifers over those timepoints. A DMR within an intron of the gene was also ranked with high importance for both machine learning algorithms (Fig. 15). This gene codes for a neuropeptide receptor and has been associated with fertility traits in cattle (Jiang, 2019).

Of all the DMG identified between pre-partum and 30 DIM samples in NoTD heifers, 502 were genes in which the DMR fell within a promoter and 256 were genes in which the DMR fell within a gene body (exon or intron). A DMR within the promoter of NMRAL1 was ranked with high importance for both machine learning algorithms (Fig. 16). NMRAL1 codes for a NADPH sensor protein. In NoTD heifers, the promoter of NMRAL1 was hypermethylated at 30DIM compared to Pre-Partum suggesting a decrease in expression in healthy cattle over this time.  NAD+/NADPH balance has been shown to be affected by ketogenic diets in humans and NADPH oxidase-ROS signaling pathway is activated in cows with subclinical ketosis (Elamin, et al., 2021; Song, et al., 2022).  It is possible the normal increase in DNAm of the NMRAL1 promoter is prevented by a ketosis event explaining why this was not identified as a DMR in TD heifers across the same timepoints.

Lastly, we wanted to compare the DMG that undergo DNAm changes in heifers that had ketosis to previously published gene lists associated with ketosis (Fig. 24A).  Our DMG significantly overlapped with the Soares study DEG’s and 1 gene, GRINA, was a ketosis candidate gene from the Yan study. The DMR near GRINA found between pre-partum and 30DIM samples of TD heifers was 32 kbp downstream of the transcription start site of the gene in an intergenic region. A DMR between TD and NoTD heifers pre-partum was also annotated to this gene but it was about 20 kbp upstream of the transcription start site. This provides even more evidence supporting the association of GRINA with ketosis in dairy cattle.   

Enrichment and Functional Analysis of inherited DNAm differences due to Dam TD:

DMRs were annotated to the nearest gene resulting in 5,130 DMGs between Dam TD and Dam NoTD calves. Gene ontology analysis was run with all of these DMGs (Fig. 25).  There are biological processes, cellular components, and molecular function GO terms related to cellular membrane transport and signaling.

Of all the DMGs, 1,036 were genes in which the DMR fell within the gene’s promoter. On chromosome 13, one DMR that was ranked as important by both machine learning algorithms falls within the promoter of the CEBPB gene (Fig. 20). The C/EBPβ transcription factor this gene codes for has been previously implicated in epigenetic memory of immune response in hematopoietic stem cells (De Laval, et al., 2020). This gene’s promoter was hypermethylated in TD calves suggesting decreased expression of the gene that could affect the future immune response of calves born to dams with a history of metritis. 

Lastly, we wanted to compare our list of DMG’s to other genes previously implicated in inheritance of maternal transition disease and metritis more broadly. Zhang, et al. previously identified 1,364 DMG’s in calves at birth associated with maternal metabolic stress pre-conception. Metabolic stress is a common cause and outcome of transition disease (2022). Brewer, et al. previously published a list of 2,211 DEG in endometrial cells of cows that had uterine infection compared to healthy cows (2020). Our DMG were significantly overlapped with both gene lists (Fig. 26B). The overlap with the Brewer gene list suggests that calves born to dams with a history of metritis may be inheriting some of the DNAm changes associated with gene expression changes seen in cows who have metritis themselves. The overlap with the Zhang study genes suggests there may be some DNAm changes associated with dam transition disease more broadly that are not disease specific.  

Discussion and Limitations:

Overall, we were able to propose three types of DNAm loci associated with transition disease in dairy cattle: potential risk loci, loci that undergo changes in DNAm around the time of disease, and loci in which dam disease status may impact offspring DNAm. There was overlap in the genes that potential risk DMR and DMR undergoing changes due to transition disease were annotated to, though the actual DMR regions were unique. GRINA is an example of one of these genes and is a promising candidate gene for ketosis in cattle. Differential methylation in or around this gene demonstrated in our study may help explain some of the functional connection between this gene and the ketosis phenotype.  

There was no overlap in the DMG identified in calves and the DMG identified in heifers in our study. This is likely due to fact that the actual transition disease captured in these populations were different. The dams of our Dam TD Calves had metritis while the TD heifers had ketosis. Nonetheless, the DMR identified between Dam TD and NoTD calves within the promoter of the CEBPB gene is very promising as a potential locus for epigenetic programming of dam disease history. The transcription factor it codes for is known to play a role in immune response and it has been implicated in the epigenetic memory of immune response in hematopoietic stem cells.

While these results are very promising, three limitations of our study should be recognized. Firstly, we only looked at DNAm within one tissue, blood, and did not account for the heterogeneity of cell types within the samples. While the majority of DNA isolated from a blood sample will be from white-blood cells, there are cell-specific DNAm patterns for different white blood cells. We did not account for this and cannot definitively say if the differential methylation we observed is simply from differing blood cell proportions in samples. Secondly, we did not intensively diagnose disease and instead relied on farm diagnosis to classify our case and control animals. This could lead to some misclassification of disease and a loss of power. Since we were able to find numerous significant DMRs, this is likely not impacting us too greatly but should be recognized.  Thirdly, due to limitation in the number of dams with transition disease we could find, we did not sample enough calves to follow their specific performance outcomes. We can only apply what we learned from the whole-herd retrospective records analysis and the known functions of DMGs to predict how the DMR identified between TD and NoTD calves may affect future performance.