Age Related Changes in the Circulating Leukocyte Transcriptome in Dairy Cows

Background: Previous studies of ageing have identi�ed many pathways which are consistently altered in humans and model organisms. 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 bene�t 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 age 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 signi�cantly correlated (r ≥ 0.25) to the trait age. 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 �ngers. 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 identi�ed 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 age in cows were identi�ed for the �rst


Background
Recent studies in human populations have highlighted age-associated changes in leucocyte functionality affecting innate and adaptive immune functions [1].The causes and consequences of ageing on the human blood transcriptome have, however, proved di cult to dissect due to interactions with environmental in uences, genetic factors and a large number of age-related diseases [2].Studies on model organisms have highlighted that ageing is characterized by many alterations at molecular, cellular and tissue level [3].Studies of the ageing transcriptome have been performed in species including C. elegans [4], ies [5], rodents [6] and humans [7].This approach has identi ed various signatures found to occur repeatedly across different tissues and organisms.Candidate genes whose expression is consistently associated with cellular ageing have been classi ed 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 [8].
Although all living creatures age at some point, our knowledge on the biology of ageing is still not su cient.The physiological process of ageing in humans is associated with a progressive loss of function and increased vulnerability to disease, frailty, and disability [9].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 [10].
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 [11,12].To maximise economic potential, heifers should rst calve at 24 months of age so becoming primiparous cows, starting their rst lactation, and beginning to pay back their rearing costs through the production of saleable milk [13].Ideally, they should continue to calve at annual intervals, with their milk production potential increasing progressively until the fourth or fth lactation [14].Cows which survive for longer achieve greater lifetime milk production associated with higher pro tability [12,15] together with the bene t of reduced greenhouse gas emissions [16].Many cows do not achieve optimum longevity due to premature culling, for which the main reasons are mastitis, infertility and lameness [17,18].There are major alterations in the metabolic pro les associated with the start of lactation [19].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 nonesteri ed fatty acids (NEFA) and beta-hydroxybutyrate (BHB) [20].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 [21].Both metabolic and immune dysfunction therefore impact on the transcriptomic changes in leukocytes in early lactation.
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 signi cant 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.
These 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 1).As we were most interested in understanding the effect of ageing, we contrasted the leukocyte expression between the two extreme age groups, PP vs MP > 3, using those classi ed as MP2-3 as background.Firstly, we performed a differential expression analysis using DESeq2.Over 5,000 genes out of the 17,216 genes initially mapped were found to be differentially expressed between age groups following Benjamini-Hochberg false discovery rate adjustment (padj < 0.1).This list was narrowed down to 2,925 differentially expressed genes (DEG) when contrasting young vs old cows (Supplementary Fig. 1).We then applied WGCNA to identify gene modules from our dataset associated to the trait of age.From an initial total of 21,207 genes, 13,769 (64.9%) genes passed the DESeq2 ltering steps.Samples from three cows were removed at this stage as outliers (two from MP3 and one from MP2-3 age group, respectively; Supplementary Fig. 2), leaving a total of 226 samples.This unsupervised technique identi ed 32 interconnected gene modules from the ltered gene list (Fig. 1A).The number of genes per module is shown in Fig. 1B.Of these, the violet and turquoise modules were signi cantly related to the trait age in both the PP and MP > 3 age groups (Fig. 1C).As expected, modules with a signi cant positive correlation in the younger cows were negatively correlated with age in the older ones, and vice versa, while no modules showed a signi cant correlation with the medium age group.Within these seven signi cant modules, 2,274 genes were more highly expressed in the leukocytes of the PP cows and 2,721 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.

