Comparison of the transcriptome in circulating leukocytes in early lactation between primiparous and multiparous cows provides evidence for age-related changes.

BACKGROUND: Previous studies have identified many immune pathways which are consistently altered in humans and model organisms as they age. Dairy cows are often culled at quite young ages due to an inability to cope adequately with metabolic and infectious diseases, resulting in reduced milk production and infertility. Improved longevity is therefore a desirable trait which would benefit both farmers and their cows. This study analysed the transcriptome derived from RNA-seq data of leukocytes obtained from Holstein cows in early lactation with respect to lactation number.
RESULTS: Samples were divided into three lactation groups for analysis: i) primiparous (PP, n = 53), ii) multiparous in lactations 2-3 (MP 2-3, n = 121), and iii) MP in lactations 4-7 (MP > 3, n = 55). Leukocyte expression was compared between PP vs MP > 3 cows with MP 2-3 as background using DESeq2 followed by weighted gene co-expression network analysis (WGCNA). Seven modules were significantly correlated (r ≥ 0.25) to the trait lactation number. Genes from the modules which were more highly expressed in either the PP or MP > 3 cows were pooled, and the gene lists subjected to David functional annotation cluster analysis. The top three clusters from modules more highly expressed in the PP cows all involved regulation of gene transcription, particularly zinc fingers. Another cluster included genes encoding enzymes in the mitochondrial beta-oxidation pathway. Top clusters up-regulated in MP > 3 cows included the terms Glycolysis/Gluconeogenesis, C-type lectin, and Immunity. Differentially expressed candidate genes for ageing previously identified in the human blood transcriptome up-regulated in PP cows were mainly associated with T-cell function (CCR7, CD27, IL7R, CAMK4, CD28), mitochondrial ribosomal proteins (MRPS27, MRPS9, MRPS31), and DNA replication and repair (WRN). Those up-regulated in MP > 3 cows encoded immune defence proteins (LYZ, CTSZ, SREBF1, GRN, ANXA5, ADARB1).
CONCLUSIONS: Genes and pathways associated with lactation number in cows were identified for the first time to date, and we found that many were comparable to those known to be associated with ageing in humans and model organisms. We also detected changes in energy utilization and immune responses in leukocytes from older cows.

Background

Longevity is an economically important trait in dairy cows, which also has welfare implications. Cattle can potentially live for over 20 years, but this is rare in practice and the average lifespan in dairy cows is currently around 4.5 to 7 years [1-3]. To maximise economic potential, heifers should first calve at 24 months of age so becoming primiparous cows, starting their first lactation, and beginning to pay back their rearing costs through the production of saleable milk [4]. Ideally, they should continue to calve at annual intervals, with their milk production potential increasing progressively until the fourth or fifth lactation [5]. Cows which survive for longer achieve greater lifetime milk production associated with higher profitability [2, 3, 6] together with the benefit of reduced greenhouse gas emissions [7]. Many cows do not achieve optimum longevity due to premature involuntary culling, for which the main reasons are mastitis, infertility and lameness [8, 9]. On the other hand, voluntary culling may be used to remove cows with low milk yields or behavioural issues and to increase the rate of genetic gain [10].

There are major alterations in the metabolic profiles associated with the start of lactation [11]. These also change with age as primiparous cows are still growing, so their energetic demands for growth compete with those of milk production, whereas in older cows the higher milk yields are associated with greater mobilisation of body tissue, leading to higher circulating concentrations of nonesterified fatty acids (NEFA) and beta-hydroxybutyrate (BHB) [12, 13] Circulating concentrations of IGF1 are also significantly lower in multiparous (MP) than primiparous (PP) cows in early lactation [14]. IGF1 is a key metabolic hormone which provides a good indication of the energy balance status of the animal and is also linked to the immune responses following calving [15]. There are therefore important metabolic, endocrine and physiological changes that occur as cows mature, stop growing and increase their milk yields. Cows with greater milk production potential face an increased risk of glucose shortages in their immune cells, which contributes to immune dysfunction in the peripartum period [16]. Both metabolic and immune dysfunction therefore impact on the transcriptomic changes in leukocytes in early lactation.

Recent studies in human populations have highlighted age-associated changes in leucocyte functionality affecting innate and adaptive immune functions [17]. The causes and consequences of ageing on the human blood transcriptome have, however, proved difficult to dissect due to interactions with environmental influences, genetic factors and a large number of age-related diseases [18]. Studies on model organisms have highlighted that ageing is characterized by many alterations at molecular, cellular and tissue level [19]. Studies of the ageing transcriptome have been performed in species including C. elegans [20], flies [21], rodents [22] and humans [23]. This approach has identified various signatures found to occur repeatedly across different tissues and organisms. Candidate genes whose expression is consistently associated with cellular ageing have been classified into six categories: i) downregulation of genes encoding mitochondrial proteins, ii) downregulation of the protein synthesis machinery, iii) dysregulation of immune system genes, iv) a reduction in growth factor signalling, v) constitutive responses to stress and DNA damage, and vi) dysregulation of gene expression and mRNA processing [24].

Although all living creatures age at some point, our knowledge on the biology of ageing is still not sufficient. The physiological process of ageing in humans is associated with a progressive loss of function and increased vulnerability to disease, frailty, and disability [25]. Because the incidence of adult diseases increases with age, a better understanding of the biology of ageing could greatly improve our efforts to elucidate the physiopathology of such conditions [26].

In this study we have compared the leucocyte transcriptome between young (primiparous, PP) and older (multiparous, MP) cows in order to shed light on the genes and related genomic pathways involved in age-related symptoms arising during the different phases of a cow’s life. This has the potential to inform both breeding and management practices, so providing a significant gain to both agricultural production and animal welfare.

Results

Weighted gene co-expression network analysis (WGCNA) to determine the relationships of leucocyte gene transcription between cow parity

