Introduction

Celiac disease (CD) is a systemic immune-mediated disorder characterized by the presence of serum autoantibodies. Prospective birth cohort studies with longitudinal screening have revealed the highest incidence between the age of one and four years, although the disease can manifest at any age after gluten introduction1. Development of CD and CD autoimmunity (CDA) requires the presence of HLA-DQ2/DQ8 genotypes2, with additional influence of at least 40 non-HLA loci3,4,5. Although the genetic risk is present in ~ 40% of White people, only a minority of the predisposed individuals are affected, and even genetically identical monozygotic twins show only partial concordance6,7. A rapid change in the incidence in many countries8, together with the varying prevalence between genetically similar populations in nearby geographical regions9, suggest a role of environmental factors in the etiology. Previous evidence indicates that e.g., the amount of gluten intake and exposure to viral infections in early childhood could constitute such factors10,11,12,13,14,15.

Further indication for the role of environment comes from the unexplained fluctuations in CD incidence depending on the season of birth (SoB), most studies finding spring and summer birth being associated with higher risk compared to fall and winter14,16,17,18,19,20. Seasonal variation has also been reported in respect to other autoimmune-mediated diseases such as Crohn’s disease and ulcerative colitis21, juvenile idiopathic arthritis22, Addison’s disease23 and type 1 diabetes (T1D)24. The SoB effect is interesting, as it likely reflects an underlying biological process. For example, Sailani et al. identified distinct patterns of seasonal fluctuations across molecules and biological pathways involved in inflammation and immunity, such as the complement system and acute phase response, as well as observed that the gut microbiome undergoes seasonal changes25. Seasonality in human infectious diseases is well-established and could impact the observed fluctuations. Accordingly, our previous observation that increased CD risk due to gastrointestinal (GI) infections is influenced by SoB and modified by HLA, gluten and breastfeeding14 imply complex gene-environment interactions. SoB-related variation in immune functions such as fluctuations in innate and adaptive cytokine responses and immune cell counts is present already in cord blood and neonatal mucosal immune system, indicating the possibility of additional prenatal effects such as maternal exposures during pregnancy26,27,28. No study has yet investigated interactions between SoB and pregnancy and/or delivery-related factors on the risk of CD/CDA.

Prospective studies focused on early life provide an opportunity to elucidate the gene-environment interactions evoking autoimmunity29. The Environmental Determinants of Diabetes in the Young (TEDDY) is a multicenter birth cohort study following children at risk for T1D and CD. Prospective data collection and sampling enabled us to study interactions between SoB, genetics and environmental factors in the risk of CDA.

Materials and methods

Participants and study design

TEDDY is an observational study following children at-risk for T1D and CD from birth to the age of 15 years6. The aim was to find environmental factors triggering autoimmunity. It includes three study centres in both the USA and Europe.

Inclusion criteria have been described elsewhere31. The reasons for exclusion were unsuccessful re-contacting or failure to schedule the first appointment, lacking HLA data, disease or birth defect that disallowed follow-up or involved a treatment that might interact with the outcomes, child being never tested for CDA, first testing occurring after the age of 36 months and ineligible HLA.

Between 2004–2010, screening for HLA risk were performed on 424,788 infants of whom 21,583 carried the haplogenotypes of interest. Altogether 8,676 families joined the study, and 6,523 children were screened for CDA and followed up until 10 years of age. Of them, 1,262 (19.4%) developed CDA. The median age for seropositivity was 3.3 years and the earliest seroconversion occurred at the age of 10 months.

The study visits occurred quarterly until the age of 4 years and biannually thereafter. Early life data, including maternal and infant illnesses, was obtained from questionnaires given 3 months postpartum. Data about later infections and gluten introduction was obtained from diaries updated by the caregivers. For illnesses, caregivers recorded the nature and date of appearance and the diagnosis.

Serological measurements and outcomes

Blood samples were drawn during every visit. Annual screening for transglutaminase antibodies (tTGA) started from 24 months of age6,32. If a child was seropositive, previous samples were tested retrospectively. CDA was the primary endpoint and was defined as positive tTGA in consecutive samples ≥ 3 months apart. Seropositive children were referred for endoscopy outside the study protocol6,33.

