Pet ownership in pregnancy and methylation pattern in cord blood.

Having pets in the house during the first years of life has been shown to protect against allergies. However, the result of different studies is heterogeneous. The aim of this study was to evaluate the methylation pattern in cord blood in relation to pet ownership during pregnancy.We investigated the methylation patterns of 96 cord blood samples, participants of the Epigenetic Hallmark of Maternal Atopy and Diet-ELMA project, born to mothers who either owned pets (n = 32) or did not own pets (n = 64) during their pregnancy. DNA from cord blood was analysed using the Infinium methylation EPIC. For statistical analysis, RnBeads software was applied.We found 113 differentially methylated sites (DMs) in the covariate-adjusted analysis (FDR p < 0.05), with small methylation differences. The top DMs were associated with genes: UBA7, THRAP3, GTDC1, PDE8A and SBK2. In the regional analysis, two promoter regions presented with significance: RN7SL621P and RNU6-211P. Cis-regulatory element analysis revealed significant associations with several immune-related pathways, such as regulation of IL18, Toll signalling, IL6 and complement.We conclude that pet exposure during pregnancy causes subtle but significant changes in methylation patterns in cord blood, which are reflected in the biological processes governing both innate and adaptive immune responses.

Background

Immunological status related to allergy could be determined as early as in utero. Environmental exposures at this particular moment in life are important for the development of a properly functioning immune system. Having pets during the first years of life has been shown to be protective against experiencing allergy at an age later in life. The effect seems to be dose-dependent, species non-specific and correlates with the number of animals present. This suggests that the effect results from extensive microbial and endotoxin exposure [1]. Having pets has been shown to be protective against asthma, allergic rhinitis, atopic dermatitis and sensitisation to aeroallergens, from both animals and pollens [2].

However, the results of studies are heterogeneous, with some showing an opposite effect. For example, household cat or dog exposure has been shown to associate with asthma in Chinese children [3]. Similarly, a meta-analysis of 11 birth cohorts did not reveal any significant association between having different kinds of pets and suffering from allergic diseases such as asthma, allergic rhinitis and atopic dermatitis, but showed protection against sensitisation to aeroallergens at age two in children exposed to furry animals at six months of age [4]. In some studies, the “pet effect” seemed to be species-specific. Cat ownership in the first year of life was revealed to be associated with a decrease in cat allergy and hay fever at age 13, while a dog in the household had no effect [5]. Different results for cats and dogs were also reported in adult populations, and the effect depended on coexistence with other risk factors [6]. A negative impact of pet keeping has been demonstrated in a few studies. An inverse correlation was described in a birth cohort where cat exposure corresponded with early onset of atopic dermatitis, with an interaction between filaggrin gene mutation (R501x) and cat exposure [7]. Epigenetic mechanism regarded as the main regulator of gene expression and the mirror of environmental exposure could be the key element in translating the effect of pet exposure early in life on the development of allergic disease.

The “pet effect” could be considered as a reflection of the effect observed in farming. Thus, the hypothesis could be proposed that household pets act as a “mini-farm”. For farming, there is strong evidence linking the methylation of different genes with protection against allergy. The top candidates are the CD14 promoter, ORMDL1, STAT6, RAD50 and IL13 [8]. There are few studies concerning prenatal exposure to pets and allergy outcomes in children. Until now, no effect was shown for either sensitisation or allergic rhinitis [9].

Previously, we have published an analysis regarding maternal atopy as the main exposure affecting offspring DNA methylation at birth [10]. Here we hypothesise that pet exposure in pregnancy influences DNA methylation in cord blood and this mechanism explains the existing pet effect on allergy occurrence.

Material and methods

Study group

We investigated methylation patterns in the cord blood of 96 newborns, which were children of women aged 25–34, living in an urban area, and participating in project Epigenetic Hallmark of Maternal Atopy and Diet (ELMA). Thirty-two of the women reported keeping animals at home while pregnant: 16 of them were keeping a cat, 14 had a dog, 10 had more than one animal, one had birds, and three had furred rodents. ELMA is an ongoing project designed to establish an epigenetic pattern related to maternal atopy and diet. The characteristics of the study group and analysis regarding maternal atopy were presented in our previous publication [10]. All women were recruited in the 3rd trimester, they were not smoking or exposed to ETS and did not report diabetes, hypertension or pre-pregnancy obesity. The study was approved by the Ethical Committee of the Wroclaw Medical University and all participants signed informed consent forms.