Whole blood transcriptomes of 229 cow samples from six dairy herds were obtained in early lactation at 14 ± 0.1 days in milk. The PP cows, as expected, produced less milk, but there was no difference in milk production between the MP2–3 and MP > 3 cows, with milk yields at the time of sample collection of PP 23.5 ± 0.88a (n = 53); MP2–3 33.8 ± 0.68b (n = 121) and MP > 3 33.0 ± 1.14b (n = 55) kg/d respectively (mean ± SEM, b > a, P < 0.001; Supplementary Table 1). There was also a minor difference in milk production between two of the herds, with ITA cows producing slightly less and DEU cows slightly more milk (Supplementary Table 1). The herd effect was, however, accounted for in the WGCNA analysis by removing the batch effect (Supplementary Fig. 1).

The leucocyte RNA samples had an average mapping rate against the reference genome of 96.2%, resulting in an average total number of reads of 34,439,525 (Supplementary Table 2). As we were most interested in understanding the effect of ageing, we contrasted the leukocyte expression between the two groups with extreme lactation numbers PP vs MP > 3, using those classified as MP2–3 as background. Firstly, we performed a differential expression analysis using DESeq2. Over 5000 genes out of the 17,216 genes initially mapped were found to be differentially expressed between lactation groups following Benjamini-Hochberg false discovery rate adjustment (padj < 0.1). This list was narrowed down to 2925 differentially expressed genes (DEG) when contrasting PP vs MP > 3 cows (Supplementary Fig. 2). We then applied WGCNA to identify gene modules from our dataset associated to the trait of lactation number. From an initial total of 21,207 genes, 13,769 (64.9%) genes passed the DESeq2 filtering steps. Samples from three cows were removed at this stage as outliers (two from MP > 3 and one from MP2–3 age group, respectively; Supplementary Fig. 3), leaving a total of 226 samples. This unsupervised technique identified 32 interconnected gene modules from the filtered gene list (Fig. 1A). The number of genes per module is shown in Fig. 1B. Of these, the violet and turquoise modules were significantly related to the trait lactation number in both the PP and MP > 3 groups (Fig. 1C). As expected, modules with a significant positive correlation in the younger cows were negatively correlated with lactation number in the older ones, and vice versa, while no modules showed a significant correlation with the medium lactation number group (MP2–3). Within these seven significant modules, 2274 genes were more highly expressed in the leukocytes of the PP cows and 2721 genes were more highly expressed in the leukocytes of the MP > 3 cows. The top 20 genes with the highest differential expression between PP and MP > 3 cows are listed in Table 1.

Fig. 1

A Cluster dendrogram of 13,769 genes (individual black lines at top) clustered by their topological overlap dissimilarity scores. The multi-coloured panel next to “Dynamic Tree Cut” shows 42 identified modules using the Dynamic Tree Cut algorithm. The second multi-coloured panel shows 32 larger modules identified after highly correlated smaller modules were merged together (agreement of 0.8). B Number of genes per module identified using WGCNA. C Module-trait correlations according to the parity of the cow: i) primiparous (PP, n = 53), ii) multiparous in lactations 2–3 (MP 2–3, n = 121), and iii) MP in lactations 4–7 (MP > 3, n = 55). Red, positive correlation, blue, negative correlation. Numbers indicate the r value: significance was taken at P = ± 0.25

Table 1

Top 20 genes with the highest differential expression between PP and MP > 3 cows. Positive fold change (FC), more highly expressed in PP cows; negative FC (shaded), more highly expressed in MP > 3 cows

Gene Symbol	Gene name	logFC	adj.p.value
NELL2	neural EGFL like	−1.581	1.40E-63
ZNF462	zinc finger protein 462	1.931	1.79E-38
PLCL1	phospholipase C like 1 (inactive)	1.042	2.45E-37
IGF2BP3	insulin like growth factor 2 mRNA binding protein 3	2.268	6.36E-37
CNTNAP1	contactin associated protein 1	1.675	2.85E-35
SOX13	SRY-box transcription factor 13	1.925	8.17E-35
BLK	BLK proto-oncogene, Src family tyrosine kinase	0.881	1.23E-32
MYBL1	MYB proto-oncogene like 1	−0.771	1.93E-32
CD163L1	CD163 molecule like 1	1.463	6.55E-32
WNT5A	Wnt family member 5A	2.673	2.56E-31
SCRN1	secernin 1	1.129	5.03E-31
LOC515828	uncharacterised	1.207	6.12E-30
CSPG4	chondroitin sulfate proteoglycan 4	2.831	8.30E-30
DBN1	drebrin 1	0.942	1.95E-29
PTPRK	protein tyrosine phosphatase receptor type K	−3.301	2.00E-29
LAMA4	laminin subunit alpha 4	2.998	4.15E-29
TCAF1	TRPM8 channel associated factor 1	1.050	1.17E-28
LOC751811	WC1 isolate CH211	1.389	2.85E-28
SOX4	SRY-box transcription factor 4	0.820	3.69E-28
LOC104969122	uncharacterised	−1.396	9.97E-28

Significant modules in leukocytes of both PP and MP > 3 cows

The violet module containing 87 genes had the highest correlations with lactation number (+ 0.59 in PP and − 0.43 in MP > 3 cows, Fig. 1C). Within this module the LAMA4 gene was the most highly correlated gene (Table 1, gene significance = 9.05 × 10− 22), with lower expression in the MP > 3 cows, representing its module membership. This is a laminin gene, part of the family of extracellular matrix glycoproteins which form the major non-collagenous constituent of basement membranes. GO enrichment analysis of the genes in this module revealed the terms receptor binding and cell periphery (Fig. 2A).

Fig. 2

GO enrichment analysis of significantly correlated modules (FDR < =0.1). A Positively correlated in PP cows: violet, darkred, midnightblue, orange, yellow; B negatively correlated in PP cows: black, turquoise

The turquoise module represents the second module which had significant correlations in opposite directions with PP and MP > 3 cows, of − 0.29 and 0.33, respectively. This module contained 1818 genes of which NELL2 was the most highly correlated with lactation number (Table 1). This was also overall the most differentially expressed gene in leucocytes with higher expression in MP > 3 cows. It was a large module containing 1818 genes with a wide variety of functions (Fig. 2B, Supplementary Table 3). The top six most significant GO terms (all < 1.36 × 10− 25 and containing between 389 and 959 genes) were extracellular region, membrane, single-organism process, vesicle and cell periphery.

Modules up-regulated in leukocytes of PP cows

