The effect of environmental enrichment on immune cell DNA methylation profiles depends on the parity of sows

The effect of environmental enrichment on immune cell DNA methylation profiles depends on the parity of sows | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results The effect of environmental enrichment on immune cell DNA methylation profiles depends on the parity of sows Mariana Mescouto Lopes , View ORCID Profile Gabriel Costa Monteiro Moreira , Aurélie Chaulot-Talmon , Anne Frambourg , View ORCID Profile Valentin Costes , View ORCID Profile Julie Demars , View ORCID Profile Elodie Merlot , View ORCID Profile Hélène Jammes doi: https://doi.org/10.1101/2025.07.01.662510 Mariana Mescouto Lopes a PEGASE, INRAE, Institut Agro , 35590 Saint-Gilles, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Gabriel Costa Monteiro Moreira b BREED, Université Paris-Saclay, UVSQ, INRAE , 78352 Jouy-en-Josas, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Gabriel Costa Monteiro Moreira Aurélie Chaulot-Talmon b BREED, Université Paris-Saclay, UVSQ, INRAE , 78352 Jouy-en-Josas, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Anne Frambourg b BREED, Université Paris-Saclay, UVSQ, INRAE , 78352 Jouy-en-Josas, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Valentin Costes b BREED, Université Paris-Saclay, UVSQ, INRAE , 78352 Jouy-en-Josas, France c ELIANCE , 149 rue de Bercy, 75595 Paris cedex 12, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Valentin Costes Julie Demars d GenPhySE, INRAE, Université de Toulouse , 31326 Castanet Tolosan, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Julie Demars Elodie Merlot a PEGASE, INRAE, Institut Agro , 35590 Saint-Gilles, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Elodie Merlot For correspondence: helene.jammes{at}inrae.fr elodie.merlot{at}inrae.fr Hélène Jammes b BREED, Université Paris-Saclay, UVSQ, INRAE , 78352 Jouy-en-Josas, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Hélène Jammes For correspondence: helene.jammes{at}inrae.fr elodie.merlot{at}inrae.fr Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract The aim of this study was to identify epigenetic markers that reflected positive affective states in multiparous pregnant sows. The animals were housed during gestation in either a conventional (C) environment (2.4 m² per sow), featuring a concrete slatted floor and minimal enrichment, or in an enriched (E) environment (3.5 m² per sow) with deep straw bedding. Peripheral blood mononuclear cells (PBMCs) were isolated at day 98 of gestation (G98) and day 12 of lactation (L12) for genomic DNA extraction and reduced representation bisulfite sequencing (RRBS). Significant effects of individual class of parity (low-parity (LP) versus high-parity (HP)) were observed with the identification of 1,358 and 680 differentially methylated cytosines (DMCs) at G98 and L12, respectively. Interestingly, some of them are organized into differentially methylated regions (DMRs). Some DMRs colocalized with or near a gene and displayed a continuous methylation distribution (44 at G98 and 15 at L12 including 4 in common): 5 targeted genes are related to epigenetic regulation ( DNMT3A , KDM8 , HDAC4 , SIRT2 and U2 ) and 12 to immune function ( CD2 , CD5 , CAMK4 , SECTM1 , URODL1 , CMIP , SEC14L1, SKI , TNFRSF1B , CCND3 , SGK1 , PACS1 ). These results suggested a true epigenetic impact of parity class on individual immunity. Considering the two parity classes separately, a minor effect of environmental enrichment was observed (at G98: 60 and 42 DMCs; at L12: 35 and 81 DMCs, in the LP and HP groups, respectively). Remarkably, some DMC-associated genes had previously been linked to affective states in humans. In conclusion, unexpected DNA methylation changes associated with parity class were identified, suggesting a genome adaptation during reproductive life and modifying the response to housing conditions. Furthermore, specific CpG sites emerged as potential biomarkers of positive affective states in pigs. Highlights This study revealed a substantial impact of the genetic background of sows on DNA methylation patterns, emphasizing the need to account for genetic factors when analysing epigenetic data. Parity, which was confounded with the age of the sows in this study, significantly influenced DNA methylation, even among multiparous sows, highlighting the importance of this factor in modulating epigenetic mechanisms in immune cells. We observed minor differences in the effects of environmental enrichment on DNA methylation profiles, which differed depending on the parity group. We suggest a short list of possible biomarkers with biological meaning responding to environmental enrichment. Introduction Affective states, or the subjective positive or negative feelings experienced by animals, play a central role in their welfare 1 – 3 . Given the inability of non-human animals to communicate their feelings verbally, scientists are reliant on indirect behavioural and physiological tests to assess their affective states 4 . While tools are now available to identify responses to acute emotions, it is more challenging to assess long-term, and usually less intense, emotions (i.e., moods). Animal welfare science would therefore benefit from the development of reliable, feasible, and practical physiological markers that reflect the affective states of animals. In mainstream pig production systems, one of the primary sources of animal stress is the barren housing environment in which they are raised. The lack of adequate stimuli limits the expression of natural behaviours and cognitive-sensory capacities, leading to signs of chronic stress and poor welfare. Environmental enrichment has been shown to positively impact pig behaviour 5 , 6 , cognitive function 7 , 8 , and physiological variables such as cortisol levels 9 , 10 . Enriched housing conditions have also been found to modulate the blood immune response and reduce disease susceptibility in pigs 11 , 12 . Blood immune cells may therefore offer valuable insights into the influence of living environments on the affective state of animals. Epigenetic mechanisms involve stable and heritable modifications during cell division, which regulate gene expression without altering nucleotide sequences 13 . These mechanisms play a crucial role in development and cell differentiation and enable the integration of environmental factors to shape cellular function. For instance, DNA methylation, which encompasses the addition of a methyl (CH3) group to a cytosine residue within a CG dinucleotide, is essential to establishing and preserving cellular identity. During haematopoiesis, from stem cells to highly differentiated immune cells, specific DNA methylation profiles are established and support the unique functions of each immune sub-type cells 14 , 15 . Despite this specificity, DNA methylation in blood immune cells can be influenced by various internal and external factors throughout an individual’s life, from early embryonic development to death. In humans and rodents, changes of DNA methylation patterns in blood immune cells have been associated with age 16 , physiological stages 17 , toxic exposure 18 and prenatal stress 19 . In pigs, specific DNA methylation profiles have been described for different immune cell subtypes 20 . Different factors such as weaning 21 , viral exposure 22 or metabolic status 23 , 24 influence DNA methylation in peripheral blood mononuclear cells (PBMC). Furthermore, DNA methylation profiles in blood immune cells have been used effectively to accurately detect and characterize age-related biomarkers in various vertebrate species, including pigs 25 . While the neuro-epigenome of brain regions in piglets has been analysed in terms of stereotypies and enriched maternal environments 26 , the direct effects of environmental enrichment on the DNA methylation of adult PBMC remained unexplored. The aim of this study was therefore to investigate the PBMC methylome of sows of different ages and parity ranks housed under two contrasted environments that could potentially have an impact on the affective states of these animals. Materials and methods Experimental design and animals The experiment was conducted at the Chambre Régionale d’Agriculture de Bretagne Experimental Farm (Saint-Nicolas-du-Pélem, France), as detailed in our published study 27 , 28 . Young nulliparous sows are randomly assigned upon purchase to one of two independent reproductive units, which feature distinct housing environments 27 , 28 . They were kept under two separate housing environments: conventional (C, n = 36) and enriched (E, n = 35) reproductive units. In the C unit, the sows were housed in a pen with a concrete slatted floor and minimal enrichment materials, with a minimum space allowance of 2.4 m 2 per sow. In the E unit, they were accommodated on deep straw bedding with more space per animal (3.5 m 2 ). Pregnant sows resulting from crosses between Landrace and Large White breeds, were housed together in groups since their previous weaning until 10 days before their expected farrowing date. Pregnant sows in both units were of different parity ranks (0 to 7 gestations). On gestation day 105 (G105; approximately 10 days before expected farrowing), sows from both units were transferred to individual farrowing pens with slatted floors. They remained there with their litter for four weeks until weaning, after which the sows returned to their previous E or C units. The sows therefore remained in their assigned housing environment throughout their reproductive life. All sows, whatever the housing conditions, were fed individually via an automatic dispensing system, with standardized management and feeding practices across units. All details regarding the animals were reported in Table 1 . The trial was carried out over the course of one year, consisting of two replicates: R1 (from November 2021 to February 2022) and R2 (from March to July 2022). In each replicate, we selected 7-8 multiparous sows per environment for blood sampling, totalling 29 animals (C, n = 7 in R1, and n = 8 in R2; and E, n = 7 in R1 and n = 7 in R2). Particular attention was paid to the sows’ genealogy. The number of individuals included in the experimental design is required for methylome analyses, given that modest environmental effects are expected. View this table: View inline View popup Download powerpoint Table 1: Physiological parameters and genealogy of sows selected for blood sampling. Body weight and age were recorded at the blood sampling time point on G98. A total of 27 different sows were used: 25 different sows (n =13 in C and n =12 in E) were included in the trial; and two additional sows (one in each environment group) were used twice successively into the trial (C_01/C_11 and E_16/E_25). Abbreviations: Rep, replicate assay; C, conventional housing; E, enriched housing. F1 to F16 and M1 to M24: identity of the fathers and mothers of the experimental sows, respectively; DF and DM: different father and different mother, respectively. PHS and MHS: paternal half-sister and maternal half-sister, respectively). Blood sampling, PBMC preparation and genomic DNA preparation The experimental timeline is summarized in Figure 1A . Blood samples (10 mL blood collection tubes at a final concentration of 1.8 mg EDTA (Ethylene Diamine Tetra Acetic acid) per mL of blood; three tubes collected per sow) were collected by an experienced veterinarian from the anterior vena cava at 98 th day of the gestation period (G98, before transfer to individual farrowing pens) and from the mammary vein at the beginning of lactation (12 th day; L12), all from the same animals. During experimental replicate R1, sows from enriched housing group were blood sampled first, and inversely, during experimental replicate R2, sows from conventional housing were collected first. These samples were transported in a thermo box with refrigerated rolls and arrived at the laboratory within three hours of the sampling procedure. Each blood sample was notified with sow identification number: a 6-digit number totally independent of age, parity or experimental treatment of the animals, enabling the samples to be processed blindly. Blood formula was immediately determined. Total numbers and percentages of lymphocytes, monocytes and neutrophils were determined with an MS-9® hematology automatic cell counter (Melet Schloesing laboratories, Osny, France; see all data in 27 , 28 ). A part of blood samples was used for the isolation of peripheral blood mononuclear cells (PBMC) according to a standard protocol 29 . Aliquots of 2 million cells were prepared with 0.2 mL lysis buffer (SLB; 10 mM Tris-HCl, pH 7.5, 10 mM EDTA, 10 mM NaCl, 0.2% SDS), thrown into liquid nitrogen for immediate freezing, and then conserved at −80°C until use. Download figure Open in new tab Figure 1. Flowchart and study overview. (A) Schematic outline of the experimental design. (B) Schematic overview of steps from PBMC preparation to the quality control of sequencing data. (C) Schematic overview of CpG background generation and analysis. Abbreviations: G, gestation day; L, lactation day; PBMC, peripheral mononuclear cells; RRBS, reduced representation bisulfite sequencing; CpG 10-500 , CpG with a minimum of 10 and maximum of 500 uniquely mapped reads; SNP, single nucleotide polymorphism. Genomic DNA (gDNA) extraction was performed following a conventional procedure that included: proteinase K digestion (2mg/mL final concentration) in SLB lysis buffer overnight at 55°C, RNA elimination (RNAse incubation; 0.2mg/mL, 37°C 2h) and two phenol/chloroform extractions. After precipitation with 0.2M NaCl in 70% ethanol (overnight at −20°C and centrifugation, 12,000 x g for 30 minutes at 4°C), the gDNA pellet was washed, dried and resuspended in TE buffer (10 mM Tris-HCl, pH 7.5, 0.1mM EDTA). Integrity (electrophoretic migration on agarose gel), purity (Nanodrop analysis; Thermo Fisher Scientific, Les Ulis, France) and concentration (QuBit fluorometer; Thermo Fisher Scientific) were measured for all gDNA ( Figure 1B ). Reduced Representation Bisulfite Sequencing (RRBS) library preparation In this study, 58 reduced representation bisulfite sequencing (RRBS) libraries were prepared on a semi-automated machine (Genomic STARlet; Hamilton) according to the procedure previously described by Gu et al. 30 and adapted by Perrier et al. 31 . gDNA (200 ng) was cleaved by MspI enzyme at CCGG sites, generating fragments with CG dinucleotides at both ends. After the end repair step and Illumina adaptor ligation, a size selection of 150-400 bp of DNA fragments (40 - 290 bp genomic DNA fragments + adapters) was performed using SPRIselect magnetic beads (Beckman-Coulter, Villepinte, France). DNA fragments were then converted twice with sodium bisulphite using the EpiTect bisulphite kit (Qiagen, Les Ulis, France) according to the manufacturer’s instructions, and amplified with Pfu Turbo Cx hotstart DNA polymerase (Agilent, Les Ulis, France) using 14 PCR cycles and Illumina index. Aliquots of PCR products were analysed by electrophoresis on a 4-20% precast polyacrylamide TBE gel (Invitrogen; ThermoFisher scientific, Les Ulis, France) and stained with SYBR green to confirm the homogeneous pattern for all libraries. The libraries were sequenced on an Illumina NovaSeq 6000 sequencer to produce a minimum of 30 million 100 bp paired-end reads (Integragen-OncoDNA, Evry, France). Data processing The data processing method is summarized in Figure 1C . Raw reads underwent quality control using FastQC (v0.12.0; https://www.bioinformatics.babraham.ac.uk/projects/fastqc ) and subsequently trimmed using Trim Galore! (v0.4.5; https://www.bioinformatics.babraham.ac.uk/projects/trim_galore ) to eliminate reads with a Phred score below 20, shorter than 20 nucleotides, and adaptor sequences. Processed reads were then aligned to the in-silico bisulphite-converted Sus scrofa reference genome (Ensembl: GCF_000003025.6, Sscrofa11.1) using Bismark v0.20.0 in the default mode with Bowtie 1.2.1 32 . On average, 40 million raw reads per library were generated, with mean total mapping rate of 70% and a mean rate of uniquely mapped reads of 53% covering the pig genome. The bisulphite conversion rate was estimated from the unmethylated cytosine added during the end-repair step. The conversion rate reached 100% in nearly all samples. An overview of sequencing parameters for all libraries is presented in Suppl. Table 1, Sheet 1. Selection process of a common CpG 10-500 background filtered for SNP Only CpGs covered by 10 to 500 uniquely mapped reads (termed CpG 10-500 ) were retained, forming the initial raw background of 3,418,048 CpGs 10-500 , which represents approximately 10% of all CpGs in the pig genome. From this set, only CpG 10-500 with a methylation value calculated for all the libraries from the G98 and L12 datasets were retained, yielding a background of 1,473,450 CpG 10-500 . To avoid the confounding effects of sequence polymorphisms and bisulfite conversion, the CpGs 10-500 co-localized with a putative sequence polymorphism affecting the C and/or G were filtered out. For this purpose, the common background (1,473,450 CpG 10-500 ) was thus subjected to further filtration using two SNP lists to identify and minimize confounding effects of the animals’ genetic background. The first list was obtained from whole genome sequencing data (accession number PRJEB51909) and genotyping calling performed in 36 Large White pigs as described in 33 . The total of 13,408,342 positions of Large White variants were compared to the positions of CpGs highlighted in the analyses to filter out common positions using bedtool v2.27.1 34 . The remaining CpGs (non-SNP CpGs), were retained for subsequent analyses, yielding 1,443,424 CpGs (30,026 CpGs removed). The second SNP list was extracted from the Farm animal Genotype-Tissue Expression (FarmGTEx) consortium that produced a comprehensive public resource of genetic regulatory variants across tissues and cells in pigs (PigGTEx 35 ). This list comprised 3,087,268 SNPs documented in three different breeds: Duroc, Landrace, and Large White. It enabled generation of a filtered background of 1,457,149 CpG 10-500 (16,301 CpGs removed). Finally, only CpG 10-500 positions that overlapped between the two filtered backgrounds were conserved (VennDiagram R package) and defined as the final filtered common SNP list consisting of 1,433,115 CpG 10-500 . Descriptive analysis of the genome-wide methylation of blood cell libraries The “methylation_extractor” function from Bismark was used to determine the methylation level at each CpG by counting C and T (methylated and unmethylated before bisulfite conversion, respectively). To monitor the genome-wide methylation pattern before and after SNP filtering, non-supervised hierarchical clustering was conducted using the FactoMineR package 36 on the two matrices of methylation percentages before (3,418,048 CpGs 10-500 ) and after (1,433,115 CpG 10-500 ) filtering. The global mean methylation levels (Suppl. Table 1, Sheet 1) for different sow groups were calculated and statistical analyses were performed using a linear mixed-effect model (lmer function of the R software lme4 package, version v. 4.1.2 (R Core Team, 2023)). The linear model included the replicate as a random factor as well as three fixed factors and their two-by-two interactions: the physiological stage (G98 or L12), the parity class (classified in two groups: low-parity (LP; ranks 2-3) and high-parity (HP, ranks 4 and higher) and the housing environment (E or C). Identification of differentially methylated cytosines/regions (DMCs/ DMRs) Each differential methylation analysis between the two groups was conducted using DSS software (Dispersion Shrinkage for Sequencing data, v2.46.0 37 ) in the default mode. Differentially methylated cytosines (DMCs) were identified considering a minimum methylation difference between groups of 15%, and the adjusted p value ≤ 0.05. This p value adjustment was performed according to the Independent Hypothesis Weighting method, using the alpha parameter set at 5% and the average methylation per group as a covariable. Moreover, DMRs (Differentially Methylated Regions) were identified and defined as regions containing at least three DMCs with a distance between each DMC of 100 bp or less. For each DMC, a particular attention was focused on individual methylation values distribution. It was observed either continue distribution (values varying from 0 to 1) or tri-modal distribution values distributed around of 0 (unmethylated), or around of 0.5 (intermediate-methylated) or around of 1 (methylated). The tri-modal distribution suggests an overlapping SNP/CpG, which can thus be misidentified as unmethylated Cs: indeed, C>T SNPs cannot be distinguished from C>T substitutions that are caused by bisulfite conversion and G>A SNPs induces a loss of methylation on C previous G. Annotation for genomic elements The final SNP-filtered common background, DMCs and DMRs were annotated based on their genomic location using reference files obtained from the Ensembl website ( ftp://ftp.ensembl.org/pub ; release101). They were classified relative to gene features that included exons, introns, transcription start sites (TSS; DMCs located within 100 base pairs up- or downstream of TSS), promoters (DMCs located within −2000 /+100 base pairs of TSS) and transcription termination sites (TTS; DMCs located within 100 base pairs up- or downstream of TTS). A DMC/DMR was considered to belong to a CpG island if there was at least 75% overlap between the site/fragment and the CpG island, and to a shore (up to 2000 base pairs from a CpG island) or shelf (up to 2000 base pairs from a shore). A DMC/DMR was considered as being overlapped by a repetitive element whatever the extent of this overlapping. In order to better characterize the DMRs detected during this study, we reassessed the initially assigned genomic features of a subset of DMRs by intersecting them with functional annotation data (Ensembl Regulation dataset): promoters, enhancers, open chromatin regions, CTCF bindings sites ( https://ftp.ensembl.org/pub/release-113/regulation/sus_scrofa/Sscrofa11.1/annotation/Sus_scrofa.Sscrofa11.1.regulatory_features.v113.gff3.gz ) and, Epigenetically Modified Accessible Regions (EMARs) ( https://ftp.ensembl.org/pub/release-113/regulation/sus_scrofa/Sscrofa11.1/annotation/Sus_scrofa.Sscrofa11.1.EMARs.v113.gff.gz ) Gene Ontology and KEGG enrichment analysis The list of genes co-localized with the resulting DMCs were subjected to enrichment analysis with the DAVID bioinformatic database v2022q2 38 , which included Gene Ontology biological process (GO-BP) terminologies and KEGG pathway categories. The whole list of genes co-localized with SNP-filtered common CpG 10-500 was used as the genome reference. Pathways with at least three genes and p-values <0.05 were considered to be enriched. We also used the gene functional annotations from GeneCards and Ensembl Regulation databases. Statistical analysis of blood cell counts and global DNA methylation level Linear mixed-effect models, were applied for statistical analysis using lmer function of the lme4 package (under R version v. 4.1.2). In these models, the sampling day as fixed effect and sow as random effect were included in order to consider the repeated measures on each sow. Factor effects were considered significant if the p-value was < 0.05. Results Blood leucocytes proportion Blood samples were collected at two time-points (G98 and L12) and blood formula was determined. No significant differences in blood cell sub-population numbers were observed between sows housed in the two environments (p > 0.05). Higher lymphocyte levels were seen in the low-parity (LP) group compared to the high-parity (HP) group of sows (see below, p < 0.05), and significant differences between sampling dates (G98 vs. L12) were observed for monocyte and neutrophil counts (p < 0.001 for both, Table 2 ). From the blood samples, PBMC preparations were performed and used to determine DNA methylation profiles. Although the lymphocytes- and monocytes-count variations are limited, it may to be considered subsequently for methylation data interpretation. View this table: View inline View popup Download powerpoint Table 2. Effect of environmental enrichment blood cell numbers. Values are presented in raw values (mean ± sem). Abbreviations: C, conventional; E, enriched; HP, high-parity; LP, low-parity; H, housing; P, parity; D, dates; G, gestation day; L, lactation day. Description of high inter-individual variability in DNA methylation profiles DNA methylation profiles were initially explored using the initial raw data (3,418,048 CpGs 10-500 before filtering; Suppl. Figure 1) and the SNP-filtered background (1,433,115 CpG 10-500 . ( Figure 2A ) . Both dendrograms displayed quite similar overall structure. Non-supervised hierarchical clustering revealed one cluster with 30 libraries and the other with 28. As expected, for all individual, libraries at G98 and L12 clustered together. Moreover, a distinct cluster of four samples (two physiological stages x two replicates) was observed for each two sows used in both replicates. The SNP-filtering reduced distance scale indicated that the suppression of SNP-overlapping CpGs contributed to a reduction in inter-individual variations. Therefore, the following analyses are only presented with methylation data for the SNP-filtered background. Download figure Open in new tab Figure 2. Hierarchical non-supervised clustering using a matrix obtained after SNP-filtering and including the methylation level of 1,433,115 CpG 10-500 . The library names are coloured according to ( A ) housing environments, ( B ) maternal genealogy ( C ). sow age range and ( D ) sow parity rank. Abbreviations: C, conventional; housing E, enriched housing; DM, different mother; MHS, maternal half-sisters; D, days of age; as described in Table 1 . Black bar underlines the libraries produced from the sows included in both replicates (8 different libraries). Surprisingly, clustering did not separate samples based on housing environment (C vs. E, Figure 2A ). To better understand the factors influencing hierarchical clustering, we first compared the two replicates of the experimental design (R1 vs. R2) and found no discernible clustering pattern (Suppl. Figure 2A), suggesting an absence of significant technical biases or seasonal effects. Sows’ genealogy was not the illustrative factor of clustering When categorized by maternal origin, we identified 19 sows with no half-sisters (38 libraries from G98 and L12) and six sows with at least one half-sister (12 libraries from G98 and L12). Maternal genealogy did not significantly influence the two branches of the dendrogram ( Figure 2B ). We then examined paternal origins, categorizing sows into nine with no half-sisters (18 libraries from G98 and L12) and 16 with at least one half-sister (32 libraries from G98 and L12; Suppl. Figure 2B). Similarly, paternal origins did not significantly influence the clustering, before or after filtering. Age and parity class clearly influenced the individual DNA methylation variability The age of sows at G98 ranged from 494 to 1,099 days. Hierarchical clustering revealed a clearer separation based on age ( Figure 2C ): the first branch (n = 26 libraries) included sows aged 400 to 599 days, while the second branch (n = 32 libraries) included sows aged 600 to 1,099 days. A similar pattern emerged when clustering was analysed by parity class ( Figure 2D ): the first branch consisted of sows with 2-3 parities, while the second branch included sows with 4 or more parities. Only four pairs of libraries did not cluster as expected: E_27 from HP group was clustered with LP group and conversely, C_04, C_05, E_20 from LP group - with HP group. Nevertheless, no clear explanation has been found. In light of these results, a new categorization of the sows based on parity class was made: those in their second or third reproductive cycle (547 ± 73 days old at G98 sampling) formed the low-parity group (LP, n = 15), while those in their fourth or subsequent cycle (860 ± 163 days old at G98 sampling) formed the high-parity group (HP, n = 14). We then investigated the effects of parity (HP vs. LP), physiological state (G98 vs. L12), and housing environment (E vs. C) on global DNA methylation percentages (Suppl. Table 2). The mean global DNA methylation of PBMCs, calculated across all filtered CpGs (1,433,115 CpGs), was 54.93% ± 0.53. No significant differences in global methylation were observed between sows from the E and C environments (p > 0.05, Figure 3A ), although limited, but significant differences were found between parity groups ( Figure 3B ) and physiological states ( Figure 3C ; p < 0.01 for both). Download figure Open in new tab Figure 3. Effects of different factors on DNA global methylation percentages in samples. (A) Comparison of conventional (C) and enriched (E) groups, (B) of low-parity (LP) and high-parity (HP) groups, and (C) gestation day 98 (G98) and lactation day 12 (L12) physiological states. ** P <0.01. Identification of differentially methylated cytosines (DMCs) confirmed the strong influence of parity When comparing HP and LP sows at day G98, we identified 1,358 DMCs, displaying a global reduction in methylation observed in HP sows (1,011 less and 347 more methylated; Suppl. Table 1, Sheet 2). The chromosomal distribution of DMCs ( Figure 4 ) indicated that chromosomes 6, 1, 12, 2, and 3 contained the highest number of DMCs, independently of their size. Chromosome 12, in particular, displayed a high density of DMCs, probably due its higher GC content (47.5% versus a mean of 42.1 ± 2.0% for all chromosomes). Except for chromosomes 10, 18, and X, DMCs were enriched in the terminal regions of most chromosomes. This pattern may partly have resulted from the RRBS method, that targets CpG-rich regions often found at chromosome ends (Suppl. Figure 3A). Download figure Open in new tab Figure 4. Genome-wide distribution of parity related -DMCs, identified at G98. Pig 18 autosomes and X chromosomes are represented by horizontal bars according to their size ( https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000003025.6/ ). DMCs identified between low and high parity on G98 are represented by coloured thin lines: blue lines for less methylated DMCs and red lines for more methylated DMCs in the HP group. The analysis of HP and LP sows at day L12 revealed 680 DMCs (552 less methylated and 128 more methylated in HP sows; Suppl. Table 1, Sheet 3). While the number of DMCs was lower at L12 than at G98, a similar trend of less methylation was observed in HP sows (81.2% of DMCs), with a comparable genomic and chromosomal distribution. However, DMC enrichment at the chromosome termini was less pronounced at L12 (Suppl. Figure 3B). While a small portion of the DMCs were shared between the two time-points (252 DMCs, accounting for 18.5% of G98 DMCs and 37% of L12 DMCs), a similar genomic distribution of DMCs was observed at both stages ( Table 2 ). A significant portion of DMCs was located in intergenic regions (30.3% and 32.3% at G98 and L12, respectively). Among genic regions, introns contained the most DMCs (40.4% and 39.4%), followed by exons (9.4% and 7.5%) and 5’UTRs (8.8% and 8.6%). Most DMCs were located in CpG island (CGi) desert (54.8% and 56.7%), with the remainder distributed across CGi, CGi shores and CGi shelves. The functional analysis revealed that the 698 genes targeted by DMCs at G98 were involved in 15 enriched biological processes (BPs), while 390 genes at L12 were involved in 10 enriched BPs ( Table 3 ). View this table: View inline View popup Download powerpoint Table 3. Genomic distribution of differentially methylated cytosines (DMCs) between high-(HP) and low-parity (LP) sows on gestation day 98 (G98) and lactation day 12 (L12); 1 Numbers and calculated percentages are indicated. Percentages were calculated relative to the total numbers of DMCs at each time point. 2 As a reference, the genomic distribution of 1,433,115 CpG 10-500 (filtered and common CpG 10-500 ) was determined. The box colour shows whether a preferentially genomic distribution of DMCs was observed compared to that of filtered and common CpG 10-500 (green = increase of % and red = decrease of %). View this table: View inline View popup Download powerpoint Table 4. Functional enrichment of DMC-associated genes linked with parity class. Further analysis looked at the distribution of methylation values for the 252 DMCs common to. Both physiological stages. Among them, 51.6% of these shared DMCs exhibited a tri-modal methylation distribution (around of 0 - unmethylated; around of 0.5 - intermediate-methylated; around of 1 - methylated), likely influenced by genetic factors, while 44% displayed a continuous distribution, indicative of a potential epigenetic effect. A small percentage (4.4%) had an unclear pattern (Suppl. data 1.pdf). Differentially methylated regions (DMRs) also highlight the impact of parity rank on DNA methylation profiles At G98, 60 DMRs were identified (58% of them being less methylated and 42% more methylated in HP), spanning 52 unique genes including 9 uncharacterized genes (Suppl. Table 1, Sheet 4 and Suppl. Table 3). Most DMRs were found in intronic (30%) and exonic (16.7%) regions, followed by promoter-associated regions (18.3% in 5’UTRs, 8.3% in promoters, and 1.7% in TSS), 3’UTRs (5%) and ± 10kb of gene (11.6%); only 8.3% of DMRs were localized in intergenic regions. At L12, only 24 DMRs were identified, evenly split between loss and gain of methylation in HP, encompassing 23 unique genes including 7 uncharacterized genes and similarly distributed on gene features than G98-DMRs (Suppl. Table 1, Sheet 4 and Suppl. Table 3). Only a few DMRs displayed a tri-modal methylation distribution (3 at G98 and two at L12 Suppl. data 2.pdf and Suppl. data 3.pdf). In particular, the MYLK DMR including seven DMCs, is in 5’UTR of gene and co-localized with a CGi shelf ( Figure 5A ). A tri-modal distribution of methylation levels was observed: among the LP group, most of them had a methylation level around 50%, while the most of HP group had methylation levels around 100%. Few individuals displayed a methylation level at 0% ( Figure 5B ). The tri-modal pattern was highly correlated between G98 and L12 in both LP ( Figure 5C , R² = 0.84) and HP ( Figure 5D , R² = 0.93) animals. These results illustrated genetics-epigenetics interactions at the MYLK 5’UTR region. Download figure Open in new tab Figure 5. MYLK -DMR. ( A ) Chromosomic position, gene organization and DMR position and sequence ((in red, CG positions identified as DMC). ( B ) Methylation levels of DMCs (named by the start genomic position) for each individual from LP and HP groups at G98; Black arrow indicates TSS. (C) Correlation between methylation levels at G98 and L12 in the LP group. (D) Correlation between methylation levels at G98 and L12 in the HP group. Nevertheless, 55 parities associated DMRs (44 at G98 and 15 at L12 including 4 in common; Suppl. Table 3) exhibited a continuous methylation distribution, suggesting an epigenetic effect of parity class. Although the methylation levels at G98 and L12 were not significantly correlated, the parity-related differential methylation pattern was preserved at both stages for CCND3 ( Figure 6A-D ) , HDAC4 , and SGK1 (Suppl. Figures 4). We re-evaluated the genomic features of the 55 DMRs (Suppl. Table 3) by intersecting them with functional annotation data (Ensembl Regulation dataset). Of these, 15 DMRs (25.86%) co-localized with enhancers, 12 DMRs (20.68%) co-localized with epigenetically modified accessible regions/open chromatin regions and, one co-localized with CTCF binding site; highlighting their potential regulatory role in gene expression control. Together the results from gene feature annotation combined to functional annotation, highlighted the fundamental regulatory role of DNA methylation changes observed with parity rank. More remarkably, DMRs-associated with DNA methylation pattern altered in HP sows encompassed genes which are putatively associated with immune response as schematically presented ( Figure 7 ). Some DMRs presented a loss of methylation in HP sows on G98, encompassing genes such as CD2 (in promoter) , CAMK4 (in 5’UTR) , SECTM1 (in 5’UTR - co-localized with an enhancer), UMODL1 (in upstream region) and CMIP (in intron – co-localized epigenetically modified accessible region). Some DMRs displayed a gain of methylation on G98, such as in intron of SEC14L1 (co-localized with an enhancer), intron of SKI (co-localized epigenetically modified accessible region), 5’UTR of PACS1 and intron of SGK1-DMR (co-localized with an enhancer and CGi-shore). Some DMRs were only observed at L12: they co-localized with intron of genes such as CD5 and TNFRSF1B and displayed a loss of methylation. One of the two SGK1 -DMRs as well as the CCND3- DMR (colocalized in 5’UTR) were found at both time-points and displayed an increase of methylation level in HP sows. Four DMRs were associated with genes involved in epigenetic processes: DNMT3A (more methylated in intron – co-localized with an enhancer) and KDM8 (more methylation in exon) at G98, SIRT2 (less methylated in promoter) at L12 and HDAC4- DMR (less methylated in intron – co-localized with an enhancer) at both time-points. Other DMRs were associated with protein coding genes with functional roles in cell interactions with the extracellular matrix ( EPB41L5, ICAM2, TNNT3, UMODL at G98), cell adhesion ( AMIGO3, THBS2 at L12), signalling cascades (DUSP6, GPR68, GRK4, RAPGEF1 at G98 and RASL10B, GRAP2 at L12), transcription, mRNA splicing and translation control ( SCAF8 , UNK , CELF5 , ZNF205 , POLA2 , FLYWCH2, TCF20, SKI ), and ubiquitin-specific processes ( USP36 , NSMCE1 , KLHL36 , RCHY1 ). Taken together, these results point to that the reproductive life, itself evolving with the age of the individual, induces changes in immune methylation profiles, which could play a central role in controlling the immune response. Download figure Open in new tab Figure 6: CCND3 -DMR. ( A ) Chromosomic position, gene organization and DMR position and sequence (in red, CG positions, * indicates them identified as DMC). ( B ) Methylation levels of DMCs (named by the start genomic position) for each individual from LP and HP groups at G98; Black arrow indicates TSS. (C) Correlation between methylation levels at G98 and L12 in the LP group. (D) Correlation between methylation levels at G98 and L12 in the HP group. Download figure Open in new tab Figure 7. Focus on genes, affected by epigenetic changes induced by parity class and hypothetic impact on immune functions. Only related parity-DMRs colocalized in promoter and/or upstream regulatory region of gene were selected: thus, seven genes were presented. They encode proteins with role in immunity based on information from GeneCards, the Human Gene database ( https://www.genecards.org/ ). The DNA methylation changes in response to housing enrichment depended on sow’s parity Given the significant impact of parity class observed on DNA methylation profiles, we investigated the effects of environmental enrichment (E vs. C) separately for each parity group, first in LP sows (E = 6, C = 9) and then in HP sows (E = 8, C = 6). At G98, we identified 60 DMCs impacted by the housing environment in LP sows (25 less - and 35 more methylated in E) and 42 in HP sows (23 less and 19 more methylated in E). At L12, 35 housing-associated DMCs were seen in LP animals (18 less - and 17 more methylated in E) and 81 in HP animals (46 less - and 35 more methylated in E). No DMCs were shared between the two parity groups whatever their physiological stages, suggesting that the reproductive life of animals may influence the epigenetic response to environmental enrichment (Suppl. Table 1, Sheets 5 and 6). No DMRs were identified as a response to environmental enrichment. As previously, the distribution of methylation profiles for each housing-associated DMC was assessed and categorized as either tri-modal or continuous (Suppl. Table 1, Sheets 5-6; Suppl data 4-7). At G98, whatever the parity group, around 53% of DMCs were classified as continuous suggesting an epigenetic effect. Around 42% were classified as tri-modal suggesting a genetic effect, while 5% exhibited an undetermined status. At L12, only 44 % of DMCs were classified as continuous and 66% tri-modal and/or with undetermined status. The genomic location of the housing-associated DMCs ( Table 5 ) was more often in CpG desert or intergenic regions, and, when found in genic regions, it was mainly situated in intronic areas. Due to the limited number of DMCs identified, pathway enrichment analysis was not feasible. Nevertheless, we were able to classify manually the annotated DMC-associated genes (Suppl. Table 4). For example, some DMCs co-localized with genes involved in immune system and inflammatory response. Two DMCs in TSS-promoter of NOS2 ( Figure 8A ), displayed a decrease of methylation in E-LP-G98 group, one DMC upstream to CHID1 gene was co-localized with an enhancer and exhibited a gain of methylation in E-HP-L12 ( Figure 8G ). Some DMCs colocalized genes involved in cell cycle control and cell adhesion/migration: one DMC in intron of MALD1L1 displayed a loss of methylation in E_LP_G98 ( Figure 8B ) such as one DMC in exon of LINGO 3 was co-localized with a CGi in E-HP-L12 ( Figure 8E ). Interestingly, some DMCs also colocalized with genes involved in neuronal functions: in E-LP-G98 group, two DMCs in intron of SHANK2 gene colocalized with a CpG island and displayed a gain of methylation ( Figure 8C ), while three DMCs in intron of SRGAP2 gene colocalized with a type I transposon/SINE and displayed a loss of methylation with enriched environment ( Figure 8D ). In E-HP-L12, one DMC in intron of SORCS2 co-localized with CGi shelves ( Figure 8F ) and one DMC colocalized in upstream region of SLC1A3 gene ( Figure 8H ) exhibited a loss of methylation with enriched environment. Interestingly, two different DMCs (related to CHILD1 and LINGO 3 ) were colocalized with regulatory elements, suggesting a putative role of methylation changes in transcription control. Download figure Open in new tab Figure 8. Methylation level of housing related DMCs targeting different genes. ( A ) NOS2 – 2 DMCs; ( B ) MALD1L1 – 1 DMC; ( C ) SHANK2 – 2 DMCs; ( D ) SRGAP2 – 3 DMCs; ( E ) LINGO3 – 1 DMC; ( F ) SORCS2 – 1 DMC; ( G ) CHID1 - 1 DMC; ( H ) SLC1A3 – 1 DMC. They were identified either in LP group at G98 (A, B, C and D) or in HP at L12 (E, F, G and H). View this table: View inline View popup Download powerpoint Table 5: Genomic distribution of differentially methylated cytosines (DMCs) in enriched (E) sows relative to conventional (C) sows in low-parity (LP) or high-parity (HP) groups on gestation day 98 (G98) and lactation day 12 (L12). 1 Numbers and calculated percentages are indicated; percentages were based on the total numbers of DMCs at each time point. 2 As a reference, the genomic distribution of 1,433,115 CpG 10-500 (CpG 10-500 background) was determined. The box colour indicates whether a preferentially genomic distribution of DMCs was observed compared to that of CpG 10-500 background (green = increase and red = decrease). Discussion The aim of this study was to investigate how environmental enrichment throughout the reproductive life of pregnant sows affects DNA methylation profiles in immune cells. To do so, we compared PBMCs from sows of different parities and involved in either the gestation or the lactation phase of their reproductive cycle. In our previous research on the same animals, we demonstrated the positive effects of environmental enrichment on cortisol secretion, explorative and social behaviours of the sows, indicating an improved welfare 27 . In the same study, some minor housing effects on gene expression in their PBMCs were reported without significant alteration of blood formula 28 . This additional study has introduced analytical approaches designed to differentiate between true environmental effects and those linked to the reproductive life of animals, while minimizing potential biases arising from individual genetic background. Putative genetic influence on individual DNA methylation profiles Some of our results suggested a genetic influence on inter-individual variations in DNA methylation profiles. In humans, growing evidence has shown that genetic variations play a role in the establishment of DNA methylation marks both in cis (same chromosome) and in trans (across chromosome 39 , 40 ). In our study, genetic diversity could be expected because the experiment was performed at a commercial farm using sows resulting from Landrace x Large White crosses and purchased from various breeder, without animal genotyping. In order to limit genetic influence, data filtering was applied using two lists of known SNPs in the porcine species 33 , 35 . This filtering only concerns the overlapping SNP/CpG. Then, we considered for each DMC (and DMR), the distribution of methylation levels between individuals. However, we identified some DMCs with tri-modal distribution of methylation level, which may be imputable from individual genotype. A marked example of a putative genetic influence on our results was DMR with a gain of methylation identified in the HP group at the 5’UTR position of the MYLK gene, associated with a type I transposon and CGi shelf . MYLK encodes non-muscle myosin light chain kinase, which facilitates the interaction between myosin and actin filaments to produce contractile activity, particularly in neutrophils 41 . In our study, we observed a tri-modal distribution of methylation levels of MYLK DMR for different animals classified either in LP or HP groups. Most LP individuals exhibited intermediate DNA methylation levels, whereas HP individuals predominantly displayed high methylation levels, suggesting an apparent genetic effect on MYLK DMR methylation. In humans, previous research had found a strong association between the methylation levels of specific CpGs and the presence of SNPs in this gene 42 . In addition, Sun et al. 43 demonstrated that the activity of MYLK promoters is influenced by factors such as LPS, TNFα, and hypoxia-inducible transcription factors, and is subject to strong epigenetic regulation, further modulated by proximal SNPs 43 . Our results suggest that genetic-epigenetic interaction at the MYLK promoter level, observed in humans, might be precisely studied in pig. Influence of parity class on individual methylome The hierarchical clustering of the filtered background suggested a significant effect of parity class on individual DNA methylation profiles. Additionally, a significant loss of global methylation in HP sows was observed. According to this data, we identified a large number of DMCs (1,358 at G98 and 680 at L12) between high- and low-parity sows with and, consistently, more DMCs with a lower methylation level in HP sows. A higher proportion of DMCs affected by parity class were located within intronic regions. This does not typically obstruct transcription but can enhance gene expression by affecting regulatory elements such as enhancers and can also affect mRNA processing/alternative splicing 44 , 45 . Interestingly, contrary to the genomic distribution of all background CpG, most parity-associated DMCs were located in CGi deserts, regions typically associated with chromatin status and important for distal regulation of gene 46 . Most DMCs displayed a continue distribution of methylation levels, suggesting that they could be considered as epigenetic biomarkers of parity. Influence of physiological stage on parity-epigenetic effect The reduced number of parity-associated DMCs observed at L12 in this study was consistent with our previous transcriptomic findings, showing a decrease of differentially expressed genes at this physiological stage when compared to the end of gestation 27 . Gestation and lactation are two periods highly hormonally different. Gestation is well known to be associated with maternal physiological adaptations to support foetal development. In particular, maternal leukocytes play a key role in managing inflammation and maintaining immune tolerance towards the foetus. The transition from gestation to lactation involves significant hormonal changes. These hormonal changes may support major evolution in immune competence and response 47 , 48 . We previously reported, using the same blood samples, a significant difference in monocyte and neutrophil numbers between G98 and L12 as well as in lymphocyte numbers in association with parity, without interactions between the two conditions (physiological stage x parity 27 , 28 ). Epigenetic data (herein investigated) such as transcriptomic data 27 were produced using the same PBMC preparations, rich in lymphocytes and monocytes. Taken together, these findings suggested that the substantial metabolic and endocrine reprogramming, as well as weak immune cell proportion variations that occurs during the gestation-lactation transition, may contribute to explain a part of DNA methylation variations. Identification of epigenetic markers associated with parity class We identified numerous DMRs when comparing high- vs. low-parity sows, suggesting a consistent epigenetic effect on the genes concerned. Many DMRs co-localized with regulatory regions suggesting their involvement in gene expression regulation. After excluding DMRs with a tri-modal DNA methylation distribution, four parity DMR-associated genes were linked to the epigenetic machinery: DNMT3A (encoding DNA methyltransferase 3A), KDM8 (encoding a histone demethylase that preferentially recognizes and cleaves monomethylated and demethylated arginine residues on histones H2, H3, and H4) at G98, SIRT2 (encoding an NAD+ (nicotinamide adenine dinucleotide-dependent deacetylase) at L12, and HDAC4 (encoding a class II histone deacetylase) common to both physiological stages. DNMT3A- DMR, located in an intron, exhibited a gain of methylation in the HP group and co-localized with an enhancer. DNMT3A, a crucial enzyme responsible for establishing de novo DNA methylation, plays a pivotal role in maintaining the delicate balance between hematopoietic stem cell differentiation and self-renewal 49 . In humans, DNMT3A function loss mediates changes in the epigenetic landscape and can promote aberrant lymphocyte differentiation and function 49 . In our study, it could be hypothesized that the methylation status of DNMT3A -DMR could induce a reduction of DNMT3A transcription, consistent with the global loss of methylation observed in the HP group. KDM8 -DMR, located in exonic region, exhibited a gain of methylation in the HP group. KDM8 functions as a transcriptional activator by inhibiting HDAC (histone deacetylase) recruitment via the demethylation of H3K36me2. Additionally, HDAC4 -DMR, located in an intronic region, exhibited a loss of methylation in the HP group and co-localized with an enhancer, potentially increasing HDAC4 gene transcription. Together, the opposite methylation patterns of KDM8 -DMR and HDAC4 -DMR in HP suggested a convergence of effects towards an increase in the expression of HDAC4 gene and potentially its activity. Thus, the increase in parity rank could result in deacetylation in favour of chromatin closure. In our sows, age and parity rank were strongly correlated. Thus, in HP sows, global methylation loss, which is a common characteristic of ageing, could be associated to the increase in chromatin instability, while transcriptional inhibition involving HDAC4 activity could be more localized at specific genes. Several parity-DMRs were associated with genes involved in the immune response, including CD2, CAMK4, CMIP, SECTM1, SEC14L1, SGK1, SKI, UMODL1, PACS1, CD5 and TNFRSF1B. Based on the genomic features, functional annotation/regulatory landscapes and DNA methylation status of the DMR, it could be hypothesized that their expression was affected by parity class. Positive transcriptional regulation could affect CAMK4, CD2, CMIP, SECTM1, UMODL1 , and TNFRSF1B in HP sows. CAMK4 belongs to the serine/threonine protein kinase family, has a limited tissue distribution including NK and T-cells. CD2 encodes a surface co-stimulatory receptor present on T-cells. Co-stimulation is needed in addition to antigen-specific signalling through the T cell receptor CD3 complex to achieve full T-cell activation. CD2 expression is increased in memory T-cells, whose numbers increase with ageing 50 . CMIP is a negative regulator of CD28, another T-cell co-stimulatory molecule, whose expression and signalling decreases with ageing in T-cells 51 . It also leads to activation of the c-Maf Th2-specific factor 52 . SECTM1 is primarily expressed in neutrophils and monocytes in peripheral blood. It acts as a ligand for the CD7 receptor, which is also a T-cell co-stimulatory receptor, and synergizes with suboptimal anti-CD28 to increase T cell functions 53 . TNFRSF1B , also known as Tumor necrosis factor (TNF) receptor 2, is expressed, among blood circulating immune cells, only by T-cells. TNFR2 is necessary for antigen-induced differentiation, for the survival of T cells, and is also involved in the regulation of the inflammatory response 54 . In contrast, negative transcriptional regulation could affect the following genes in HP sows: SEC14L1, SGK1, SKI, and PACS1. SEC14L1 encodes a protein that negatively regulates RIG-I 55 , an intra-cytosolic receptor expressed ubiquitously that recognizes RNA viruses and triggers a signalling cascade resulting in the production of type I interferons and pro-inflammatory cytokines 56 . SGK1 is a serine/threonine protein kinase that activates certain potassium, sodium and chloride channels. It is upregulated by cell stress, inflammation as well as a variety of hormones including glucocorticoids 57 . In T-cells, SGK1 acts for differentiation towards the pro-inflammatory Th17 lineage and inhibits development towards the T regulatory phenotype 58 . SKI is a nuclear proto-oncogene that represses TGF-beta signalling, a ubiquitous cytokine acting as a negative growth regulator and playing a major role in the negative regulation of immune responses 59 . PACS1 , a highly conserved cytosolic protein that facilitates the trafficking of cargo between membrane-bound compartments. It is known to be expressed in human peripheral blood leukocytes and is important for maintaining peripheral T and B-cell populations and favours lymphocyte quiescence and survival 60 . Comparisons with the transcriptomes of the same PBMCs used in the present study 27 revealed 28 genes affected at both the epigenetic and transcriptome levels at G98 ( AOC1, ANXA2, BCL3, CHST1, CUX2, DIP2C, DNMT3A, DSP, EPB41L5, EXOC7, GRAMD4, IL9R, JUP, LDLRAD4, MPV17, PRDM11, PRKAR1B, RCHY1, RUNX2, SLC22A15, SLC6A17, SORCS2, SOX13, SYTL3, TATDN3, TNS1, TNFRSF4, and URB1 ). Both transcriptome and methylome data agreed regarding the influence of parity class on biological processes (BP) such as cell-cell adhesion functions (notably with the ANXA2 gene affected in both datasets), the ERK1/ERK2 cascade, and cell surface signalling pathways (particularly those involving G protein-coupled receptors in both datasets). The effects of parity on T-cell differentiation (with the RUNX2 gene affected in both datasets), positive regulation of INF-γ (notably with the BCL3 gene in both datasets) and the positive or negative regulation of transcription also appeared in both datasets, either during gestation or lactation. Previous works in the pig had demonstrated immune differences between primiparous and multiparous animals 61 , 62 , but the effects of successive reproductive cycles in multiparous females were not investigated. Taken together, our results indicate that among multiparous sows, parity class remains an important factor that should be considered when exploring immune-related characteristics. The parity rank was highly correlated with the sows’ age. In humans, age is a well-known factor that induces a loss of DNA methylation across various cell types and tissues 63 . Epigenetic aging can be assessed using various epigenetic clues. Horvath proposed a method to measure the cumulative effect of an epigenetic maintenance system throughout an individual’s life (from various tissues, including PBMCs) 64 . A similar study was conducted in pigs 25 . In our study, we identified 13 genes that were common to both Horvath’s list and our list containing the 1,786 DMC-associated genes affected by parity at both time-points ( DST, FES, GPR68, GRIN2C, LOXL2, MN1, NDUFA13, RAPGEF1, RXRA, SCD5, SLC28A2, ST3GAL4 and UROS ). However, we believe that the differences observed between parity classes could not solely be linked to the aging process. The sows in this study could not be considered as old because their ages ranged from 16 to 36 months, while the natural average lifespan of a pig is approximately 15 years. The possibility of a parity effect (i.e. of accumulated reproductive cycles and not of aging) on DNA methylation is supported by data in humans. Some studies in women have suggested that parity either slows down blood epigenetic aging 65 – 68 or accelerates cell senescence, epigenetic ageing and a shortening of telomeres 69 – 71 . Epigenetic marks of environmental enrichment Despite a clear effect of environmental enrichment on behaviour and salivary cortisol levels in the sows used in this study 27 , 28 , suggesting their improved welfare, the effects of housing environment on DNA methylation were modest. This weak epigenetic effect was in line with our previous transcriptomic observations in PBMC, with few genes affected by the housing environment 27 . Environmental enrichment generated DMCs in genomic regions similar to those seen for parity, i.e. mainly in intronic regions and CpG desert, which might influence chromatin structure and accessibility. Our results agreed with those of a previous study which showed that environmental enrichment increased chromatin accessibility and transcription factor binding in the brains of young mice 72 . In both parity groups, environmental enrichment modulated the methylation of a few DMC-associated genes involved in the emotion functions or disorders resulting from exposure to severe stressors in humans, making these genes relevant candidates for assessing affective states in pigs. In the low-parity group, sows housed in the enriched environment exhibited a loss of methylation in two DMCs in the intronic regions of the SHANK2 gene located inside a CpG island, and in one DMC in the intronic region of the MAD1L1 gene. SHANK2 encodes a synaptic protein. In human, its aberrant expression has been implicated in phenotypes such as intellectual disability and depression 73 . Its role in immune cells is unknown but gain of methylation of the SHANK2 promoter region has been observed in the blood of humans exposed to severe life stressors 74 . In pig, five different SHANK2 transcripts were previously described with variations in coding exons leading to different protein isoforms. Further molecular investigations are needed to be explore the functional epigenetic change induced by environmental enrichment. In humans, specific DNA methylation profiles of MAD1L1 were shown to be linked to reward system functioning, anxiety regulation and associated psychiatric diseases 75 , 76 . In human immune cells, methylation variations of this gene throughout life were also associated with depression 77 , and a higher gene expression was reported in self-reported highly stressed individuals 78 . In the high-parity group, environmental enrichment was associated with a loss of methylation compared to conventional control sows in an exonic CpG in LINGO3 gene, and an intronic CpG in SORCS2 gene. LINGO3 encodes a protein involved in the myelinisation process and has been linked to depression 79 . SORCS2 codes for a receptor of brain-derived neurotrophic factor (BDNF), which is involved in synaptic plasticity 80 and the regulation of mood disorders 81 . In mice, environmental enrichment increased BDNF levels in the brain through epigenetic processes 82 . In pigs, environmental enrichment increased BDNF plasma levels, which have been proposed as an indicator of animal welfare 83 . Three other genes for neurological functioning, regardless of parity class and time point, were identified ( SRGAP2, GRIK4, SLC1A3 ). SRGAP2 , coding SLIT-ROBO Rho GTPase–activating protein 2, is involved in controlling human brain plasticity 84 , but also in cell polarization, regulating neutrophils firm adhesion to endothelial cells 85 , and consequently playing a critical role in inflammatory responses. The loss of methylation observed in an intronic CpG of SRGAP2 may be associated with a downregulation of its expression. GRIK4 encodes a protein that belongs to the glutamate-gated ionic channel family. Glutamate functions as the major excitatory neurotransmitter in the central nervous and glutamate dysfunction is associated mood disorders 86 . SLC1A3, encodes a member of a high affinity glutamate transporter family which plays role in the termination of excitatory neurotransmission. It is also associated with major depressive disorders in humans 87 . Thus, the loss of methylation of a CpG detected herein in the 5’UTR region of GRIK4 gene, as well as in the CpG detected upstream SLC1A3 gene, may be associated with up-regulation of their expression, putatively affecting the levels of the excitatory signal of glutamate in the sows. Collectively, our results highlight some interesting biomarkers of environment effect, occurring in peripheral blood cells, but with putative impact on brain activity and mood regulation. The positive effects of environmental enrichment on brain structure and function have been demonstrated extensively in the literature on rodents 88 . Some studies have also explored the relationship between DNA methylation in brain regions and environmental enrichment. For example, Zocher et al. compared aged and young mice housed in enriched conditions and found that environmental enrichment preserved a methylation pattern in aged mice similar to that seen in younger individuals 89 . Conclusion This study suggested a marked influence of animal genetic background on individual DNA methylation profile, emphasizing the need for a cautious interpretation of the results. Despite this putative genetic influence, a significant effect of parity was observed, even among multiparous sows. This highlights the importance of reproductive life history as a key factor in modulating epigenetic mechanisms related to immune cell function in female animals. Pigs could therefore be considered an effective model for investigating the impact of successive gestation and lactation cycles on maternal DNA methylation profiles and for exploring potential shifts in the epigenetic clock. Furthermore, our findings indicate that positive environmental enrichment exerts a modest effect on DNA methylation profile. This suggests that while environmental enrichment may have beneficial effects on behavioural and physiological aspects, these benefits may be translated into limited molecular changes. Nevertheless, this study provides an original set of biomarkers that need to be further validated and could be routinely followed in further studies with diversified environments in pig farms. Funding details This work was supported by INRAE under the SANBA (Health and Welfare of Farmed Livestock) programme, and the Brittany Region. A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant’s open access conditions. Disclosure statement The authors report there are no competing interests to declare. Ethics of Experimentation The experiment was carried out at the Chambre Régionale d’Agriculture de Bretagne (CRAB) Experimental Farm in Crécom (Saint - Nicolas du Pélem, France), and complied with EU Directive 2010/63/EU on the protection of animals used for scientific purposes. It was approved by the French Regional Ethics Committee on Animal Experiments no 007 and by the French Ministry of Higher Education and Research (APAFIS#27038-2020090315423669 v2). Moreover, this study adheres to ARRIVE guidelines. Data availability The data supporting the findings of this study will be available from the corresponding author upon reasonable request until December 31, 2025. Thereafter, the raw and processed RRBS datasets will be publicly accessible at the European Nucleotide Archive (ENA) under project accession PRJEB90105 (secondary accession ERP173119). Supplementary materials availability All supplementary figures, tables and data sets are upload to a public repository and available using the link “, https://doi.org/10.57745/0GWNEV , Recherche Data Gouv. and referenced under “MERLOT, Elodie; JAMMES, Hélène, 2025, “BEPPI project – data, tables and figures on sow peripheral white blood cell epigenome”. Authors contribution M. M. Lopes: Data curation, formal analysis, Investigation, validation, visualization, writing original draft. G. Monteiro Moreira: Resources, visualization, writing review editing. A. Chaulot. Talmon.: investigation, resources. A. Frambourg.: Resources, software, writing review editing. V. Costes: Resources, software, writing review editing. J. Demars: Resources, writing review editing. E. Merlot: conceptualization, funding acquisition, investigation, project administration, supervision, writing review and editing. H. Jammes: conceptualization, supervision, visualization, writing review and editing Supplemental material Suppl. Figure 1. Hierarchical non-supervised clustering using a matrix obtained before SNP filtering and including the methylation level of 3,418,048 CpGs 10-500 for all libraries (n = 58). The library names are coloured according to ( A ) housing conditions and to ( B ) maternal genealogy (Abbreviations: C, Conventional housing; E, enriched housing; DM, different mother; MHS, maternal half-sisters, as described in Table 1 ). The bar underlines the libraries produced from the sows included in the two replicates (n=8). Suppl. Figure 2. Hierarchical non-supervised hierarchical clustering using a matrix obtained after SNP filtering and including the methylation level of 1,433,115 CpGs 10-500 for all libraries (n = 58). The library names were colored according to (A) assay design and to (B) paternal genealogy. Abbreviations: R, replicate; DF, different father; P-HS, paternal half-sisters - described in Table 1 ). The bar underlines the libraries produced from the sows included in the two replicates (n=8). Suppl. Figure 3. ( A ) Genome-wide distribution of 1,433,115 CpGs 10-500 , analysed by RRBS, SNP filtered and common at all libraries. ( B ) Genome-wide distribution of parity related-DMCs at L12. Pig 18 autosomes and X chromosomes are represented by horizontal bars according to their size ( https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000003025.6/ ). DMCs identified between low and high parity on L12 are represented by coloured thin lines: blue lines for less methylated DMCs and red lines for more methylated DMCs in the HP group. Suppl. Figure 4. For each DMR, the chromosomic position and sequence of DMR and gene organization were represented. The individual methylation rates for each CpG (named by the start genomic position) are shown as a bar chart for each individual from LP (in pink) and HP (in red) groups at G98 and L12. The correlations between methylation levels at G98 and L12 were presented in the LP group (in pink) and in the HP group (in red). (A) SGK1 – DMR1 (located on chromosome 1, 3 DMCs), (B) SGK1 – DMR2 (located on chromosome 1, 5 DMCs) and (C) HDACA – DMR (located on chromosome 15, 3 DMCs). Supplementary Table 1_Sheet 1 : Overview of sequencing parameters for all samples Supplementary Table 1_Sheet 2 : Differentially methylated cytosines (DMCs) affected by parity on gestation day 98 (G98) Supplementary Table 1_Sheet 3 : Differentially methylated cytosines (DMCs) affected by parity on lactation day 12 (L12) Supplementary Table 1_Sheet 4 : Differentially methylated regions (DMRs) affected by parity on gestation day 98 (G98) and lactation day 12 (L12) Supplementary Table 1_Sheet 5 : Differentially methylated cytosines (DMCs) affected by environmental enrichment in low-parity sows on gestation day 98 and lactation day 12 Supplementary Table 1_Sheet 6 : Differentially methylated cytosines (DMCs) affected by housing conditions in high-parity sows on gestation day 98 and lactation day 12 Supplementary Table 2 : Global methylation analysis using SNP_filtered 1 433 115 CpGs Supplementary Table 3 : Functional Annotation of parity related DMRs Supplementary Table 4 : Functional Annotation of housing related DMCs Supplementary data 1 : Methylation distribution_common_parity_DMCs Supplementary data 2 : Methylation distribution_parity_gestation_DMRsall Supplementary data 3 : Methylation distribution_parity_lactation_DMRsall Supplementary data 4 : Methylation distribution_housing_LP_G98_DMCs Supplementary data 5 : Methylation distribution_housing_LP_L12_DMCs Supplementary data 6 : Methylation distribution_housing_HP_G98_DMCs Supplementary data 7 : Methylation distribution_housing_HP_L12_DMCs Acknowledgments The authors are grateful to C Gerard and N Villain, engineers from Brittany’s Regional Agricultural Chamber (CRAB), and C Connan and P Lirzin, staff at the CRAB Experimental Farm, for their help in monitoring the experiment, sampling and data collection. We also thank C Trombani, veterinarian (Breizh Pig, France), for his assistance during blood sampling procedures. We would also like to acknowledge A Vincent, F Thomas, L Le Normand, and P Touanel for their technical assistance (PEGASE, INRAE, Saint-Gilles, France). Our thanks also go to the sows who took part in this study. Funder Information Declared INRAE - SANBA metaprogram (Health and Welfare of Farmed Livestock) Brittany Region - France Footnotes ↵ * Last co-authors https://www.ebi.ac.uk/ena/browser/home https://doi.org/10.57745/0GWNEV References 1. ↵ Mendl , M. , Burman , O. H. P. & Paul , E. S . An integrative and functional framework for the study of animal emotion and mood . Proceedings of the Royal Society B: Biological Sciences 277 , 2895 – 2904 ( 2010 ). OpenUrl CrossRef PubMed 2. Mendl , M. & Paul , E. S. Animal affect and decision-making . Neurosci Biobehav Rev 112 , 144 – 163 ( 2020 ). OpenUrl CrossRef PubMed 3. ↵ Reimert , I. , Webb , L. E. , van Marwijk , M. A. & Bolhuis , J. E . Review: Towards an integrated concept of animal welfare . Animal 17 , ( 2023 ). 4. ↵ Kremer , L. , Klein Holkenborg , S. E. J. , Reimert , I. , Bolhuis , J. E. & Webb , L. E . The nuts and bolts of animal emotion . Neurosci Biobehav Rev 113 , 273 – 286 ( 2020 ). OpenUrl CrossRef PubMed 5. ↵ Bolhuis , J. E. , Schouten , W. G. P. , Schrama , J. W. & Wiegant , V. M . Behavioural development of pigs with different coping characteristics in barren and substrate-enriched housing conditions . Appl Anim Behav Sci 93 , 213 – 228 ( 2005 ). OpenUrl CrossRef Web of Science 6. ↵ Luo , L. , Reimert , I. , Middelkoop , A. , Kemp , B. & Bolhuis , J. E . Effects of Early and Current Environmental Enrichment on Behavior and Growth in Pigs . Front Vet Sci 7 , 525815 ( 2020 ). OpenUrl 7. ↵ Grimberg-Henrici , C. G. E. , Vermaak , P. , Elizabeth Bolhuis , J. , Nordquist , R. E. & van der Staay , F. J . Effects of environmental enrichment on cognitive performance of pigs in a spatial holeboard discrimination task . Anim Cogn 19 , 271 – 283 ( 2016 ). OpenUrl CrossRef PubMed 8. ↵ Luo , L. , Reimert , I. , de Haas , E. N. , Kemp , B. & Bolhuis , J. E . Effects of early and later life environmental enrichment and personality on attention bias in pigs (Sus scrofa domesticus) . Anim Cogn 22 , 959 – 972 ( 2019 ). OpenUrl CrossRef PubMed 9. ↵ Munsterhjelm , C. et al. Environmental enrichment in early life affects cortisol patterns in growing pigs . Animal 4 , 242 – 249 ( 2010 ). OpenUrl CrossRef PubMed 10. ↵ Merlot , E. , Calvar , C. & Prunier , A . Influence of the housing environment during sow gestation on maternal health, and offspring immunity and survival . Anim Prod Sci 57 , 1751 – 1758 ( 2016 ). OpenUrl 11. ↵ Merlot , E. et al. Health and immune traits of Basque and Large White pigs housed in a conventional or enriched environment . Animal 6 , 1290 – 1299 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 12. ↵ Van Dixhoorn , I. D. E. et al. Enriched housing reduces disease susceptibility to co-infection with porcine reproductive and respiratory virus (PRRSV) and actinobacillus pleuropneumoniae (A. Pleuropneumoniae) in young pigs . PLoS One 11 , ( 2016 ). 13. ↵ Cavalli , G. & Heard , E . Advances in epigenetics link genetics to the environment and disease . Nature 571 , 489 – 499 ( 2019 ). OpenUrl CrossRef PubMed 14. ↵ Álvarez-Errico , D. , Vento-Tormo , R. , Sieweke , M. & Ballestar , E . Epigenetic control of myeloid cell differentiation, identity and function . Nat Rev Immunol 15 , 7 – 17 ( 2015 ). OpenUrl CrossRef PubMed 15. ↵ de la Calle-Fabregat , C. , Morante-Palacios , O. & Ballestar , E . Understanding the relevance of DNA methylation changes in immune differentiation and disease . Genes (Basel) 11 , ( 2020 ). 16. ↵ Salameh , Y. , Bejaoui , Y. & El Hajj , N . DNA Methylation Biomarkers in Aging and Age-Related Diseases . Front Genet 11 , ( 2020 ). 17. ↵ Lin , M. W. et al. Comparison of DNA Methylation Changes Between the Gestation Period and the After-Delivery State: A Pilot Study of 10 Women . Front Nutr 9 , ( 2022 ). 18. ↵ Ruiz-Hernandez , A. et al. Environmental chemicals and DNA methylation in adults: A systematic review of the epidemiologic evidence . Clin Epigenetics 7 , ( 2015 ). 19. ↵ Cao-Lei , L. et al. DNA methylation signatures triggered by prenatal maternal stress exposure to a natural disaster: Project ice storm . PLoS One 9 , ( 2014 ). 20. ↵ Corbett , R. J. et al. Assessment of DNA methylation in porcine immune cells reveals novel regulatory elements associated with cell-specific gene expression and immune capacity traits . BMC Genomics 23 , ( 2022 ). 21. ↵ Corbett , R. J. et al. Weaning Induces Stress-Dependent DNA Methylation and Transcriptional Changes in Piglet PBMCs . Front Genet 12 , ( 2021 ). 22. ↵ Wang , H. et al. Genome-wide DNA methylation and transcriptome analyses reveal genes involved in immune responses of pig peripheral blood mononuclear cells to poly I:C . Sci Rep 7 , ( 2017 ). 23. ↵ Yang , Y. et al. Comparative analysis of DNA methylome and transcriptome of skeletal muscle in lean-, obese-, and mini-type pigs . Sci Rep 7 , ( 2017 ). 24. ↵ Jacobsen , M. J. et al. Epigenetic and Transcriptomic Characterization of Pure Adipocyte Fractions From Obese Pigs Identifies Candidate Pathways Controlling Metabolism . Front Genet 10 , ( 2019 ). 25. ↵ Schachtschneider , K. M. et al. Epigenetic clock and DNA methylation analysis of porcine models of aging and obesity . Geroscience 43 , 2467 – 2483 ( 2021 ). OpenUrl CrossRef PubMed 26. ↵ Tatemoto , P. et al. An enriched maternal environment and stereotypies of sows differentially affect the neuro-epigenome of brain regions related to emotionality in their piglets . Epigenetics 18 , ( 2023 ). 27. ↵ Lopes , M. M. et al. Immune transcriptomic profile in adult female pigs: dominance status has more influence than environmental enrichment . BMC Genomics 25 , 1211 ( 2024 ). OpenUrl CrossRef PubMed 28. ↵ Lopes , M. M. et al. Data paper: Dataset describing the effects of environmental enrichment and sows’ characteristics on the peripheral blood mononuclear cell transcriptome . Animal - Open Space 3 , 100078 ( 2024 ). OpenUrl CrossRef 29. ↵ Couret , D. , Otten , W. , Puppe , B. , Prunier , A. & Merlot , E . Behavioural, endocrine and immune responses to repeated social stress in pregnant gilts . Animal 3 , 118 – 127 ( 2009 ). OpenUrl CrossRef PubMed 30. ↵ Gu , H. et al. Preparation of reduced representation bisulfite sequencing libraries for genome-scale DNA methylation profiling . Nat Protoc 6 , 468 – 481 ( 2011 ). OpenUrl CrossRef PubMed 31. ↵ Perrier , J. P. et al. A multi-scale analysis of bull sperm methylome revealed both species peculiarities and conserved tissue-specific features . BMC Genomics 19 , ( 2018 ). 32. ↵ Krueger , F. & Andrews , S. R . Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications . Bioinformatics 27 , 1571 – 1572 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 33. ↵ Boitard , S. et al. Whole-genome sequencing of cryopreserved resources from French Large White pigs at two distinct sampling times reveals strong signatures of convergent and divergent selection between the dam and sire lines . Genetics Selection Evolution 55 , ( 2023 ). 34. ↵ Quinlan , A. R. & Hall , I. M . BEDTools: A flexible suite of utilities for comparing genomic features . Bioinformatics 26 , 841 – 842 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 35. ↵ Teng , J. et al. A compendium of genetic regulatory effects across pig tissues . Nat Genet 56 , 112 – 123 ( 2024 ). OpenUrl CrossRef PubMed 36. ↵ Lê , S. , Josse , J. & Husson , F . FactoMineR: An R Package for Multivariate Analysis . J Stat Softw 25 , 1 – 18 ( 2008 ). OpenUrl CrossRef PubMed 37. ↵ Feng , H. , Conneely , K. N. & Wu , H . A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data . Nucleic Acids Res 42 , ( 2014 ). 38. ↵ Huang , D. W. et al. DAVID Bioinformatics Resources: Expanded annotation database and novel algorithms to better extract biology from large gene lists . Nucleic Acids Res 35 , ( 2007 ). 39. ↵ Hawe , J. S. et al. Genetic variation influencing DNA methylation provides insights into molecular mechanisms regulating genomic function . Nat Genet 54 , 18 – 29 ( 2022 ). OpenUrl CrossRef PubMed 40. ↵ Villicaña , S. et al. Genetic impacts on DNA methylation help elucidate regulatory genomic processes . Genome Biol 24 , ( 2023 ). 41. ↵ Xu , J. et al. Nonmuscle myosin light-chain kinase mediates neutrophil transmigration in sepsis-induced lung inflammation by activating β2 integrins . Nat Immunol 9 , 880 – 886 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 42. ↵ Szilágyi , K. L. et al. Epigenetic contribution of the myosin light chain kinase gene to the risk for acute respiratory distress syndrome . Translational Research 180 , 12 – 21 ( 2017 ). OpenUrl CrossRef PubMed 43. ↵ Sun , X. et al. Genetic and epigenetic regulation of the non-muscle myosin light chain kinase isoform by lung inflammatory factors and mechanical stress . Clin Sci 135 , 963 – 977 ( 2021 ). OpenUrl CrossRef 44. ↵ Borsari , B. et al. Enhancers with tissue-specific activity are enriched in intronic regions . Genome Res 31 , 1325 – 1336 ( 2021 ). OpenUrl Abstract / FREE Full Text 45. ↵ Jones , M. J. , Goodman , S. J. & Kobor , M. S . DNA methylation and healthy human aging . Aging Cell 14 , 924 – 932 ( 2015 ). OpenUrl CrossRef PubMed 46. ↵ Skinner , M. K. & Guerrero-Bosagna , C . Role of CpG deserts in the epigenetic transgenerational inheritance of differential DNA methylation regions . BMC Genomics 15 , ( 2014 ). 47. ↵ Clodi , M. et al. Oxytocin alleviates the neuroendocrine and cytokine response to bacterial endotoxin in healthy men . Am J Physiol Endocrinol Metab 295 , ( 2008 ). 48. ↵ Costanza , M. , Binart , N. , Steinman , L. & Pedotti , R . Prolactin: A versatile regulator of inflammation and autoimmune pathology . Autoimmun Rev 14 , 223 – 230 ( 2015 ). OpenUrl CrossRef PubMed 49. ↵ Belizaire , R. , Wong , W. J. , Robinette , M. L. & Ebert , B. L . Clonal haematopoiesis and dysregulation of the immune system . Nat Rev Immunol 23 , 595 – 610 ( 2023 ). OpenUrl CrossRef PubMed 50. ↵ Linton , P. J. et al. Intrinsic versus environmental influences on T-cell responses in aging . Immunol Rev 205 , 207 – 219 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 51. ↵ Goronzy , J. J. & Weyand , C. M . Mechanisms underlying T cell ageing . Nat Rev Immunol 19 , 573 – 583 ( 2019 ). OpenUrl CrossRef PubMed 52. ↵ Ho , I. C. , Lo , D. & Glimcher , L . H. c-maf Promotes T helper cell type 2 (Th2) and attenuates Th1 differentiation by both interleukin 4-dependent and -independent mechanisms . Journal of Experimental Medicine 188 , 1859 – 1866 ( 1998 ). OpenUrl Abstract / FREE Full Text 53. ↵ Wang , T. et al. K12/SECTM1, an interferon-γ regulated molecule, synergizes with CD28 to costimulate human T cell proliferation . J Leukoc Biol 91 , 449 – 459 ( 2011 ). OpenUrl CrossRef PubMed 54. ↵ Cabal-Hierro , L. & Lazo , P. S . Signal transduction by tumor necrosis factor receptors . Cell Signal 24 , 1297 – 1305 ( 2012 ). OpenUrl CrossRef PubMed 55. ↵ Li , M.-T. et al. Negative Regulation of RIG-I-Mediated Innate Antiviral Signaling by SEC14L1 . J Virol 87 , 10037 – 10046 ( 2013 ). OpenUrl Abstract / FREE Full Text 56. ↵ Thoresen , D. et al. The molecular mechanism of RIG-I activation and signaling . Immunol Rev 304 , 154 – 168 ( 2021 ). OpenUrl CrossRef PubMed 57. ↵ Lauro , D. et al. Role of Serum and Glucocorticoid-Inducible Kinase (SGK)-1 in Senescence: A Novel Molecular Target Against Age-Related Diseases . Curr Med Chem 22 , 3765 – 3788 ( 2015 ). OpenUrl CrossRef PubMed 58. ↵ Wu , C. et al. SGK1 Governs the Reciprocal Development of Th17 and Regulatory T Cells . Cell Rep 22 , 653 – 665 ( 2018 ). OpenUrl CrossRef PubMed 59. ↵ Liu , X. , Sun , Y. , Weinberg , R. A. & Lodish , H. F. Ski/Sno and TGF-β signaling . Cytokine Growth Factor Rev 12 , 1 – 8 ( 2001 ). OpenUrl CrossRef PubMed Web of Science 60. ↵ Nair-Gill , E. , et al. Calcium flux control by Pacs1-Wdr37 promotes lymphocyte quiescence and lymphoproliferative diseases . EMBO J 40 , ( 2021 ). 61. ↵ Choi , S. P. et al. Monitoring mRNA transcription of genes involved in early pregnancy from endometrium and peripheral blood mononuclear cells of pregnant pigs with different parity . Reproduction in Domestic Animals 53 , 1594 – 1599 ( 2018 ). OpenUrl CrossRef PubMed 62. ↵ De Conti , E. R. et al. Effect of sow mass vaccination against Mycoplasma hyopneumoniae on the humoral immune response of newborn piglets . Trop Anim Health Prod 54 , ( 2022 ). 63. ↵ Jones , P. A . Functions of DNA methylation: Islands, start sites, gene bodies and beyond . Nat Rev Genet 13 , 484 – 492 ( 2012 ). OpenUrl CrossRef PubMed 64. ↵ Horvath , S . DNA methylation age of human tissues and cell types . Genome Biol 14 , ( 2013 ). 65. ↵ Russo , J. , Moral , R. , Balogh , G. A. , Mailo , D. & Russo , I. H . The protective role of pregnancy in breast cancer . Breast Cancer Research 7 , 131 – 142 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 66. Barha , C. K. et al. Number of children and telomere length in women: A prospective, longitudinal evaluation . PLoS One 11 , ( 2016 ). 67. Mehta , D. et al. Cumulative influence of parity-related genomic changes in multiple sclerosis . J Neuroimmunol 328 , 38 – 49 ( 2019 ). OpenUrl CrossRef PubMed 68. ↵ Ross , K. M. et al. Immune epigenetic age in pregnancy and 1 year after birth: Associations with weight change . American Journal of Reproductive Immunology 83 , ( 2020 ). 69. ↵ Kresovich , J. K. et al. Reproduction, DNA methylation and biological age . Human Reproduction 34 , 1965 – 1973 ( 2019 ). OpenUrl CrossRef PubMed 70. Ryan , C. P. et al. Reproduction predicts shorter telomeres and epigenetic age acceleration among young adult women . Sci Rep 8 , ( 2018 ). 71. ↵ Ryan , C. P. et al. Pregnancy is linked to faster epigenetic aging in young women . Proc Natl Acad Sci U S A 121 , ( 2024 ). 72. ↵ Espeso-Gil , S. et al. Environmental Enrichment Induces Epigenomic and Genome Organization Changes Relevant for Cognition . Front Mol Neurosci 14 , ( 2021 ). 73. ↵ Ciuculete , D. M. et al. MeQTL and ncRNA functional analyses of 102 GWAS-SNPs associated with depression implicate HACE1 and SHANK2 genes . Clin Epigenetics 12 , 1 – 14 ( 2020 ). OpenUrl CrossRef 74. ↵ Saadatmand , F. , Gurdziel , K. , Jackson , L. , Kwabi-Addo , B. & Ruden , D. M . DNA methylation and exposure to violence among African American young adult males . Brain Behav Immun Health 14 , ( 2021 ). 75. ↵ Trost , S. et al. Investigating the impact of a genome-wide supported bipolar risk variant of MAD1L1 on the human reward system . Neuropsychopharmacology 41 , 2679 – 2687 ( 2016 ). OpenUrl CrossRef PubMed 76. ↵ Snijders , C. et al. Longitudinal epigenome-wide association studies of three male military cohorts reveal multiple CpG sites associated with post-traumatic stress disorder . Clin Epigenetics 12 , ( 2020 ). 77. ↵ Starnawska , A. et al. Epigenome-wide association study of depression symptomatology in elderly monozygotic twins . Transl Psychiatry 9 , ( 2019 ). 78. ↵ Le-Niculescu , H. et al. Towards precision medicine for stress disorders: diagnostic biomarkers and targeted drugs . Mol Psychiatry 25 , 918 – 938 ( 2020 ). OpenUrl CrossRef PubMed 79. ↵ Park , C. et al. Stress, epigenetics and depression: A systematic review . Neurosci Biobehav Rev 102 , 139 – 152 ( 2019 ). OpenUrl CrossRef PubMed 80. ↵ Glerup , S. et al. SorCS2 is required for BDNF-dependent plasticity in the hippocampus . Mol Psychiatry 21 , 1740 – 1751 ( 2016 ). OpenUrl CrossRef PubMed 81. ↵ Jin , Y. , Sun , L. H. , Yang , W. , Cui , R. J. & Xu , S. B . The role of BDNF in the neuroimmune axis regulation of mood disorders . Front Neurol 10 , ( 2019 ). 82. ↵ Costa , G. A. et al. Environmental Enrichment Increased Bdnf Transcripts in the Prefrontal Cortex: Implications for an Epigenetically Controlled Mechanism . Neuroscience 526 , 277 – 289 ( 2023 ). OpenUrl CrossRef PubMed 83. ↵ Rault , J. L. , Lawrence , A. J. & Ralph , C. R . Brain-derived neurotrophic factor in serum as an animal welfare indicator of environmental enrichment in pigs . Domest Anim Endocrinol 65 , 67 – 70 ( 2018 ). OpenUrl CrossRef PubMed 84. ↵ Ho , N. T. T. , Kutzner , A. & Heese , K . Brain plasticity, cognitive functions and neural stem cells: A pivotal role for the brain-specific neural master gene |-SRGAP2-FAM72-| . Biol Chem 399 , 55 – 61 ( 2017 ). OpenUrl CrossRef PubMed 85. ↵ Ren , C. et al. Leukocyte Cytoskeleton Polarization Is Initiated by Plasma Membrane Curvature from Cell Attachment . Dev Cell 49 , 206 – 219 .e7 ( 2019 ). OpenUrl CrossRef PubMed 86. ↵ Knight , H. M. et al. GRIK4/KA1 protein expression in human brain and correlation with bipolar disorder risk variant status . American Journal of Medical Genetics, Part B: Neuropsychiatric Genetics 159 B , 21 – 29 ( 2012 ). OpenUrl CrossRef PubMed 87. ↵ Zhong , X. et al. Identification of Potential Biomarkers for Major Depressive Disorder: Based on Integrated Bioinformatics and Clinical Validation . Mol Neurobiol 61 , ( 2024 ). 88. ↵ Macartney , E. L. , Lagisz , M. & Nakagawa , S . The relative benefits of environmental enrichment on learning and memory are greater when stressed: A meta-analysis of interactions in rodents . Neurosci Biobehav Rev 135 , ( 2022 ). 89. ↵ Zocher , S. , Overall , R. W. , Lesche , M. , Dahl , A. & Kempermann , G . Environmental enrichment preserves a young DNA methylation landscape in the aged mouse hippocampus . Nat Commun 12 , ( 2021 ). View the discussion thread. Back to top Previous Next Posted July 04, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following The effect of environmental enrichment on immune cell DNA methylation profiles depends on the parity of sows Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share The effect of environmental enrichment on immune cell DNA methylation profiles depends on the parity of sows Mariana Mescouto Lopes , Gabriel Costa Monteiro Moreira , Aurélie Chaulot-Talmon , Anne Frambourg , Valentin Costes , Julie Demars , Elodie Merlot , Hélène Jammes bioRxiv 2025.07.01.662510; doi: https://doi.org/10.1101/2025.07.01.662510 Share This Article: Copy Citation Tools The effect of environmental enrichment on immune cell DNA methylation profiles depends on the parity of sows Mariana Mescouto Lopes , Gabriel Costa Monteiro Moreira , Aurélie Chaulot-Talmon , Anne Frambourg , Valentin Costes , Julie Demars , Elodie Merlot , Hélène Jammes bioRxiv 2025.07.01.662510; doi: https://doi.org/10.1101/2025.07.01.662510 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genomics Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41913) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13372) Ecology (19889) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40365) Molecular Biology (17163) Neuroscience (88540) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15130) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9818) Zoology (2269)