Infections

Maternal infections during pregnancy and child’s infections were categorized as upper/lower respiratory tract infections, GI infections and other34. The child’s infections after the first 3 months were collected at each visit and further translated into ICD-10 codes. The infectious episodes were examined up to age of 12 months.

Statistical analyses

The seasons were defined as spring (March–May), summer (June–August), fall (September–November) and winter (December–February). The association of SoB with risk of CDA were tested using Proportional Hazard models. Children who were negative for CDA were censored on the time of the last tTGA test result or before the age of 10 years. All models were adjusting for HLA haplogenotypes, sex and country. The strength of relations is given by hazard ratios (HR) and 95% confidence intervals (CI). Interactions were evaluated on the ratio scale (multiplicative interaction), or by including a cross-product term in the hazard model, or on the difference scale (additive interaction), by approximating relative excess risk due to interaction (RERI) by adding in the model a 4-category variable describing the occurrence (1) and absence (0) of specific time of birth and genetic factor and estimating the RERI as HR11 – HR10—HR01 + 1 (25). The RERI and 95% CI assess the additional risk due to interaction with RERI > 0 indicating synergistic interaction and RERI < 0 implying antagonistic interaction. The strongest interactions were considered to occur on both the additive and multiplicative scale. P-values < 0.05 were considered significant.

Genetic analyses