The modules midnightblue, darkred, and orange were all significantly positive in PP cows with correlations spanning from 0.27 to 0.25, implying lower expression in leukocytes as lactation number increased (Fig. 2A, Supplementary Table 3). The most significant GO terms in the midnightblue module were all biological processes: protein autophosphorylation, immune system process, cell activation and regulation of cell-cell adhesion. The darkred module had enrichment of biological regulation and the molecular function term cell periphery, while the orange module was mainly enriched with genes involved in cilium and cilium organisation; RNA processing and regulation of gene expression, and microtubule organizing center.

Black module: down-regulated in leukocytes of PP cows

The black module was down-regulated in leukocytes of PP (correlation of − 0.25) and almost significantly up-regulated in MP > 3 (correlation of 0.24). Genes were therefore more highly expressed in the older cows with higher lactation numbers. GO enrichment analysis of genes in this module (Fig. 2B, Supplementary Table 3) showed that the most significant terms all related to cellular components (membrane, vesicle, cytoplasm, extracellular region and cell). Within the biological process category single-organism process, monocarboxylic acid metabolic process, small molecule catabolic process and negative regulation of cellular process were significant.

Yellow module: down-regulated in leukocytes of MP > 3 cows

Lastly, the yellow module was significantly down-regulated in leukocytes as lactation number increased (MP > 3; correlation − 0.27). In this module IGF2BP3 (insulin like growth factor 2 mRNA binding protein 3) was the most highly correlated gene (Table 1). Within this module the two cellular component terms binding and nucleic acid binding transcription factor activity were most highly represented, together with three biological processes (regulation of gene expression, regulation of nitrogen compound metabolic process and regulation of cellular macromolecule biosynthetic process) (Fig. 2A, Supplementary Table 3).

Functional annotation cluster analysis

The next stage of the analysis involved pooling all the genes from the modules which were more highly expressed in the PP cows (yellow, orange, midnightblue, darkred, violet) and those which were more highly expressed in the MP > 3 cows (black and turquoise). These two gene lists were then subjected to David functional annotation cluster analysis [27]. All the results are provided in Supplementary Table 4. The top six most significant clusters in each category are summarised in Table 2, with the genes in each of these clusters listed in Supplementary Tables 5A and B.

Table 2

Summary of the six most significant clusters obtained from DAVID functional annotation cluster analysis of the DEG in the modules which were significantly positively or negatively correlated to lactation number in primiparous (PP) cows, indicating either lower or higher gene expression in older multiparous MP > 3 cows

Positive Modules	Negative Modules
Cluster 1: Enrichment Score: 16.18	Cluster 1:Enrichment Score: 29.73
• Krueppel-associated box • Zinc finger C2H2-type/integrase DNA-binding domain • Zinc finger, C2H2 • Nucleic acid binding • Zinc finger, C2H2-like • KRAB • ZnF_C2H2 • Regulation of transcription, DNA-templated • Metal ion binding	• Membrane • Transmembrane helix • Transmembrane • integral component of membrane
Cluster 2: Enrichment Score: 7.01	Cluster 2: Enrichment Score: 19.87
• Zinc-finger • Metal ion binding • Zinc • Metal-binding	• Signal • Disulfide bond • Glycoprotein • Signal peptide • Glycosylation site: N-linked (GlcNAc...) • Disulfide bond
Cluster 3: Enrichment Score: 6.83	Cluster 3: Enrichment Score: 7.32
• Transcription regulator SCAN • Retrovirus capsid, C-terminal • SCAN	• Transmembrane region • Topological domain: Cytoplasmic • Topological domain: Extracellular
Cluster 4: Enrichment Score: 5.05	Cluster 4: Enrichment Score: 3.94
• Nucleotide binding • Nucleotide-binding, alpha-beta plait • RNA recognition motif domain • RRM	• Glycolysis / Gluconeogenesis • Glycolysis • Glycolytic process • Carbon metabolism • Biosynthesis of amino acids
Cluster 5: Enrichment Score: 2.02	Cluster 5: Enrichment Score: 3.70
• BTB/POZ fold • BTB/POZ-like • BTB	• C-type lectin fold • C-type lectin-like • C-type lectin • CLECT • C-type lectin, conserved site
Cluster 6: Enrichment Score: 1.88	Cluster 6: Enrichment Score: 3.28
• Immunoglobulin-like domain • Immunoglobulin-like fold • Immunoglobulin subtype 2 • Immunoglobulin I-set • Immunoglobulin subtype • IGc2 • IG	• Immunity • Innate immunity • Innate immune response

Clusters 1, 2 and 3 from the modules which were positively correlated in the PP cows were all involved in the regulation of gene transcription. Cluster 1 consisted almost entirely of 56 genes encoding zinc finger proteins. The most significant term in Cluster 1 was Krueppel-associated box. All but five of the Cluster 1 genes were included in the much larger Cluster 2 with 142 genes, for which the most significant terms were Zinc-finger and metal ion binding. Cluster 2 also contained all 18 genes from Cluster 3, for which the three significant terms were Transcription regulator SCAN, Retrovirus capsid, C-terminal and SCAN. Cluster 4 contained 61 genes involved in nucleotide binding. Cluster 5 contained 31 genes for which the most significant term was BTB/POZ fold. Finally, Cluster 6 contained 50 genes with the most significant term being Immunoglobulin-like domain. Overall, this analysis shows that the most significant functions which were more highly expressed in the first lactation cows all related to the regulation of gene transcription.

Turning to the genes which were more highly expressed in the MP > 3 cows, Clusters 1, 2 and 3 were all large clusters containing 471, 435, and 254 genes, respectively, with quite general functions relating to the cell surface, transport across it and signalling. In Cluster 1 the top three terms were Membrane, Transmembrane helix and Transmembrane, in Cluster 2 they were Signal, Disulfide bond and Glycoprotein and in Cluster 3 transmembrane region, topological domain: Cytoplasmic and topological domain: Extracellular. The top term in Cluster 4 was Glycolysis/Gluconeogenesis and this contained 27 genes, nearly all of which encoded enzymes involved in the glycolytic pathway (Fig. 3).

Fig. 3