Methylome analysis

DNA from 96 cord blood samples was analysed using the Infinium MethylationEPIC Kit (Illumina, San Diego, CA) according to the standard protocol. All the procedures of sampling, DNA isolation, quality control are described in our previous study [10].

The applied Infinium MethylationEPIC array covers over 850,000 methylation sites scattered across the human genome. Apart from sites within known CpG islands, it includes probes identifying CpG sites outside CpG islands, non-CpG methylated sites identified in human stem cells (CHH sites), differentially methylated sites identified in tumour versus normal cells, FANTOM5 enhancers, ENCODE open chromatin and enhancers, DNase hypersensitive sites, and sites located in miRNA promoter regions.

Data analysis

Data analysis was performed RnBeads 2.4.0 software [11], as in our previous study [10] (and included probes filtering (based on missing data, detection p-value, location on sex chromosomes), normalisation (BMIQ method), bath effect removal (surrogate variable analysis) and covariates analysis. The only difference was the covariates used, which included the mother’s atopy, gestational age and sex. All covariates were included in a final analysis in which the refFreeEWAS method [12] was used for beta-value (methylation value) corrections regarding cell-type heterogeneity and the inclusion of covariates. The refFreeEWAS method allows reference-free deconvolution, providing proportions of putative cell types defined by their underlying methylomes and allowing explicit quantitation of the mediation of phenotypic associations with DNA methylation by cell composition effects. Finally, both FDR-corrected p-values from the refFreeEWAS model (statistical significance) and a RnBeads software score were used to independently detect two sets of sites that were differentially methylated (DM) between the groups considered. The RnBeads score which is originally designed and recommended by the software authors combines initial statistical testing (statistical significance) with a priority ranking based on the absolute and relative effect size of the differences between groups, assigning a combined rank score for differential DNA methylation to each analysed CpG site and genomic region [13].

The detailed analysis of gene content was performed for sites (CpGs) that were differentially methylated between groups with refFreeEWAS model FDR < 0.05. Additionally, to make use of the RnBeads rank (score), the top 100 (approximately top 0.01% of all sites) best-ranked sites were analysed, as well as the top 100 differentially methylated regions (DMRs). This cutoff was much more restrictive than the one proposed in the RnBeads software (the auto.select.rank.cut software function, which selected 126,646 sites and 1008 DMRs). DMRs were defined only within promoters as having the greatest potential to regulate genes expression through DNA methylation. Ensembl gene definitions were used for that, defining promoter regions as the 1500 bases upstream (TSS1500) and 500 bases downstream of the transcription start sites of corresponding genes. The genes associated with DM sites and DMRs were analysed using the WEB-based GEne SeT AnaLysis Toolkit [14] to identify enriched biological processes, pathways and human disease phenotypes. Overrepresentation tests were performed with respect to all known genes using an FDR correction for multiple testing [15]. Additional analysis was done using GREAT 4.0.4 software [16], which assigns biological meaning to a set of non-coding genomic regions by analysing the annotations of the nearby genes (with default settings for region association to genes). Thus, it is particularly useful for studying cis functions of sets of non-coding genomic regions (cis-regulatory elements). The resulting genes were tested for overrepresentation using standard settings of GREAT software, according to GO categories as well as human phenotype annotations.

Results

The comparison of the basic characteristics of the study groups is presented in Table 1.

Table 1

Comparison of groups regarding pet keeping.

	Pet-keeping n = 32 (%)	No pets n = 64 (%)	p*
Maternal atopy	11 (34.38)	39 (60.94)	0.01
Maternal asthma	1 (3.13)	6 (9.38)	0.42
Maternal age (years)	29.37 (±2.37)	29.72 (±2.08)	0.63
Healthy diet index	5.91 (±1.67)	5.77 (±1.96)	0.80
BMI (kg/m2)	20.98 (±3.08)	20.99 (±2.41)	0.49
GWG (kg)	13.95 (±4.00)	13.87 (±4.97)	0.89
Education—tertiary	31 (96.86)	61 (95.31)	0.71
Parity—first child	32 (100.00)	55 (85.94)	0.03
Smoking ever before pregnancy	14 (43.75)	24 (37.50)	0.48
Vitamin supplementation	29 (90.63)	56 (87.50)	0.75
Multivitamin supplementation	24 (0.75)	44 (68.75)	0.64
Folic acid separate supplementation	5 (15.63)	8 (12.50)	0.45
Gestational age (weeks)	39.66 (±1.26)	39.55 (±1.21)	0.58
Birth weight (g)	3448.03 (±475)	3444.1 (±546)	0.30
Caesarean section	13 (40.63)	30 (46.88)	0.56
Apgar score—median	10	10	0.70
Sex
Girls	11 (34.38)	26 (40.63)	0.55
Boys	21 (65.62)	38 (59.38)

Continuous data are given as mean ± SD.

*chi2/Fisher or Mann–Whitney U test.

Analysis of global methylation profile

As a final outcome of the filtering procedures, 22,423 probes and none of the samples were removed (96 samples and 827,101 probes were retained). Initial analysis detecting associations between principal components of the methylation profile and covariates showed that global methylation profile was only affected by gestational age and sex. Maternal atopy did not show a significant effect on global methylation profile, but it was associated with subtler methylation changes, as we previously described [10]. All confounding factors were used for β-value correction in the final statistical analysis. Additionally, a correction for cell-type heterogeneity was implemented as proposed by Houseman et al. (2014) [17]. The performed PCA showed a rather low level of global methylation profile variation among the study groups, with a portion of the samples from both groups having clearly distinct methylation profiles. Additional hierarchical clustering based on Euclidean distance confirmed the rather poor agreement of global methylation profile clustering and sample division into study groups, suggesting more discrete methylation changes caused by pet-keeping during pregnancy (Supplementary file 1).

Differential methylation analysis of single CpG sites

Differential methylation analysis of specific CpG sites revealed that, at the genome-wide level (after FDR correction for multiple testing), β-values of 113 sites differed significantly (FDR < 0.05) between the study groups (Table 2, Supplementary file 2). Analysis of the distribution of the sites across the genome showed that most (65%) were located in non-coding genome regions and outside annotated CpG islands (78%). The remaining sites were predominantly located in gene bodies (17%) or TSS1500 (11%). Of the sites located in the annotated CpG islands (n = 25), 13 (52%) were located in island cores, while the remaining sites were found in island shelves and shores. The average absolute delta-β for the differentially methylated (DM) sites was low with a value of 0.0044. Of the DM sites, 56 were hypomethylated in samples obtained from pet owners children and 57 were hypermethylated. The average delta-beta was slightly lower for hypomethylated sites (0.0037) than for hypermethylated sites (0.0052). The distribution of hyper- and hypomethylated DM sites in different functional elements did not differ significantly between groups (Supplementary file 3).

Table 2

Top 30 DM sites with p-adjusted FDR < 0.05.

Probe ID	Gene	Context	Chr	Mean methylation difference	p.adj. FDR	RnBeads combined rank
cg09579833	UBA7	Body	chr3	0.0023	0.008449078	615076
cg23691858	THRAP3	TSS1500	chr1	0.0017	0.012506292	591286
cg21346980	GTDC1	–	chr2	−0.0077	0.012506292	308246
cg26896818	PDE8A	–	chr15	0.0071	0.012506292	364252
cg23774667	SBK2	TSS1500	chr19	−0.0019	0.012506292	641075
cg00601174	–	–	chr5	−0.0080	0.013286039	294928
cg23181172	LSP1	–	chr11	−0.0026	0.018547607	561243
cg11373423	SNCAIP	–	chr5	−0.0039	0.020878939	471246
cg09182208	–	–	chr8	0.0040	0.020878939	507305
cg10261205	ATHL1	Body	chr11	−0.0055	0.020878939	353086
cg21373811	CDC42BPB	Body	chr14	0.0078	0.020878939	339250
cg26953462	FHOD3	TSS1500	chr18	0.0021	0.020878939	530093
cg26764250	–	–	chr8	0.0054	0.020896712	211010
cg27577234	PCGF1	Body	chr2	−0.0235	0.021600796	92777
cg25102065	GPRC6A	–	chr6	−0.0028	0.021686754	527512
cg16889669	–	–	chr6	0.0119	0.023289415	254307
cg04858110	SCRN1	–	chr7	−0.0115	0.023289415	204253
cg26185731	AMPH	–	chr7	0.0019	0.023289415	555903
cg17001135	AMBRA1	–	chr11	−0.0077	0.023289415	304030
cg04629501	VEGFA	–	chr6	0.0026	0.024046468	460838
cg17939295	–	–	chr13	0.0047	0.024046468	259422
cg19981515	MIR99AHG	–	chr21	−0.0024	0.024046468	605108
cg07398561	AP4B1	–	chr1	0.0027	0.024539895	445446
cg11315500	TRIM25	TSS200	chr17	0.0026	0.024539895	467521
cg02637222	C20orf196	–	chr20	0.0009	0.024539895	717852
cg05295810	CASP6	–	chr4	0.0038	0.024792058	326477
cg02559773	SLC7A11	Body	chr4	0.0044	0.024792058	491416
cg11447849	–	–	chr15	0.0056	0.02497826	425725
cg10040521	CLCA4	–	chr1	0.0037	0.025060386	333449
cg26907429	ALDH9A1	1st Exon	chr1	0.0039	0.02506039	322528

Additional analysis was dedicated to the most differentially methylated sites according to RnBeads software rank (score) (Table 3, Supplementary file 4). This was due to the fact that none of the top-ranked sites were identified as significant after FDR correction. Analysis of the distribution of top-ranked sites across the genome revealed that most of them were located in non-coding genome regions (63%) and outside CpG islands (72%). Of all sites, 21% were located in gene bodies, 14% in gene promotes and 13% within islands cores. The average absolute difference in methylation between study groups for ranked sites was 0.0427, tenfold higher than for FDR-detected sites. Most of the sites (60%) were hypomethylated in the group of pet owners’ children (delta-beta = 0.042). The hypermethylated sites had a similar difference in methylation to the hypomethylated sites. Statistical analysis of the distribution of hyper- and hypomethylated CpG sites showed that they were differentially distributed (p = 0.009) within CpG islands. Among the hypomethylated sites, there were more sites located outside known islands (83% vs. 55%), while hypermethylation was more common in the islands cores (27% vs. 3%) (Supplementary file 3). The mutual relationship between FDR-detected and rank-detected DM sites can be seen in Fig. 1.

Table 3

Top 30 DM sites according to RnBeads score.

Probe ID	Gene	Context	Chr	Mean methylation difference	p.adj.FDR	diffmeth. p	RnBeads combined Rank
cg06653848	–	–	chr11	−0.1054	0.122031518	0.0001818	1226
cg10681981	–	–	chr5	−0.0743	0.12702147	0.000236689	1535
cg14912045	RNFT2	–	chr12	0.0748	0.121427824	0.000158968	2106
cg01550393	LOC102467224	Body	chr5	−0.0338	0.124513233	0.000202637	3288
cg14089881	ICMT	TSS1500	chr1	−0.0481	0.150292714	0.000616767	3392
cg22256607	ANK1	–	chr8	−0.0818	0.153641308	0.000695921	3745
cg18190829	TRABD2B	Body	chr1	0.0577	0.120037772	0.000141243	3895
cg13575925	LOC144571	TSS200	chr12	−0.0636	0.162086886	0.000844093	4304
cg04344225	LINC00709	–	chr10	−0.0758	0.162407247	0.000851302	4334
cg23719516	–	–	chr2	0.0301	0.168156681	0.001034821	5089
cg17118775	–	–	chr8	0.0352	0.168489099	0.001039735	5103
cg22974952	FAM66B	Body	chr8	−0.0268	0.165638137	0.000919948	5274
cg24730756	LOC100292680	Body	chr12	0.11818	0.170099505	0.001099439	5345
cg21550804	–	–	chr8	0.02663	0.162292635	0.000849158	5372
cg05079227	ADAMTS17	Body	chr15	0.06015	0.170654206	0.001123272	5439
cg27131486	ATP9B	–	chr18	−0.0414	0.144064576	0.000484628	5604
cg01655658	HLA-L	Body	chr6	0.06373	0.172013605	0.001170114	5626
cg25277809	–	–	chr2	0.02790	0.172218332	0.001181436	5674
cg27543291	BTN3A2	–	chr6	−0.0327	0.172882576	0.001212748	5802
cg04034998	ADAM12	–	chr10	0.05475	0.121888543	0.000171399	5867
cg20495738	CACNA1C	–	chr12	0.03925	0.062335675	1.05 × 10−5	6567
cg26007606	–	–	chr4	0.08330	0.18186689	0.001543645	7020
cg19021236	MICAL3	–	chr22	0.03998	0.146775623	0.000557573	7065
cg04155231	LOC144571	TSS1500	chr12	−0.0640	0.182684404	0.001574384	7128
cg03876418	–	–	chr11	0.04886	0.186880497	0.001783942	7894
cg07812827	–	–	chr8	0.02809	0.177558793	0.001353962	7906
cg04888234	FCRLA	TSS1500	chr1	0.03067	0.165645598	0.000921328	7942
cg23795048	LOC144571	TSS1500	chr12	−0.0601	0.187057682	0.00180595	7979
cg04488036	–	–	chr11	−0.0224	0.148181979	0.000576978	8146
cg00570269	C13orf33	Body	chr13	−0.0224	0.174226104	0.001260299	8228

Fig. 1

Volcano plot for differential methylation quantified by various metrics.

Color scale according to combined ranking.The mutual relationship between FDR-detected and rank-detected DM sites.

Differentially methylated regions

Regional analysis was restricted to gene promoters in this study. An FDR-based threshold detected only two significant (FDR < 0.05) genome regions, located on chromosomes 5 and 8 (Supplementary file 5). One of the regions was hypo- and another was hypermethylated. The regions were associated with the putative promoters of RN7SL621P (hypomethylated in pet owners) and RNU6-211P (hypermethylated). The rank-based analysis allowed us to analyse the DMRs, of which 47% were hypermethylated in pet owners’ children. The average methylation difference for those sites was 0.018. The delta-beta for hypomethylated sites was similar (0.016) (Supplementary file 6).

Analysis of genes associated with DM sites and DMRs

The 113 DM sites detected by FDR were annotated to 82 different genes (Supplementary file 2). The genes did not enrich any biological processes after FDR but showed pointwise enrichments for processes associated with positive regulation of B cell-mediated immunity and positive regulation of immunoglobulin-mediated immune response (C17orf99, PTPRC, TP53BP1). Pointwise enrichments were also detected for some KEGG pathways, such as ascorbate and aldarate metabolism, cortisol synthesis and secretion, and Fc gamma R-mediated phagocytosis. Among detected disease phenotypes, there were none with significance after FDR and they were difficult to associate with the expected phenotype (Supplementary file 2). Additional GREAT analysis of cis-regulatory elements did not reveal any significantly enriched GO biological processes or phenotypes, however, it showed pointwise significance for epithelial structure maintenance (MKS1, CROCC, SOX9, MUC4) and negative regulation of ERAD pathway (USP25, UBXN1), amongst others (Supplementary file 2).

Sites detected as DM by RnBeads rank were annotated to 59 different genes (Supplementary file 4). The analysis of their associated biological processes did not show significant enrichment of any categories, however, some enrichment was found in processes such as positive regulation of amyloid precursor protein catabolic process, regulation of amyloid-beta formation (GSAP, SPON1), membrane lipid catabolic process (ENPP7, PNLIPRP2) and regulation of angiogenesis (ADAM12, IL17F, PRKCB, PTPRM) (Supplementary file 4). The associated KEGG pathways also did not show statistical significance after correction for multiple testing, however, point enrichments were found for insulin secretion, aldosterone synthesis and secretion, steroid hormone biosynthesis (CACNA1C, CREB5, PRKCB), cortisol synthesis and secretion (CACNA1C, CREB5), and others (Supplementary file 4). However, the GREAT analysis revealed FDR-corrected significant biological processes such as negative regulation of complement activation, lectin pathway (A2M) (FDR = 1.10 × 10−8), negative regulation of immune response (A2M, IL1RL1, TRIM27, FOXF1, HLA-G, DUSP22) (FDR = 4.92 × 10−5), cellular defense response (KLRG1, IL1RL2, HLA-G) (FDR = 5.32 × 10−5), and primary lung bud formation (RDH10) (FDR = 4.67 × 10−5), amongst others (Supplementary file 4).

Region-based analysis using FDR showed only two associated genes, RN7SL621P, and RNU6-211P. Thus, enrichment analysis was applied to the top DM regions classified by RnBeads rank. These regions were associated with 62 different genes. No significant enrichment of biological processes was found for the genes while using FDR, however, uncorrected p-values were significant for processes associated with regulation of leucocyte tethering or rolling (CCL25, CCR2), regulation of cytokine production (CCR2, HLA-DPA1, MAPK13, MAST2, SMAD3, SULF1) and regulation of leucocyte adhesion to vascular endothelial cell (CCL25, CCR2) (Supplementary file 6). Point significance was also found for KEEG pathways such as Th17 cell differentiation (HLA-DPA1, MAPK13, SMAD3), the intestinal immune network for IgA production (CCL25, HLA-DPA1), inflammatory bowel disease (HLA-DPA1, SMAD3), and Th1/Th2 cell differentiation (HLA-DPA1, MAPK13). Cis-regulatory element analysis by GREAT software showed some other immune-related biological processes associated with the analysed regions, including regulation of interleukin-18 biosynthetic process (FDR = 1.18 × 10−54), negative regulation of Toll signalling pathway (NLRP12) (FDR = 1.18 × 10−51), negative regulation of interleukin-6 biosynthetic process (NLRP12, GHRL, FOXJ1, GHSR, INPP5D) (2.88 × 10−24), and antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-independent (B2M, HLA-B, HLA-F, HLA-G, LNPEP) (FDR = 9.77 × 10−16), amongst others. The genes associated with DMRs (by GREAT analysis) also showed enrichment of human phenotypes connected with elevated C-reactive protein level (NLRP12) (FDR = 2.40 × 10−51), recurrent aphthous stomatitis (IL17RC, ITK, NCF4, NLRP12, NLRP3, ORAI1, PRKDC) (FDR = 1.19 × 10−29), progressive gait ataxia (C10orf2, POLG, SACS, SETX, UBE3A) (FDR = 8.00 × 10−26), and sleep–wake cycle disturbance (APOE, DHCR7, PER2, PSEN2, UBE3A) (FDR = 5.33 × 10−25) (Supplementary file 6).

Discussion

Environmental exposures during pregnancy, among them pet exposure are believed to have an impact on offspring DNA methylation patterns, further affecting the development of a child’s immune system. This is the perspective relates to foetal programming of allergic diseases. There is evidence that animal exposure could protect against allergy, specifically if it occurs early in life. The mechanism of such protection is poorly understood but epigenetics could be the key in that process. We have studied global DNA methylation patterns in cord blood of neonates born to mothers who either owned or did not pets during pregnancy and recognised multiple DMs in relation to this exposure. According to the standard procedure, we have identified 113 DMs and two DMRs with FDR-adjusted p-values < 0.05. Overall, the methylation differences were small. The top DMs corresponded to genes that were revealed to be associated with asthma in other EWAS study, such as GTDC1, SNCAIP, CDC42BPB, AMPH, AMBRA1, MIR99AHG, SLC7A11—revealed as DMs in lung tissues. The probes which were assigned as DMs in our analysis were with a different location than those described in other studies, though annotated to the same genes [18, 19]. Among top DMs, there was also UBA7 which is related to the pathways of the innate immune system [20], and THRAP3 is linked to circadian rhythm genes [21]—component of this molecular clock has been shown to be disrupted in asthma. Other gene from the top 30 —PDE8 is expressed in the smooth muscle of the airway and regulates the response to beta-adrenergic receptors, so it has been proposed as a novel therapeutic target for asthma remodelling [22]. VEGFA gene has been revealed as associated with postnatal lung development. VEGF itself has been implicated in remodelling and asthma development [23]. In the regional analysis, we identified two promoter regions with significant methylation differences (FDR-adjusted). Both are related to pseudogenes with unknown functions. In pathway analysis, the genes corresponding to the DMs and DMRs did not enrich significantly any biological processes. However, a trend was observed for regulation of B cell immunity, immunoglobulin production and epithelial structure (in the analysis of individual sites), and cytokine production, Th17 differentiation, B cell immunity and Th1/Th2 differentiation (in regional analysis).

All of these processes may be related to the early events in immune switching occurring as a consequence of pet exposure. Specifically, Th17 may play a role. Th17 cells have the same precursor as Tregs, which are known for their induction of immune tolerance. While Tregs inhibit inflammation, Th17 cells activate processes involved in the immune response against bacteria and fungi. Segmented filamentous bacteria induce Th17, while Bacteroides and Clostridia induce Tregs [24]. The Th17/Treg balance may be important in the context of pet exposure, which is proposed to be a so-called “mini-farm” environment, with a corresponding diversity of microbial stimuli. The same process has been observed in animal models where specific microbial exposure has been shown to be connected to prenatal priming of the acquired asthmatic resistant phenotype [25]. We observe only a trend here, which requires further evaluation in the future. Also, cis-regulatory element analysis did not reveal any significant association for the DMRs and DMs identified with an FDR-corrected p-value cutoff.

According to the RnBeads scores, the top-ranking sites were CpG annotated to the genes: RNTF2, LOC102467224, ICMT, ANK1, and TRABD2B. RNTF2 is of special interest. This gene has been shown to be associated with innate immune processes and IL-3 responses related to inflammation. Priming of RNTF2 with LPS (lipopolisacharyde), in mice exposed to bacteria (Pseudomonas), results in receptor abundance, causing lung inflammation. Conversely, overexpression of RNTF2 is linked to reductions in these processes [26]. ICMT and ANK1 are proteins with functional connections to cell membranes, making them potential candidates important for the integrity of epithelial barriers. The disruption of this integrity is believed to be the first step in allergen sensitisation. TRABD2B is a gene encoding a metalloproteinase involved in Wnt signalling, which regulates processes including fibrosis and smooth muscle hypertrophy—phenomena observed in asthma [27]. The 100 top-ranking (by RnBead score) DMs did not enrich significantly any biological process. However, some were identified in another EWAS study regarding respiratory allergy and the methylation profile of saliva: cg13575925, cg04155231, cg23795048 and cg14089881(p < 0.05, |Δβ| > 0.2) [28]. The first three map to the gene LOC14457, encoding a long non-coding RNA, and the last maps to ICMT [28].

Cis-regulatory element analysis of the 100 top-ranking DMs (by RnBead score), showed significant enrichment for important biological processes potentially associated with allergy, including negative regulation of complement activation connected to A2M gene, negative regulation of immune response related to genes (A2M, IL1RL1, TRIM27, FOXF1, HLA-G, DUSP22), negative regulation of cellular defence (KLRG1, IL1RL2, HLA-G), and lung formation (RDH10). KLRG1 was the most significant gene in this analysis.

The complement system constitutes a link between innate and adaptive immune responses. Complement components C3a and C5a have been shown to take part in type I hypersensitivity responses [29]. 4,5-fold increase in expression of the A2M protein has been revealed by proteome analysis in the nasal mucosa of patients with allergic rhinitis [30]. Additionally, elevated serum A2M levels have been seen in allergic rhinitis patients, correlating with IL17 level and nasal congestion [31]. For IL1RL1, there are some studies concerning genetic variants and methylation in regard to asthma, however, with rather opposite results [32, 33]. Also, some polymorphisms in this gene seem to be associated with specific types of food allergy in children (peanut and egg/chicken) [34]. Nevertheless, IL33/IL1RL1 is thought to be a possible therapeutic target for asthma treatment [35]. Another gene product, TRIM27, takes part in the negative regulation of mast cells, crucial for the allergic response [36]. HLA-G could also suppress allergic inflammation [37], but soluble particles of HLA-G act as a modulator of allergic responses [38] KLRG1 is a marker of ILC cells, which are important components of the innate immune response. ILC cells presenting KLRG1 predominantly differentiate into ILC2-type cells, which exhibit functional characteristics of the type 2 response (independently of Th2), and are also a cellular source of IL5 and IL13. The role of these cells has been demonstrated in the HDM model of chronic asthma and IL9-derived, allergen-induced inflammation [39].

In the regional analysis, LINC00612, A2M-AS1, SULF1, MIR4766 and RNA55P253 were the top-ranking with regard to the RnBeads score. LINC00612 is a non-coding RNA, taking part in apoptosis, inflammation and oxidative stress, and also in the pulmonary microvascular response to ETS [40]. A2M-AS1 is a long non-coding RNA with a known function in chronic obstructive pulmonary disease (COPD) [41]. SULF1 is associated with lung development and endothelial function linked to inflammatory processes in the lungs [42]. The genes corresponding to the 100 top-ranking DMRs by RnBeads score did not enrich significantly any biological processes; the only trend observed was for the intestinal network for IgA production, and Th17, Th1 and Th2 cell differentiation—processes with possible links to food allergy.

In cis-regulatory element analysis for top-ranking DMRs by RnBeads score, we have found significant enrichment in several biological processes which possibly link the environmental exposure to the switch in immune status. Among these processes were: negative regulation of IL18, negative regulation of Toll signalling (NLRP12), and negative regulation of IL6 (NLRP12, GHRL, FOXJ1, GHBR, INPP4D). IL18 belongs to the IL1 family and its primary function is linked to an increase in the production of IFN gamma and Th1 maturation [43]. However, in some other studies, IL18 seems to cooperate with IL5 to produce eosinophilic inflammation in asthma [44]. The exact function of this cytokine may depend on the presence of other cytokines.

Negative regulation of Toll signalling is of special interest regarding exposure to the “pet environment”. This process seems to be directly connected to microbial stimuli. NLRP12 belongs to the family of NOD-like receptors, functioning as an attenuating factor of inflammation by suppressing inflammatory responses in activated monocytes. The exact role of NOD receptors is not well understood. NOD-1 and NOD2 activate Th1 and Th17 cells with increased production of TNF and IL1. SNPs within the corresponding genes have been shown to be associated with childhood asthma, while SNPs within NLRP12 are related to atopic dermatitis. Toll-like receptors play a dual role. They suppress allergy on activation by microbes, but they can also activate the response to allergens [45, 46]. Their presence and activity may be modulated by “animal” allergens [47]. Both NOD signalling and Toll signalling are part of the innate immune response to microbes. This constitutes the link between pet exposure (and the associated increased exposure to microbes) and the change in immune status.

Another component with regard to cis-regulatory elements analysis—IL6 is a typical inflammatory cytokine, taking part in the immune response to viruses, with enhanced production seen in COVID-19 disease Its role in asthma is not well established, but elevated levels have been seen in severe cases. ORMDL3, which has a known role in asthma pathogenesis, activates the IL6 pathway [48]. Among the genes affecting IL6 is FOXJ1, which has a role related to cilia function in the respiratory tract, and is also associated with primary ciliary dyskinesia (PCD) [49].

The evidence for the possible impact of pet-keeping during pregnancy on the methylation pattern in the newborn is scarce, with only a few studies addressing this problem. In one recent study, constant pet ownership from birth to secondary school was associated with changes in methylation status at cg03565274 and was shown to be protective against the combined phenotype asthma and/or rhinitis. According to eQTM analysis, methylation at this site corresponded to ZYMND10 gene expression. This gene is primarily responsible for PCD, a disease related to the abnormal structure of the cilia and recurrent respiratory tract infections. The connection with asthma or rhinitis from a pathomechanistic point of view is not clear [50]. Similarly, methylation status at CD14 has been reported to be associated with pet-keeping throughout childhood. CD14 is known for its role in the response to the microbiome, with altered expression in children raised on a farm [51].

Atopic parents tend to avoid having pets at home. Thus, reverse causality is sometimes considered as an explanation for why pet-keeping seems to protect against allergy. We addressed this problem by adjusting for maternal atopy in the current analysis.

Our study has some limitations. The sample size is small and we did not replicate the study in another cohort. However, the design of this study was based on a unique group with strictly defined risk factors, thus making finding an available cohort for replication difficult. There is not any similar study performed in 850K platform.

Conclusion

We concluded that pet exposure during pregnancy causes subtle but significant changes in methylation patterns, which are reflected in changes in the biological processes governing both the innate (Toll signalling, ILC, complement) and adaptive (IL18, IL6) immune responses. The functional relationship of this alternation requires further research, however, the switch in immunological status appears to be a plausible link between the pet exposure in pregnancy and allergy outcome in offspring

Supplementary file 1

Supplementary file 2

Supplementary file 3

Supplementary file 4

Supplementary file 5

Supplementary file 6