Signi cant modules in leukocytes of both PP and MP > 3 cows
The violet module containing 87 genes had the highest correlations with age (+ 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 signi cance = 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).
The turquoise module represents the second module which had signi cant correlations in opposite directions with PP and MP > 3 cows, of -0.29 and 0.33, respectively.This module contained 1,818 genes of which NELL2 was the most highly correlated with age (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 1,818 genes with a wide variety of functions (Fig. 2B, Supplementary Table 2).The top six most signi cant 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 signi cantly positive in PP cows with correlations spanning from 0.27 to 0.25, implying lower expression in leukocytes as cows aged (Fig. 2A, Supplementary Table 2).The most signi cant GO terms in the midnightblue module were all biological processes: protein autophosphorylation, immune system process, cell activation and regulation of cellcell 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 signi cantly up-regulated in MP > 3 (correlation of 0.24).Genes were therefore more highly expressed in the older cows.GO enrichment analysis of genes in this module (Fig. 2B, Supplementary Table 2) showed that the most signi cant 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 signi cant.
Yellow Module: down-regulated in leukocytes of MP > 3 cows Lastly, the yellow module was signi cantly down-regulated in leukocytes as cows aged (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 2).

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 [22].All the results are provided in Supplementary Table 3.
The top six most signi cant clusters in each category are summarised in Table 2, with the genes in each of these clusters listed in Supplementary Tables 4A and 4B.
Table 2 Summary of the six most signi cant clusters obtained from DAVID functional annotation cluster analysis of the DEG in the modules which were signi cantly positively or negatively correlated to age 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 • Krueppel-associated box The most signi cant term in Cluster 1 was Krueppel-associated box.All but ve of the Cluster 1 genes were included in the much larger Cluster 2 with 142 genes, for which the most signi cant terms were Zinc-nger and metal ion binding.Cluster 2 also contained all 18 genes from Cluster 3, for which the three signi cant 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 signi cant term was BTB/POZ fold.Finally, Cluster 6 contained 50 genes with the most signi cant term being Immunoglobulin-like domain.Overall, this analysis shows that the most signi cant functions which were more highly expressed in the younger animals 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, Disul de 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).
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 age was Cluster 16 (Supplementary Table 3).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.

Comparison with other species
Finally, we investigated the overlap between genes signi cantly associated with age in leucocytes in our study and 170 candidate genes identi ed from previous studies of ageing in humans and model organisms (mainly mice and C. elegans [2]).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 signi cantly correlated modules, indicating reduced or increased expression in older cows.Of these 16 of the 38 DEG were previously identi ed as being age-related in the human blood transcriptome ( [2]; Table 3).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 [23].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 bene t 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 [17,18].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 pro tability of the dairy industry.For example, one genome wide association study (GWAS) which investigated longevity found genes such as NPFFR2, previously identi ed as a candidate gene for mastitis resistance and two zinc nger proteins (ZNF613, ZNF717), which have been associated with calving di culties [24].Information about the age-related morbidities and causes of death that a ict 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 rst, to our knowledge, to explore age-related changes in a cattle population by measuring global gene expression in their leukocytes.For this we compared rst 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 [19,25] 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.Despite the many differences in physiology, lifestyle and lifespan between species, many of the genes and pathways identi ed in cows were the same as those highlighted in previous studies of ageing in other species [2][3][4]7].
All the samples analysed were whole blood leukocytes.Many aspects of the immune system alter during ageing, eventually leading to immunosenescence [8].During this process different immune cell subsets are affected in different ways [1].There is also an increase in baseline systemic in ammation with age, termed "in ammaging" [26,27].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 [2,28,29].The leukocyte transcriptome provides a re ection 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.The age-related changes reported here will, to some extent, re ect altered abundance of different cell types as well as their changing expression patterns.This is partly mitigated by the use of WGCNA, which identi es networks of co-expressed genes whose expression is highly correlated.The use of whole blood also avoids the pitfall of potential artefacts which can be induced during cell separation procedures.
Previous work based on studies of human populations and model organisms has identi ed a number of transcriptional signatures for cellular ageing which occur repeatedly across different tissues and organisms and which segregate into six main groups [8,30].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 for all of these 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 [31].Oxidative metabolism causes endogenous production of free radical molecules and oxidative damage accumulates in multiple tissues and species with age [30].
Ageing in the model species Drosophila melanogaster was associated with down-regulation of numerous mitochondrial genes, including those involved in the electron-transport-chain and mitochondrial metabolism [28,32].Knockdown of mitochondrial ribosomal proteins was shown to reduce mitochondrial respiration and activate the mitochondrial unfolded protein response [33] while accumulated mutations in somatic mitochondrial DNA (mtDNA) and respiratory chain dysfunction were associated with ageing in mice [34].In our study, genes encoding proteins involved in fatty acid beta-oxidation were more highly expressed in leukocytes in PP cows.There were also age-related 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 [30].In our study, GO enrichment analysis of all the genes in the violet module, which were more highly expressed in the younger cows, revealed their involvement in protein autophosphorylation and catalytic activity (enzyme activity).Autophosphorylation is a type of posttranslational modi cation 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 identi ed in the yellow module as being down-regulated in the older MP > 3 cows.This is in line with the expectation that leukocytes in ageing cows might have an overall downregulation of protein synthesis.

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 [9].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 [1].Our results are in accord with the study by Peters et al [2], who investigated leucocytes from ageing human populations, in nding 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 signi cant overlap between the DEG with immune functions and candidate age-related genes in the human transcriptome [2].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 ADARB1encodes 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 older cows are more actively engaged in combatting disease pathogens through activation of the innate immune system and also support a higher level of in ammation with ageing.Parturition is itself an in ammatory process and older cows are also more vulnerable to metabolic disorders after calving.Our previous study showed that metabolic disorders led to prolonged uterine in ammation by up-regulating the genes and pathways related to immune and in ammatory processes [35].

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 [3].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 [36].The turquoise module was the second highest correlated module with ageing, 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 3).The latter pathway is interlinked with metabolic pathways to ensure coordinate regulation and ne-tuning of cellular metabolic responses in line with cellular energy status, nutrient availability and hormonal/growth factor signalling input [36].