David functional annotation Cluster 4 contained many genes (shown in blue) in the Glycolysis/Gluconeogenesis KEGG pathway which were more highly expressed in MP > 3 vs PP cows (figure drawn by the authors). This would potentially maximise the use of glucose as a precursor of acetyl CoA for use in the TCA cycle

The top term in Cluster 5 was C-type lectin fold, and 11 of the 25 genes in this cluster encoded C-lectins. Finally, Cluster 6 contained 45 genes involved in Immunity and Innate immunity. These are illustrated in Fig. 4, which shows that they mainly encoded proteins for pattern recognition receptors (PRR), the bovine MHC complex, the complement system and several with antimicrobial activity. Another noteworthy cluster which was negatively related to lactation number was Cluster 16 (Supplementary Table 4). The top pathway in this was Fatty acid degradation (P < 0.0013). This cluster contained a number of genes encoding enzymes in the mitochondrial beta-oxidation pathway, which converts fatty acids into acetyl CoA to enter the TCA cycle (ACSL1, ACSL4, ACSL5, ACSL6, CPT1, CPT2, ACOX1, ACADVL, ECHS1, ACAD5, ACAA1, ACAT2, ACADS, ACADVL). This implies greater use of this pathway to supply energy in the leukocytes of the PP cows.

Fig. 4

David functional annotation Cluster 6 contained genes involved in immunity which were more highly expressed in MP > 3 vs PP cows. This included genes encoding proteins for pattern recognition receptors (PRR), the bovine MHC complex, the complement system and several with antimicrobial activity. NRROS and TNFAIP8L2 can both act as negative regulators of TLR signalling mediated via inhibition of NF-kappa-B activation

Comparison with other species

Finally, we investigated the overlap between genes significantly associated with lactation number in leucocytes in our study and 170 candidate genes identified from previous studies of ageing in humans and model organisms (mainly mice and C. elegans [18]). Almost one quarter of these candidate genes (38/170, 22.3%) were present in our list of DEG, of which 14 and 24 respectively were from positive and negative significantly correlated modules, indicating reduced or increased expression in older cows. Of these 16 of the 38 DEG were previously identified as being age-related in the human blood transcriptome ([18]; Table 3).

Table 3

List of genes overlapping candidate age-related genes found in humans and with differential expression between PP (highly expressed) and MP > 3 (highly expressed) cows

	Gene name	Module	Description	Reference
PP	CCR7	darkred	Control the migration of memory T cells to inflamed tissues, as well as stimulating dendritic cell maturation	[28]
	CD27	midnightblue	Required for generation and long-term maintenance of T cell immunity	[28]
	FOXO1	midnightblue	Transcription factor which may play a role in myogenic growth and differentiation	[29]
	IL7R	midnightblue	Plays a critical role during lymphocyte development (age and longevity expression association)	[30]
	CAMK4	midnightblue	Implicated in transcriptional regulation in lymphocytes, neurons and male germ cells (longevity genetics candidate)	[31]
	CD28	midnightblue	Essential for T-cell proliferation and survival, cytokine production, and T-helper type-2 development	[28]
	WRN	yellow	RecQ helicase important in the maintenance of genome stability. Mutations cause Werner syndrome (premature ageing)	[32]
	MRPS27	yellow	Mitochondrial ribosomal protein assisting protein synthesis within the mitochondrion (mitochondrial ageing)	[33]
	MRPS9	yellow	Mitochondrial ribosomal protein assisting protein synthesis within the mitochondrion (mitochondrial ageing)	[33]
	MRPS31	orange	Mitochondrial ribosomal protein assisting protein synthesis within the mitochondrion (mitochondrial ageing)	[33]
MP > 3	SREBF1	black	Transcription factor involved in sterol biosynthesis	[29]
	ANXA5	turquoise	Potential role in cellular signal transduction, inflammation, growth and differentiation	[34]
	LYZ	turquoise	Antimicrobial peptide	[34]
	CTSZ	turquoise	Lysosomal cysteine protease with role in immune defense	[34]
	ADARB1	turquoise	Enzyme in the A (adenosine) to I (inosine) RNA editing pathway, involved in a general maintenance of cellular health (longevity genetics candidate)	[35]
	GRN	turquoise	Granulin precursor protein which regulates cell growth. Also, co-ordinates CpG trafficking to TLR9, promoting activation of signalling cascade	[34]

Those up-regulated in PP cows were mainly associated with T-cell development and function (CCR7, CD27, IL7R, CAMK4, CD28) and protein synthesis within the mitochondria (MRPS27, MRPS9, MRPS31), and also included WRN, a gene associated with premature ageing in humans [36]. Those up-regulated in MP > 3 cows included genes encoding proteins involved in immune defence (LYZ, CTSZ, SREBF1, GRN, ANXA5 and ADARB1).

Discussion

Livestock species are raised primarily for their economic benefit to humans. Most dairy cows are culled before they reach the end of their potential lifespan due to poor milk production or fertility and/or an increased prevalence of diseases such as those causing mastitis or lameness [3, 8, 9]. Various cow longevity indexes have been defined, some of which also take account of lifetime milk production, which is in turn affected by both the milk yield capacity and the number of lactations achieved [3]. The main focus of previous studies into longevity in cows has been to increase the average survival time in the milking herd, so improving the profitability of the dairy industry. For example, one genome wide association study (GWAS) which investigated longevity found genes such as NPFFR2, previously identified as a candidate gene for mastitis resistance and two zinc finger proteins (ZNF613, ZNF717), which have been associated with calving difficulties [37]. Information about the age-related morbidities and causes of death that afflict cattle due to natural ageing is, however, limited. In contrast, there is a growing body of previous work into the underlying causes of cellular ageing which has been based on studies of human populations and model organisms.