SNPs were genotyped using the ImmunoChip and/or the T1DExomeChip (an array with > 90,000 custom content SNPs added to the Infinium® CoreExome-24 v.1.1 BeadChip (Illumina, CA)) at the Center for Public Health Genomics at the University of Virginia. Quality control measures involved ignoring subjects with a call rate < 95% or discordance with reported sex, and SNPs with a call rate < 95%, deviating from Hardy–Weinberg equilibrium in controls with the European ancestry, or with concordance < 99% in duplicates, as described previously35. To investigate interactions outside the HLA region, 87 non-HLA SNPs listed in TEDDY5 and available from ImmunoChip/TEDDY-T1DExomeChip data35 were studied for their association with CDA. The Eagle2 software Version 2.4.1 (https://alkesgroup.broadinstitute.org/Eagle/) (Loh P-R—PMID: 27,694,958) was used to estimates haplotype phase of 4 SNPs (rs840016, rs864537, rs2056626, and rs1723016). Associated SNPs (p-value < 0.05) were selected for interaction testing with spring/summer birth on both additive and multiplicative scale at 5/10 years of age accounting for population stratification (ancestral heterogeneity) (utilizing the top 3 principal components approximated from the T1DExomeChip). Haploview (Version 4.2) was used to generate the linkage disequilibrium (LD) plot.

To study for possible functional role of genes interacting with SoB, mRNA expression from whole-blood samples using RNA-Seq were examined up to age 12 months. The data were available from two TEDDY case–control studies and comprised of 12% of the participants36. No difference in SoB was observed between children with and without expression data (Supplementary Table 1). Differences in expression across CD247 genotypes were tested at 3, 6, 9 and 12 months of age with adjustment for country, sex, HLA, first-degree relative with T1D and presence of islet autoantibodies.

The TEDDY RNA-sequencing (RNA-seq) samples were prepared using Illumina’s TruSeq Stranded mRNA Sample Prep Kit. Sequence were created on the Illumina HiSeq4000 platform using paired-end 2 × 101 bp reads with a targeted 50 million reads per sample by the Broad Institute (Cambridge, MA). The raw sequences were aligned against the Gencode Genome Reference Consortium Build 38 (GRCh38.p12) Release 31 reference genome37 and comprehensive gene annotation. Reads were aligned using the STAR aligner (version 2.6.1d)38 based on the Trans-Omics for Precision Medicine GTEx RNA-Seq analysis pipeline39. Ribosomal RNA (rRNA) was identified and removed using BBSplit version 38.86 (https://sourceforge.net/projects/bbmap/) prior to alignment. Ensembl release 9740 was used as reference.

Gene expression was calculated using RNA-SeQC (version 2.1.0)41 and normalized using the Total Ubiquitous (TU) method from the NormExpression R package (version 0.1.0)42. The ubiquitous genes were selected by the Total Ubiquitous method in the NormExpression R package. The code for this can be found in Sect. 1D of the documentation for the package (https://github.com/cran/NormExpression/blob/master/inst/doc/readme.pdf). For the TEDDY RNA-seq data, it was run using a pre ratio = 1, lower_trim = 0.4, and upper_trim = 0.95. To account for changes in sampling protocol over time, for each gene an RNA Quality Score (RQS) correction was applied to non-zero TU normalized levels by fitting a linear spline with a knot at the 7.5 RQS value. To ensure the RQS adjusted values for samples with RQS >  = 7.5 remained nonnegative and were similarly centered to the original TU normalized data, the original TU normalized values were centered at the 7.5 RQS value before running the linear spline and all adjusted RQS values were shifted to align with the minimum of the original TU normalized value for the samples with RQS ≥ 7.5. Any adjusted RQS values less than zero for samples with RQS < 7.5 were set to zero and were assumed to be too degraded to be considered having any gene expression signal or level. Expression levels were provided as transcripts per million.

Statistical analyses were performed using SAS Software Version 9.4 (https://support.sas.com/software/94/) (SAS Institute, Cary, North Carolina, U.S.A.) and R (Version 4.0.2).

Data availability

Data (https://doi.org/https://doi.org/10.58020/y3jk-x087) reported here can be obtained for request at the NIDDK Central Repository website, Resources for Research (R4R), https://repository.niddk.nih.gov/. The TEDDY RNA Sequencing data that provide support for the findings of the present study have been added in NCBI’s database of Genotypes and Phenotypes (dbGaP) with access code phs001442.v4.p3.

Results

Season of birth and risk of CDA

The number of births across months of the year ranged from 504–586 (Table 1). The 0–9-year incidence of CDA ranged from 20.1–31.9 per 1000 person-years but clustered with high incidence between April and July (Table 1). Spring/summer SoB was associated with an increased likelihood of CDA by 10 years of age (HR 1.29, 95% CI 1.09–1.53, p = 0.004) and less so by 5 years of age (1.18, 1.03–1.34, p = 0.01) after adjusting for HLA, sex or country. No significant interactions by the adjusting factors were observed with SoB (data not shown).

Table 1 Risk of celiac disease autoimmunity (CDA) by 10 years according to season of birth, country of residence, sex and HLA haplogenotypes.
Full size table

Interactions between SoB and non-HLA variants

After accounting for country, sex and HLA, a total of 16 SNPs previously listed in TEDDY4 were associated (p < 0.05, false discovery rate < 0.20) with CDA at 5 and 10 years of age (Supplementary Table 2). Of these, rs864537 in CD247 exhibited a significant interaction with SoB, both on the additive (p = 0.003) and multiplicative scale (p = 0.001), on the risk of CDA. The association of a spring/summer birth with CDA risk was significantly different between children with and without the minor G allele (Fig. 1). A positive correlation was seen only for children having CD247-AA genotype and an excess risk was seen only for children with this genotype and born in spring/summer (RERI additive scale − 0.46, 95% CI − 0.76 to − 0.16 and multiplicative scale 0.69, 0.55–0.87).

Fig. 1
figure 1

Absolute incidence of CDA among children with and without CD247 polymorphism as shown by season of birth (panel A), birth cohort (panel B) and age of CDA development (panel C and D).

Full size image

An examination of the absolute risk of CDA by MoB and year of birth as well as by age of seroconversion showed an excess risk of CDA only for children carrying the CD247-AA genotype born in spring/summer compared to children born in winter/fall or carrying the G allele (Fig. 1A). The spring/summer effect among children with CD247-AA was consistent by year of birth (Fig. 1B) and age of seroconversion (Fig. 1A, C and D) although absolute risk changed by age and birth cohort (data not shown). The risk of CDA for children with CD247-AA was consistently higher for every MoB in spring/summer compared to fall/winter (Supplementary Fig. 1A/B). There were no significant associations between MoB or SoB and CDA among children with CD247-AG (Supplementary Fig. 1C/D) and GG (Supplementary Fig. 1E/F).

Next, the SoB by CD247 within different risk groups identified in TEDDY for CDA was examined. The relative HR of CDA for children born in the spring/summer compared to fall/winter was consistently higher for those with CD247-AA, i.e. a spring/summer birth was associated with an increased risk of CDA in every subgroup except for those who started gluten after 6 months of age (Fig. 2). SoB showed no association with CDA in any of the stratified groups of children who carried the CD247-G allele.

Fig. 2
figure 2

Hazard ratios (± 95% confidence intervals) of spring/summer season of birth on risk of CDA by 10 years of age shown by CD247-G polymorphism and adjusted a for HLA haplogenotypes, country of residence, sex of child and age when child started gluten. a = included in model if not stratified.

Full size image

Other SNPs in CD247

The ImmunoChip contained 35 SNPs in CD247 in addition to rs864537 (Supplementary Fig. 2A). Four of these were excluded from further analysis as the minor allele frequency was less than 2.5%. Of the remaining SNPs, 5 out of 31 (rs840016, rs2056626, rs1723016, rs2995093, rs2982479) showed evidence of interaction with SoB on risk of CDA (Supplementary Fig. 2B), of which 3 (rs840016, rs2056626 and rs1723016) were in linkage disequilibrium (LD) with rs864537. A haplotype analysis (rs840016, rs864537, rs2056626, rs1723016) showed that SoB was associated with CDA risk (HR 1.55, 95% CI 1.27–1.89, p < 0.0001) only if the haplotype comprising of major alleles from all four SNPs was present as homozygous (CATT/CATT, Supplementary Table 3). The interaction of SoB with rs2995093 on CDA risk was no longer significant when adjusting for SoB and rs864357 interaction, whereas the interaction of SoB with rs2982479, further away from rs864537, remained associated. Further analysis revealed that the minor G allele frequency of rs2982479 was low (6.3%), but the allele was associated with a higher risk of CDA (1.51, 1.12–2.04, p = 0.007) for children born in the fall and not in other seasons.

mRNA expression

Data of mRNA expression from whole blood during first year of life were available to analyze CD247 expression in 777 children. The distribution of CD247 variants and SoB was similar for children with and without available expression data (Supplementary Table 1). The expression of CD247 was found to be lower until 6 months for children with the AA genotype in rs864537 compared to those with the G allele (Fig. 3A). Furthermore, haplotype analysis of rs840016, rs864537, rs2056626 and rs1723016 showed that the association was stronger for children with homozygous CATT haplotype (major alleles only) than those with TGGC minor allele only haplotype (Fig. 3B). Differences in expression were most prominent between 3–6 months of age (Fig. 3).

Fig. 3
figure 3

CD247 mRNA gene expression Log2 TPM levels (± 95% confidence intervals) during first year of life as shown by CD247-G polymorphism of SNP-rs864537 (Panel A) and haplogenotypes CATT/CATT or haplogenotype TGGC of SNPs rs840016, rs864537, rs2056626, rs1723016 (Panel B). Differences in means at each time point tested with adjustment for number of time points (p < 0.01) and for country of residence, sex, family history with type 1 diabetes (matching factors in the source Nested Case Control study) as well as HLA and positive for islet autoantibodies (outcome in source Nested Case Control Study).

Full size image

Gene-environment interactions

Febrile respiratory and GI infectious episodes between 3–6 months of age interacted both on additive and multiplicate scale with CD247 on CDA risk by age of 5 years (additive p = 0.008; multiplicative, p = 0.007; Supplementary Table 4) and less so by age 10 years (both additive and multiplicative, p = 0.03). Other infections during pregnancy interacted with CD247-AA genotype on CDA risk by age 5 years (both additive and multiplicative p = 0.02) and 10 years (both p = 0.02), maternal urinary system infections and cold sores explaining these interactions (Supplementary Table 4). Additionally, febrile GI and respiratory infectious episodes in winter/spring between the age of 3–6 months were associated with lower incidence of CDA among CD247-G (Fig. 4), but they alone did not explain the SoB correlation with CDA risk (Fig. 4).

Fig. 4
figure 4

Absolute incidence of CDA from 0 to 4 years of age by CD247 genotype and seasonal infections reported between age 3 and 6 months of age.

Full size image

After adjusting for country of residence, sex, HLA and family history of CD at age 9 months and for non-HLA risk SNPs (Supplementary Table 2), a spring or summer birth remained associated with CDA risk at 5 years of age (HR 1.38, 95% CI 1.13–1.70, p = 0.002) among children with CD247 AA. After additional adjustment for seasonal respiratory or GI infections between 3–6 months of age, the SoB association with CDA risk by 5 years of age was reduced (1.29, 0.97–1.71, p = 0.08). Respiratory or GI infections reported between 3–6 months of age, specifically in August (1.59, 1.08–2.34, p = 0.02) and September (1.58, 1.13–2.21, p = 0.007), remained associated with CDA by 5 years of age among children with CD247-AA. A similar association was seen at 10 years (data not shown).

Discussion

To date, no study has explicitly explained why SoB affects the risk for CD/CDA16,17,19,20,43. Although the differences were relatively small, it was shown that the SoB effect on the risk of CDA is synergistically linked with rs864537 in CD247 as well as with three additional SNPs rs840016, rs2056626 and rs1723016, all in LD. Particularly, children homozygous for the haplotype comprising solely of major alleles of SNPs were at a higher risk of CDA when born during spring/summer (March-August) compared to other seasons/months, whereas this was not observed among those with a haplotype with minor alleles. We also found some variation in the mRNA expression levels of CD247, which could indicate differences in early immune functions depending on distinct genetic features.

CD247 encodes for T-cell surface glycoprotein CD3ζ, a component of the TCR-CD3-complex, which is essential in the adaptive immune response coupling antigen recognition to intracellular signal-transduction pathways30. Aberrant T-cell responses have been linked to CD, and defective expression of CD247 as well as CD3ζ polymorphisms or downregulation may contribute to autoimmunity in T1D, systemic lupus erythematosus (SLE) and thyroidal disease30,44,45,46,47. Furthermore, CD3ζ may play a role in the thymic T-cell selection48. Based on murine models, mutations in CD247 can also impair thymic development of T-cells, revealing the importance of CD3ζ in T-cell development and activation49. Recent findings in vitro suggest that mutations in CD247 can impair proper TCR assembly and lead to defective function of the receptor signaling that might be associated with immunodeficiency and autoimmunity50. Pathways including CD247 have also been linked to antitumor immunity51.

Previous evidence indicates that the mRNA expression of CD247 is lower in children with CD52, further directing towards its potential role in autoimmunity. The expression has also been reported to be reduced in individuals with AA genotype in rs864537 and the minor allele G to be less abundant in celiac patients compared to A allele4. In line, here children carrying the AA genotype had lower levels of mRNA expression than those with AG/GG genotypes. Regarding the SNPs in LD with rs864537, rs840016 polymorphisms have been linked to methylation changes and associated with T1D and dermatomyositis53 and rs2056626 has been associated with systemic sclerosis54,55. Here, a haplotype analysis revealed that SoB was associated with CDA risk only if the haplotype with major alleles existed as homozygous (CATT/CATT), indicating that fine mapping studies in CD247 region are needed to identify the causal variant.

Of note, neonate gene expression is distinct from that of adults and is specific to a given cell type and stimulatory condition. There appears to be an overlap of variants linked to immune-mediated conditions and those which impact on gene expression in neonatal immune cells. Based on Mendelian randomization, changes in neonatal gene expression may have causal effects on autoimmunity risk56. Details of the variations in the expression profiles of CD247 and their possible influence in the evolvement of CD remain to be deciphered, but for example in SLE, defective expression of CD247 and a subsequent T-cell-dysfunction has been suggested to promote autoimmunity due to failure in immune tolerance55.

Maternal infections and childhood febrile respiratory/GI-infections were also associated with CDA depending on the presence of the rs864537 G-allele. Moreover, the incidence of CDA in children reporting febrile infections at 3–6 months of age varied according to CD247 rs864537 genotype. Notably, genotype-dependent differences in the expression of CD247 were most prominent at 3–6 months of age in children with CD247-AA compared to other rs864537 genotypes and, similarly, in subjects with homozygous haplotype comprising of only major alleles for rs840016, rs864537, rs20256626 and rs1723016 (children with CATT/CATT).

The association between childhood and maternal infections and CD is controversial10,11,12,13,14,57,58,59,60,61,62,63,64. Recent evaluation of information gained from prospective studies on genetic and environmental risk of CD points especially to the role of enterovirus and rotavirus infections1. Besides proposed direct and indirect effects of viral or bacterial agents64, season-dependent epigenetic programming could occur during pregnancy. For instance, according to the Developmental Origins of Health and Disease hypothesis, fetal environment has an important role in development of diseases in postnatal life. The programming could happen through epigenetic changes such as DNA-methylation, histone modification and activation or silencing of noncoding RNA-associated genes. SoB has already been associated with DNA-methylation enriched among networks relating to development, the cell cycle and apoptosis65,66. Temporal patterns have also been observed in some placenta-mediated pregnancy complications such as hypertensive disorders during pregnancy67. The role of seasonal fluctuations in immune function is supported by the variable upregulation of genes involved in proinflammatory responses and fluctuating number of immunologically important blood cells across seasons68. Moreover, SoB-related variation in neonatal immunity such as season-dependent overall leukocyte count and cell-subset variation as well as alternating cytokine and chemokine levels has been observed26.

The association between SoB and CDA was present only if gluten was introduced before the age of 7 months. While age of first gluten exposure was not associated with CDA in TEDDY69, high intake during first years of life predisposed to increased risk15. Furthermore, concomitant exposures to gluten and infections together with early cessation of breastfeeding seem to have synergistic effects14,70. The latter might be due to the beneficial effects of breastfeeding on microbiota65, although it was not an independent protective factor in TEDDY71,72. Notably, significant variations in the abundance of particular taxa and microbiome diversity across seasons have been observed, possibly due to changes in the dietary nutrient and fiber composition72,73. The observation that HLA alleles influence the microbiota74 further encourages investigation of gene-environment interactions in autoimmunity.

Strengths and limitations

Our main strength is the prospective design with an exceptionally large birth cohort. The results were also consistent between countries, and the long follow-up time enabled to explore both early and delayed exposures. Inclusion of expression analysis extended the scope beyond genetic associations, although this was possible only for a small subgroup. Self-reporting of early infections and gestational factors is a limitation; however, it might even better cover episodes not requiring healthcare contact. Further limitation is that the study comprised of at-risk children of predominantly white background. An additional factor when considering the generalizability of the results is that the inclusion criteria of the TEDDY study are initially designed to evaluate the risk of T1D, CD and CDA being secondary outcomes. Consequently, the study design does not capture a part of the known HLA-DQ risk alleles for CD. However, early screening for CD autoantibodies is systematically included in the study protocol and thus strengthens the design compared to other prospective studies with T1D as the primary outcome.

Conclusions and potential directions for future research

In conclusion, there was an increased risk of CDA in children born in spring/summer (March–August). This association was present only in children with CD247-AA genotype of the rs864537 variant in CD247 encoding for CD3ζ chain of TCR-CD3-complex and in children homozygous for a haplotype comprised of only major alleles of three additional SNPs rs840016, rs2056626 and rs1723016 in LD. Children with the rs864537 AA genotype and homozygous for the CATT haplotype also expressed reduced levels of CD247 at 3–6 months of age. Information obtained from the present study can contribute to providing a foundation for risk modelling and biomarker discovery that could be exploited in designing future studies aiming to improve our understanding of CD and other autoimmune diseases.