Stress and DNA damage
Accumulation of genetic damage represent 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 [37] leading to accumulation of mutations in the genome aggravated by loss of capacity in the DNA repair systems [38].DNA damage is tightly linked to various ageing stresses, such as oxidative stress, telomere shortening, in ammation, irradiation, exposure to chemicals, and mitotic stress [39].In our dataset, NELL2 was the most highly correlated gene with age 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 mitogenactivated protein kinase (MAPK) pathway has been relatively well studied in neural tissues [40,41].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 [42].
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 [43,44].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 identi ed a region on Bta16 containing MTOR [45].In our study FOXP1 expression was positively related to increasing age, 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 [46].LAMTOR1, -2 and − 3 all also featured in the list of DEG which were negatively related to age 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 [47].

Regulation of gene expression
The control of gene expression becomes more dysregulated with cellular ageing.large number of genes and pathways identi ed in this study are involved in regulation of gene expression.The violet module was the most highly correlated with age and contained a group of genes involved in transcriptional regulation including ZNF462, SOX13, and SOX4, all of which featured 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 identi ed 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 [48].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 older 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 [23].Expression of WRN was similarly reduced in leucocytes of older humans [2].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 identi ed in a study of men aged 90-119 years, and is also associated with longevity in C. elegans [49].

Metabolism
One key difference between PP and MP > 3 is that milk production potential in dairy cows increases with age [14].The liver coordinates the extensive metabolic changes required for milk production and these are re ected in circulating metabolite concentrations [19,25].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 insu ciency during early lactation, the time period we investigated [21].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 signi cantly higher in older cows [20].BHB, NEFA and glucose concentrations can all in uence leukocytes directly.Immune cells require an adequate supply of nutrients including glucose to mount an effective immune defense [50].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 [51].This suggested that neutrophil survival time was longer in the circulation when exposed to pro-in ammatory 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 [52].
In this study genes in both black and turquoise modules were negatively signi cantly correlated to the trait age in PP 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 older 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 [2].
Other key age-related genes 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 [2] 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 [2].This gene encodes Rho GTPase activating protein 15 which is involved in T-and B-cell signalling and promotes an increase in actin stress bres and cell contraction [53].

Conclusions
The samples collected for this study provided the opportunity to analyse the transcriptomic pro le of blood leukocytes in a large number of early lactating cows and we were thus able to capture the complex and temporally dynamic biological pathways potentially involved in the ageing process in this species.
We have, we believe for the rst time, identi ed genes and pathways associated with age 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 older animals.In addition, we found age-associated changes in mitochondrial and metabolic function, ribosome biogenesis, transcriptional regulation and DNA replication, elongation and repair.Pathways supplying leukocytes with energy changed with age, with higher 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 bene tting agricultural production and animal welfare.

GPlusE study sample collection
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 Scienti c 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 are provided in Krogh et al. [54].
A total of 229 adult Holstein Friesian cows ranging between parities 1 and 7 were recruited.These were subsequently divided into three age groups for analysis: i) primiparous (PP, n = 53), ii) multiparous in

RNA isolation, library preparation, and sequencing
The individuals performing the sample were blinded to the age 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).Library preparation was conducted at the University of Liege (GIGA Research Facility, Liege, Belgium), using 0.75 ug 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 [55].The quality of raw and trimmed FASTQ les 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 [56].Next, SAM les were converted to BAM les and coordinate sorted with SAMtools [57].BAM les 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 [58].
Differential leukocyte gene expression analysis between the three age groups was conducted with the package DESeq2 [59].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; [60] 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 ltered 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 identi es 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" ( rst principle component), using all samples' gene expression values for all genes in that module.A module eigengene is considered a summarized expression pro le representative of that module for all samples.Finally, each module eigengene was tested for statistical association to the phenotypic trait of age (PP, MP2-3 and MP > 3).

••
Zinc nger C2H2-type/integrase DNA-binding domain 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 nger proteins.

Figures
Figures

Figure 2 GO
Figure 2

Figure 3 David
Figure 3

Table 1
Top 20genes 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.

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.