This study is the first, to our knowledge, to assess changes in the global gene expression in leukocytes of dairy cows associated with increasing lactation number. For this we compared first lactation PP cows versus older multiparous MP > 3 cows using WCGNA analysis in order to identify potential genes and related pathways involved in the ageing process. The blood samples were all collected in early lactation, a time when lactating cows are placed under metabolic stress due to the nutrient requirements of milk production [11, 38] but none are pregnant. These aspects improved comparability between samples but meant that we did not include any data from younger growing heifers of < 2 years. The general aim of most dairy farmers is to calve animals for the first time at 2 years of age and then to achieve an annual calving interval, although in practice these targets are rarely met [39]. Although we did not record the exact ages of the cows in this study, the expectation would be for the PP cows to have been between 2 and 2.5 years old, the MP2–3 between 3.5–6 years old and the MP > 3 (which were in lactations 4–7) would have been between 6 and 10 years old. In most countries worldwide, dairy cow longevity has declined over the past 50 years and has been negatively related to the rise in milk yields achieved over the same period [3]. Despite the many differences in physiology, lifestyle and lifespan between species, many of the genes and pathways identified in our study as being associated with increased lactation number in dairy cows were nevertheless the same as those highlighted in previous studies of ageing in other species [18–20, 23].

All the samples analysed were whole blood leukocytes. In dairy cows the average leucocyte population is between 5 and 12 cells per ml of blood, predominantly consisting of lymphocytes and neutrophils [40]. The blood does, however, also contain natural killer cells, platelets, peripheral blood mononuclear cells (PBMC), eosinophils, and basophils. These different cell types all vary in both their basal gene expression and their transcriptional amplification responses to particular stimuli [41]. Furthermore, many aspects of the immune system alter during ageing, eventually leading to immunosenescence [24]. During this process different immune cell subsets are affected in different ways [17]. There is also an increase in baseline systemic inflammation with age, termed “inflammaging” [42, 43]. This may arise as a result of multiple mechanisms, including the accumulation of misfolded proteins, impaired clearance of senescent cells and obesity. Most individual leukocytes have a short lifespan in the circulation, measured in days or weeks, before they are destroyed by the lymphatic system, although there is a small pool of long-lived T and B lymphocytes which can survive for years, providing immunological memory. In humans, the relative abundance of naive T-cells decreases with chronological age while the population of memory T-cells increases [18, 44, 45]. Another possible source of variation in our analysis was the potential differences between herds, but this was accounted for by removing the batch effect in the WGCNA. Two of the six herds also fed three different lactational diets. Although these did influence cow metabolism this was not yet apparent at Day 14 when the samples used in the present study were collected. A comparison of the transcriptome in PBMC collected from the same UK herd failed to detect any DEG attributable to diet at this time [46].

The leukocyte transcriptome thus provides a reflection of the basic functions required for cell survival together with the various responses of the different cell types to the changing environment within the body, which alters with key factors such as disease exposure and nutrition. Cells may proliferate, undergo apoptosis or migrate to or from tissues into blood in response to different signals [47]. The changes reported here which were associated with lactation number will, to some extent, reflect altered abundance of different cell types as well as their changing expression patterns. Despite this caveat, many previous studies have now reported that transcriptional signatures of whole blood can reliably differentiate individuals with a variety of infections, for example Johne’s disease (Mycobacterium avium subsp. Paratuberculosis) in cattle [48]. Our analysis also benefitted from the use of WGCNA, which identifies networks of co-expressed genes whose expression is highly correlated. The use of whole blood also avoided the pitfall of potential artefacts which can be induced during cell separation procedures and provided greater consistency in a multi-site study.

Previous work based on studies of human populations and model organisms has identified a number of transcriptional signatures for cellular ageing which occur repeatedly across different tissues and organisms and which segregate into six main groups [24, 49]. In brief these are: i) downregulation of genes encoding mitochondrial proteins; ii) downregulation of the protein synthesis machinery (including ribosome biogenesis); iii) dysregulation of immune system genes, immune senescence; iv) reduction in growth factor signalling; v) constitutive responses to stress and accumulated DNA damage, and vi) dysregulation of processes regulating gene expression and mRNA processing (transcription and translation). We have obtained evidence that all of these were associated with lactation number within our population of cows.

Mitochondria and oxidative stress

Mitochondria regulate a multitude of different metabolic and signalling pathways and also play an important role in programmed cell death [50]. Oxidative metabolism causes endogenous production of free radical molecules and oxidative damage accumulates in multiple tissues and species with age [49]. For example, accumulated mutations in somatic mitochondrial DNA (mtDNA) and respiratory chain dysfunction were associated with ageing in mice [51]. In our study, genes encoding proteins involved in fatty acid beta-oxidation were more highly expressed in leukocytes in PP cows. There were also changes in expression of genes encoding mitochondrial ribosomal proteins, with higher expression of MRPS9, MRPS27 and MRPS31 in leucocytes of PP cows whereas MRPL17 and MRPL38 showed enriched expression in MP > 3 cows.

Protein synthesis

Protein homeostasis is essential to maintain protein structure and function, but the control of this process declines during ageing [49]. In our study, GO enrichment analysis of all the genes in the violet module, which were more highly expressed in the first lactation cows, revealed their involvement in protein autophosphorylation and catalytic activity (enzyme activity). Autophosphorylation is a type of post-translational modification of proteins. In eukaryotes, this process occurs by the addition of a phosphate group to serine, threonine or tyrosine residues within protein kinases, normally to regulate the catalytic activity. Genes involved in regulation of nitrogen compound metabolic process were identified in the yellow module as being down-regulated in the older MP > 3 cows. Maintenance of the proteome is essential to enable cells to respond appropriately to their environment. This requires correct synthesis and assembly of proteins and is controlled by molecular chaperones and clearance mechanisms which help to prevent protein misfolding and the associated accumulation of toxic aggregates. The efficiency of this process declines with age and has previously been associated with both metabolic and immunological diseases [49, 52].

Immune system

As we were studying a leukocyte population, changes in gene expression relating to immune function were expected. As animals age they are exposed to an increasing variety of disease causing microorganisms, while a progressive loss of function of the immune system increases their vulnerability to infection [25]. Notable age-related changes within the immune cell population include reduced cytokine signalling, diminished production of nitric oxide and peroxide, decreased phagocytic ability and reduced ability of dendritic cells to migrate and process antigens [17]. Our results are in accord with the study by Peters et al. [18], who investigated leucocytes from ageing human populations, in finding pathways which were either up- or down-regulated with increasing age. The GO term immune system process was enriched in both PP and MP > 3 cows, with genes involved in adaptive immunity up-regulated in the PP cows (FYN, ITK, LCK, etc) and genes related to innate immunity up-regulated in MP > 3 cows (CTSS, CTSH, IRAK2, TLR2, etc). The black module was down-regulated in the PP cows and contained genes involved in the regulation of cytokine secretion (IL10, CD14, FGR, IL17RC). Expression of genes in the turquoise module increased in older cows, containing genes involved in neutrophil degranulation and innate immune system. The darkred and midnightblue modules were both more highly expressed in the PP cows and contained the terms disease and immune system.

There was also a significant overlap between the DEG with immune functions identified in our population of cows and candidate age-related genes in the human transcriptome [18]. Those up-regulated in PP cows were mainly associated with T-cell development and function (CCR7, CD27, IL7R, CAMK4, CD28) while those up-regulated in MP > 3 cows included genes encoding proteins involved in immune defence. Of these LYZ encodes lysozyme, an antimicrobial peptide, CTSZ encodes cathepsin Z a lysosomal cysteine protease with multiple roles in host immune defense mechanisms, SREBF1 encodes a transcription factor involved in TLR4 signalling, while GRN encodes granulin, involved in TLR9 signalling, ANXA5 is involved in T-cell activation and ADARB1 encodes a deaminase enzyme with A-to-I RNA editing activity, which is important for the maintenance of cellular health but may also play a role in response to viral infection. The results from our study therefore suggest that the higher lactation number cows are more actively engaged in combatting disease pathogens through activation of the innate immune system and also support a higher level of inflammation with ageing. Parturition is itself an inflammatory process and, in addition to higher rates of infectious diseases including mastitis and metritis, older cows are also more vulnerable to metabolic disorders including milk fever, ketosis and displaced abomasum after calving. The prevalence of all of these metabolic diseases increased significantly in cows which were ≥ 5 lactations [53]. Our previous study showed that metabolic disorders led to prolonged uterine inflammation by up-regulating the genes and pathways related to immune and inflammatory processes [15].

Growth factor signalling

The relationship between ageing and metabolic regulation is bidirectional, as ageing impairs the activity of key metabolic signalling pathways and the ensuing metabolic dysregulation results in accelerated ageing [19]. Cell signalling pathways that sense the availability of nutrients and the energy status of the cells communicate with hormonal and growth factor signalling pathways to co-ordinately regulate whole body metabolic homeostasis. Ageing results in a gradual deterioration of various cellular functions including metabolic regulation [54]. The turquoise module was the second highest correlated module with lactation number, containing genes more highly expressed in MP > 3 cows. The insulin-like growth factor binding protein pathway was indeed enriched with genes like INSR, IGF1R, IGF1, LDLR, HTRA1, IGFBP7 etc. (Cluster 8; Supplementary Table 4). The latter pathway is interlinked with metabolic pathways to ensure coordinate regulation and fine-tuning of cellular metabolic responses in line with cellular energy status, nutrient availability and hormonal/growth factor signalling input [54]. We have shown previously that circulating IGF1 concentrations are significantly lower in PP compared with MP cows, falling to a lower nadir in the first week after calving [14].

Stress and DNA damage

Accumulation of genetic damage represents one of the major contributions to ageing of cells and organisms. Cellular DNA is constantly exposed to exogenous and endogenous DNA-damaging agents like reactive oxygen species, nitric oxide metabolites, and alkylating agents [55] leading to accumulation of mutations in the genome aggravated by loss of capacity in the DNA repair systems [56]. DNA damage is tightly linked to various ageing stresses, such as oxidative stress, telomere shortening, inflammation, irradiation, exposure to chemicals, and mitotic stress [57]. This is supported by a recent study which took repeated measurement of the relative leukocyte telomere length in a dairy herd. Higher rates of telomere attrition in individual cows was predictive of a shorter productive lifespan, suggesting a link between telomere loss and health [58].

In our dataset, NELL2 was the gene most highly correlated with lactation number in the turquoise module and was also overall the most differentially expressed gene between the PP and MP > 3 cows, with greater leucocyte expression in the older animals. The encoded protein Neural EGFL Like 2 is highly conserved in mammals and is a glycoprotein containing several von Willebrand factor C and epidermal growth factor (EGF)-like domains. This has a variety of possible roles but amongst these a cell survival-promoting effect mediated by an intracellular mitogen-activated protein kinase (MAPK) pathway has been relatively well studied in neural tissues [59, 60]. NELL2 is also important in protecting cells from death caused by endoplasmic reticulum (ER) stress resulting in the accumulation of unfolded proteins which trigger the unfolded protein response. Within this context, overexpression of NELL2 decreased expression of ER stress-induced C/EBP homologous protein (CHOP) and the pro-apoptotic caspases 3 and 7 while increasing expression of ER chaperones and anti-apoptotic Bcl-xL [61].

Another relevant gene to consider is MTOR, which encodes Mechanistic Target of Rapomycin, a protein belonging to a family of phosphatidylinositol kinase-related kinases which mediate cellular responses to stresses such as DNA damage and nutrient deprivation. This gene is a central regulator of metabolic homeostasis and is associated with lifespan in many species [62, 63]. MTOR is a component of two distinct complexes of which mTORC1 controls protein synthesis, cell growth and proliferation. Genes which are part of the MTOR pathway include the transcription factors FOXO1 and FOXP1. A GWAS for longevity, based on a population of Holstein cows, previously identified a region on Bta16 containing MTOR [64]. In our study FOXP1 expression was positively related to increasing lactation number, whereas FOXO1 expression was negatively related. Also, of potential relevance here is the opposing roles of these two FOX genes in the regulation of glucose homeostasis, through competition in binding to the insulin response element in gene promoters. Up-regulation of FOXP1 in mice inhibited the hepatic expression of key gluconeogenic genes, including PGC-1α, PEPCK and G6PC [65]. LAMTOR1, − 2 and − 3 all also featured in the list of DEG which were negatively related to lactation number in our study. These genes encode late endosomal/lysosomal adaptor, MAPK and MTOR activator-1, − 2 and − 3 respectively, all subunits of the Ragulator complex. This functions as a lysosome anchor, which recruits Rag GTPase and its associated mTORC1 complex to the lysosomal surface prior to MTOR activation [66].

Regulation of gene expression

The control of gene expression becomes more dysregulated with cellular ageing. A large number of genes and pathways identified in this study are involved in regulation of gene expression. The violet module was the most highly correlated with lactation number and contained a group of genes involved in transcriptional regulation including ZNF462, SOX13, and SOX4, all of which featured individually in the top 20 genes whose expression was most highly correlated with age (Table 1). Genes in the orange module were more highly expressed in PP cows and this module was enriched with genes involved in gene expression (CHTOP, CPSF6, THOC1, UPF3B, etc.) and metabolism of RNA (AMDHD1, SRSF2, SRSF4, SRSF5, etc.). In the yellow module IGF2BP3 was down-regulated in the older cows and was the most highly correlated gene. IGF2BP family members were initially identified as post-transcriptional regulators of IGF2. They are RNA-binding proteins which direct nuclear RNA export and translation/degradation rates, so playing a major role as regulators of the RNA life cycle. IGF2BP3 has recently risen to prominence as a potential oncogene [67]. Other genes in this module were also grouped under regulation of gene expression (TGFB3, NEUROG2, ZNF554, etc.). The yellow module also contained the gene WRN, which exhibited lower transcript abundance in MP > 3 cows. A mutation in this gene in humans causes Werner Syndrome, an autosomal recessive disorder characterized by the premature development of ageing features. The encoded protein is a member of the RecQ family of proteins and is involved in DNA replication and repair, and telomere maintenance, so playing a crucial role in genome stability [36]. Expression of WRN was similarly reduced in leucocytes of older humans [18]. ADARB1 was another candidate gene from previous studies which was up-regulated in the older MP > 3 cows. It encodes a deaminase enzyme with A-to-I RNA editing activity, was previously identified in a study of men aged 90–119 years, and is also associated with longevity in C. elegans [35].

Metabolism

One key difference between PP and MP > 3 is that milk production potential in dairy cows increases with age [5]. The liver coordinates the extensive metabolic changes required for milk production and these are reflected in circulating metabolite concentrations [11, 38]. Milk synthesis has a high requirement for glucose. In ruminants this demand is met almost exclusively through hepatic gluconeogenesis and cows are at risk of glucose insufficiency during early lactation, the time period we investigated [16]. NEFA are released from lipid stores as an alternative energy source and are either used by the udder to provide milk triglycerides, fully oxidized in the liver to provide energy, or partially oxidised resulting in the production of ketone bodies, in particular BHB. Circulating BHB concentrations are thus an index of fatty acid oxidation and concentrations are significantly higher in older cows [12]. BHB, NEFA and glucose concentrations can all influence leukocytes directly. Immune cells require an adequate supply of nutrients including glucose to mount an effective immune defense [68]. Neutrophils from cows with more elevated NEFA and BHB concentrations after calving had reduced expression of genes important for granulocyte recruitment, IFN signaling and apoptosis [69]. This suggested that neutrophil survival time was longer in the circulation when exposed to pro-inflammatory conditions. Another study of the circulating leucocyte transcriptome in early lactation found that expression of genes in KEGG pathways relating to DNA replication, cell cycle, homologous recombination, base excision repair, and valine, leucine, and isoleucine biosynthesis were all inhibited as plasma BHB increased, whereas genes involved with vitamin metabolism, the endocrine system, signalling molecules and the immune system were activated [70].

In this study genes in both black and turquoise modules were negatively significantly correlated to cow parity and David functional annotation cluster analysis found enrichment of the terms fatty acid metabolism and glycolysis. Interestingly, we found that the majority of 27 genes in Cluster 4, with higher expression in MP > 3 cows, were involved in the Glycolysis/Gluconeogenesis pathway (Fig. 3). Although the underlying metabolic background is very different, up-regulation of a small cluster of genes relating to “Fatty acid metabolism, peroxisome activity” was also associated with ageing in the human blood transcriptome [18].

Other key genes associated with lactation number

Within the violet module LAMA4 was the most highly correlated gene overall, with greater expression in the PP cows. This encodes a subunit of laminin, part of the family of extracellular matrix glycoproteins which form the major non-collagenous constituent of basement membranes. Laminin is thought to mediate the attachment, migration and organization of cells into tissues by interacting with other extracellular matrix components, suggesting that these activities may be reduced in older cows. There is a second MTOR complex mTORC2, which acts as a regulator of the actin cytoskeleton. This is a network of actin and actin binding proteins which are important for a range of essential cellular processes including organelle transport, cell migration, phagocytosis, and cell cycle progression. This is in agreement with the blood transcriptomic study in humans [18] in which expression of a cluster of genes relating to the actin cytoskeleton and focal adhesion also increased with ageing (ACTA2, ACTN4, CSRP1, ILK, LPP, TAGLIN, TLN1, VCL and WDR1). In contrast, ARHGAP15 was more highly expressed in younger cows (this study) and in humans [18]. This gene encodes Rho GTPase activating protein 15 which is involved in T- and B-cell signalling and promotes an increase in actin stress fibres and cell contraction [71].

Conclusions

The samples collected for this study provided the opportunity to analyse the transcriptomic profile of blood leukocytes in a large number of early lactation cows and we were thus able to capture the complex and temporally dynamic biological pathways which alter as cows age. We have, we believe for the first time, identified genes and pathways associated with increasing lactation number in cows and found that many were comparable to those known to be associated with ageing in humans and model organisms. Immune-related pathways were differentially expressed between primiparous and older cows, including genes involved in innate and adaptive immunity, with many immune defence genes being more highly expressed in MP > 3 animals. In addition, we found changes in mitochondrial and metabolic function, ribosome biogenesis, transcriptional regulation and DNA replication, elongation and repair. Pathways supplying leukocytes with energy changed in cows with different lactation numbers, with increased expression of genes encoding enzymes involved in beta-oxidation of fatty acids in the PP cows whereas genes involved in glycolysis were up-regulated in the older cows. These changes may be particularly relevant to understanding how dairy cows respond to metabolic stress during early lactation, when they are short of energy and there is competition for nutrient supply between the demands of milk production and the need for immune defence. An improved understanding of these processes may help dairy farmers to improve both genetic selection and nutritional management to increase longevity, so benefitting agricultural production and animal welfare.

Methods

GPlusE study design and sample collection

The samples used in this study were collected as a part of the Genotype plus Environment (GplusE) FP7-Project (http://www.gpluse.eu). The animals were located on six experimental dairy farms belonging to members of the consortium (Supplementary Table 1). One 3 ml blood sample from each cow was drawn from the tail vein into a Tempus blood collection tubes (Thermo Fischer, UK) in early lactation, at around 14 days after calving and stored at − 20 °C. This sample was taken with the informed consent and ethical approval of the organisations involved and complied with the relevant national and EU legislation under the European Union (Protection of Animals used for Scientific Purposes) Regulations 2012 (S.I. No. 543 of 2012). All cows were subsequently released back into the herd. Details of the management of each herd and their average milk yields over the first 50 days in milk are provided in Krogh et al. [72].

A total of 229 adult Holstein Friesian cows ranging between parities 1 and 7 were recruited. These were subsequently divided into three lactation number groups for analysis: i) primiparous (PP, n = 53), ii) multiparous in lactations 2–3 (MP 2–3, n = 121), and iii) MP in lactations 4–7 (MP > 3, n = 55).

RNA isolation, library preparation, and sequencing

The individuals performing the sample processing were blinded to the lactation groups. Total RNA was isolated from whole blood leukocytes using the Tempus Spin RNA isolation Kit (Thermo Fischer, Loughborough, UK) following the manufacturer’s instructions. Sample purity and RNA quantity were measured using both a NanoDrop 1000 (Thermo Fischer, UK) and an Agilent BioAnalyzer 2000 using the Agilent RNA 6000 Nano Kit (Agilent, Cheadle, UK; Supplementary Table 6). No samples were excluded based on RNA quality. Library preparation was conducted at the University of Liege (GIGA Research Facility, Liege, Belgium), using 0.75 μg of total RNA with the Illumina TruSeq Stranded Total RNA Ribo-Zero Gold Sample Preparation kit (Illumina, San Diego, USA) and sequenced on the Illumina NextSeq 500 sequencer, producing on average 30 million single end reads of 75 nucleotides length per sample.

Transcriptomic analyses

Reads with poor quality were trimmed or removed according to base quality using Trimmomatic v.0.36 [73]. The quality of raw and trimmed FASTQ files was assessed with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The latest version of the Bos taurus assembly (ARS-UCD1.2), and its corresponding gene set, was used as reference to map reads using the splice aware aligner HISAT2 [74]. Next, SAM files were converted to BAM files and coordinate sorted with SAMtools [75]. BAM files were processed further with Picard Tools (http://picard.sourceforge.net/) in order to remove PCR duplicates, add read group information, sort by chromosome and create indexes. Reads per gene were counted with StringTie [76].

Differential leukocyte gene expression analysis between the three lactation number groups was conducted with the package DESeq2 [77]. Herd effect was considered and removed using limma remove batch effect (limma::removeBatchEffect; Supplementary Fig. 3). Weighted gene co-expression network analysis (WGCNA, R package [78]; was used to construct a co-expression network on the DESeq2 normalized data. WGCNA follows a 6-step process to predict which genes are connected to each other, cluster them into gene networks and test which gene networks are associated with phenotypic status, leading to the selection of hub genes.

Genes with variance < 0.05 were filtered out, and the results (total of 13,769 genes) were used as input to the signed WGCNA network construction. We generated a “signed weighted” adjacency matrix for downstream analyses in which the direction of a pair of genes’ correlation is preserved, and a positive correlation may indicate “activation” whereas a negative correlation may indicate “repression”. The adjacency matrix reported a correlation value between every pair of gene expression measurement across all 229 samples. The next step was to raise the adjacency matrix to a software-determined exponential power, thereby reducing noise by pushing weaker pairwise connection values closer to zero relative to stronger values. The exponential power used is the lowest value needed to ensure the network approximates scale-free topology and was set to seven. The adjacency matrix was then transformed into a “topological overlap” matrix by calculating topological overlap (TOM) scores for each gene. This score accounts for each pair of genes’ connection strength (adjacency value) to each other as well as their connection strengths (adjacency values) to every other gene in the adjacency matrix. WGCNA identifies gene co-expression networks via average linkage hierarchical clustering using a TOM-based dissimilarity measure as a measure of distance. The resultant dendrogram of clustered genes is segregated into individual modules with at least 35 genes using WGCNA’s dynamic tree-cutting algorithm. WGCNA calculates each module’s “eigengene” (first principle component), using all samples’ gene expression values for all genes in that module. A module eigengene is considered a summarized expression profile representative of that module for all samples. Finally, each module eigengene was tested for statistical association to the phenotypic trait of lactation number (PP, MP2–3 and MP > 3).

Cattle genes based on the Ensembl v.100 database were retrieved using the Ensembl BioMart online tool [79] and gene ontology (GO) analysis was performed using genes found in significantly correlated modules to trait. Moreover, a David functional annotation cluster analysis [27] was performed; all 13,769 genes were used as background against all significantly correlated modules.

Additional file 1: Additional tables supporting the main text of the manuscript. Referred to as Supplementary Tables 1–6. Supplementary Table 1. Milk yield data (kg/d) for the six dairy herds used in the study, organised by herd and parity. Supplementary Table 2. Mapping rate of all blood samples against the Bos taurus reference ARS-UCD1.2 and summary statistics. Supplementary Table 3. Gene ontology enrichment analysis of genes in significantly correlated modules. Supplementary Table 4. David functional cluster analysis of genes in significantly correlated modules. Supplementary Table 5A, 5B. (i) Top six most significant clusters; (ii) Number of genes highly expressed in the PP/MP3 cows per cluster. Supplementary Table 6. RNA concentration measured with Nanodrop, Bioanalyzer, and RIN values.

Additional file 2: Supplementary Fig. 1. Additional figures supporting the main text of the manuscript. Referred to as Supplementary Figs. 1–3. Supplementary Fig. 1. Volcano plot of differentially expressed genes (padj < 0.01, LFC > 1). Supplementary Fig. 2. Outliers detected by sample clustering (WGCNA). Supplementary Fig. 3. Principal component analysis of herd effect. A) with herd effect, B) herd effect removed.