Genomic imprinting analyses identify maternal effects as a cause of phenotypic variability in type 1 diabetes and rheumatoid arthritis

Imprinted genes, giving rise to parent-of-origin effects (POEs), have been hypothesised to affect type 1 diabetes (T1D) and rheumatoid arthritis (RA). However, maternal effects may also play a role. By using a mixed model that is able to simultaneously consider all kinds of POEs, the importance of POEs for the development of T1D and RA was investigated in a variance components analysis. The analysis was based on Swedish population-scale pedigree data. With P = 0.18 (T1D) and P = 0.26 (RA) imprinting variances were not significant. Explaining up to 19.00% (± 2.00%) and 15.00% (± 6.00%) of the phenotypic variance, the maternal environmental variance was significant for T1D (P = 1.60 × 10−24) and for RA (P = 0.02). For the first time, the existence of maternal genetic effects on RA was indicated, contributing up to 16.00% (± 3.00%) of the total variance. Environmental factors such as the social economic index, the number of offspring, birth year as well as their interactions with sex showed large effects.

theory A unique mixed model, previously used on animal data, was applied to investigate the existence of imprinting [36][37][38][39][40]45 . The advantage this model confers is that it is able to simultaneously consider all kinds of imprinting (i.e. maternal, paternal, full, and partial imprinting) in its analyses 17 , ultimately separating maternal imprinting effects from maternal "non-imprinting" effects (e.g., maternal environmental and maternal genetic effects). This was not possible with previous population-scale imprinting analyses models, for example, that of Engellandt and Tier 46 . Our imprinting model estimates two parental gametic variances and one covariance simultaneously. It is written as: where Y is a vector of the response variable; b is a vector of fixed effects; g s is the vector of random gametic effects under a paternal expression pattern; g d is the vector of random gametic effects under a maternal expression pattern; X, Z s , and Z d are the corresponding incidence matrices; and e is the vector of random residuals. The variance-covariance structure is: where σ 2 s and σ 2 d are the gametic variances and σ sd is the covariance. Matrix G is the gametic relationship matrix reflecting the relationships between the gametes of all individuals in a pedigree. It is therefore twice the size of the number of individuals included in the analysis 47,48 . The symbol ⊗ denotes the Kronecker product. The imprinting effect is defined as the vector of differences (g s − g d ) and the corresponding variance of differences is σ 2 i = σ 2 s + σ 2 d −2σ sd , which represents the imprinting variance. Where no imprinting is observed, σ 2 s = σ 2 d = σ sd and σ 2 i = 0.

Results
Parent-of-origin effects. Type 1 diabetes. Genomic imprinting. Using a REML log-likelihood ratio test (RLRT), the significance of the imprinting variance was tested by comparing the logarithmic value of the restricted maximum likelihood (REML log-likelihood) of the linear imprinting model to the REML log-likelihood outcome of a corresponding linear Mendelian model (equivalent null model that assumes the non-existence of imprinting). At a 5% significance level, the analysis revealed that imprinted genes did not significantly contribute to the total genetic variance in T1D susceptibility in the Swedish population data (P = 0.18).
Maternal effects. Initially, data were analysed using linear models in order to test the significance of the variance components. First, a linear model that ignored maternal effects was applied, i.e. only the genetic effect of the individual was included in the model (Mendelian model 1). This led to a T1D heritability estimate (h 2 ) of 0.19 (± 0.1 × 10 −1 ), i.e. 19% of the phenotypic variation in T1D is due to the variation in genetic factors in the analysed population ( Fig. 1). Adding a maternal environmental effect to the model (Mendelian model 2) revealed significant maternal environmental variance with P = 1.60 × 10 −24 . The relative maternal T1D environmental variance was 0.19 (± 0.2 × 10 −1 ), i.e. 19% of the phenotypic variance in T1D is due to the variation in maternal environmental effects (Fig. 1). Heritability was reduced to 0.10 (± 0.1 × 10 −1 ). Augmentation of Mendelian model 2 by the maternal genetic effect (Mendelian model 3) did not change the REML log-likelihood or variance component ratios ( Fig. 1). More detailed information on the variance component estimates in T1D and REML log-likelihood models is provided in Supplementary Table S1. In addition to the linear models, threshold models were applied to account for the binary nature of the phenotypic traits. However, each of the threshold models could only pick up one variance component, i.e. with the addition of parameters, the same amount of variation was explained by additive genetic effects, then by maternal environmental effects, and then by maternal genetic effects (Table 1).
Rheumatoid arthritis. Genomic imprinting. As maternally derived environmental and genetic effects could not be unambiguously disentangled, the imprinting model was applied in two forms: (a) with only the maternal environmental effect in addition to the two parental gametic effects, and (b) with only the maternal genetic effect in addition to the two parental gametic effects. The RLRT of model version (a) did not indicate significant imprinting variance (P = 0.26). Model version (b) led to a REML log-likelihood of 8,408.70, which was slightly smaller than the REML log-likelihood obtained from the corresponding null model containing a gametic effect and a maternal genetic effect (8,408.88). Because the addition of a parameter to a model should result in an REML log-likelihood value either being equal to or larger than that found here, these results could indicate a flat likelihood surface or numerical inaccuracies.
Maternal effects. Maternal effects were initially ignored (Mendelian model 1), which resulted in an h 2 value of 0.10 (± 0.3 × 10 −1 ; Fig. 1). Following the inclusion of a maternal environmental effect (Mendelian model 2), the h 2 value was reduced to 0.85 × 10 −1 (± 0.3 × 10 −1 ). The corresponding maternal variance component was significant at a 5% significance level (P = 0.02). The relative maternal environmental variance was 0.15 (± 0.6 × 10 −1 ; Fig. 1). While the REML log-likelihood value was not significantly altered upon addition of the maternal genetic effect (P = 0.21; Mendelian model 3), the relative maternal RA genetic variance estimate was 0.14 (± 0.4 × 10 −1 ; Fig. 1), the h 2 estimate dropped to zero and the relative maternal environmental variance was reduced from 0.15 (± 0.6 × 10 −1 ) to 0.49 × 10 −1 (± 0.7 × 10 −1 ). To investigate the importance of maternal genetic effects in more detail, Figure 1. Phenotypic variance of type 1 diabetes (a) and rheumatoid arthritis (b) is partitioned into the residual variance (gray), additive genetic variance (green), maternal environmental variance (blue), and maternal genetic variance (red). The variance components were estimated using models with a gametic effect (g), a maternal environmental effect (c), a maternal genetic effect (m), and a residual (e). Standard errors are indicated by error bars. . As an equal additive genetic variance, and thus the same h 2 , was found using the threshold version of Mendelian model 2, maternally derived environmental factors appeared not to play a role (Table 1) in RA. However, when the maternal genetic effect was added (Mendelian model 3) the h 2 value was 0.65 × 10 −2 (± 0.1), while the maternal environmental variance remained zero and the relative maternal genetic variance was 0.20 × 10 −1 (± 0.1; Table 1).
Environmental and sex effects. Type 1 diabetes. Birth year. Across all models (including linear and threshold models), the effects of the year of birth (ranging from 1944 to 2012) were shown to differ significantly (P < 1.00 × 10 −3 ; Table 2; Supplementary Table S3). Effects increased until the end of the 1950s and started declining slightly at the beginning of the 1960s. The effects increased after 1972 until the mid-1980s, declined again until the mid-1990s, and have been increasing ever since (Fig. 2). With the exception of a strong increase of effects and standard errors in 1992 when applying the threshold models (data not shown), trends observed and effects generated under the threshold models were in accordance with those observed for the linear models.
Social economic index of the mother. When considering the social economic index (SEI), analyses revealed that the effects of the mother's SEI differed significantly for T1D with P values ranging from 2.56 × 10 −12 to 1.13 × 10 −9 across all models ( Table 2; Supplementary Table S3). Although small, the largest effect (0.02; ± 3.00 × 10 −3 ) was found for the intermediate group of non-manual employees (code 4). The lowest effect (− 2.00 × 10 −3 ; ± 5.00 × 10 −3 ) was found for professionals as well as higher civil servants and executives (code 5).
Medical region. To investigate the effect of geographical location on T1D susceptibility, medical regions were used. Using linear and threshold models, effects differed significantly for T1D across medical regions with P values ranging from 2.12 × 10 −156 to 4.60 × 10 −104 (  (Table 2; Supplementary Table S3). Negative effects were observed from 1939 with the lowest point obtained in 1958. Since then, RA susceptibility has increased with a positive effect being observed in 1963; a trend that continued until the end of the 1970s. The trend remained constant for approximately 10 years, with a slight increase noticeable towards the end of the 1980s. Increasing standard errors must however be noted (Fig. 2). While a similar trend was observed for the threshold models (data not shown), effects increased in 1974 and remained constant before decreasing in 1991. Large standard errors were observed for effects after 1973. Table 1. Heritability (h 2 ), relative maternal environmental variance (c 2 ), and relative maternal genetic variance (m 2 ) for type 1 diabetes (T1D) and rheumatoid arthritis (RA) estimated using threshold models with a gametic effect (g), a maternal environmental effect (c), a maternal genetic effect (m), and a residual effect (e). Standard errors are in parentheses. for linear models and P = 0.01 for threshold models ( Table 2; Supplementary Table S3). Effect estimates varied little across models. The largest effect (0.11; ± 0.03) was found for unskilled or semi-skilled workers, while the lowest effect (0.02; ± 8.00 × 10 −3 ) was observed for foremen in industrial production and assistant non-manual employees.
SEI of the mother. Maternal SEIs had a small but significant effect on RA susceptibility in offspring under both the linear and threshold models. P values ranged from 0.01 to 0.02 (Table 2; Supplementary Table S3). Effects varied little across models with similar estimates being calculated. The lowest effect was found for foremen in industrial production and assistant non-manual employees (− 0.03; ± 0.01), followed by skilled manual workers (− 0.01; ± 0.01). Except for the unknown SEI group, the highest effect was found for the group of professionals as well as higher civil servants and executives (4.00 × 10 −3 ; ± 0.01).
Number of offspring. In the dataset, women had an average number of 1.94 children (ranging from zero to 11 children; sd = 1.25), while men had an average number of 2.04 children (ranging from zero to 11 children; sd = 1.37). The number of offspring affected RA development significantly across all models (P < 1.00 × 10 −3 ). An inverse and almost linear relationship between RA susceptibility and the number of children is depicted in Fig. 4. While effects greater than zero were estimated for individuals with zero, one or two children, decreasing effects were observed below zero for individuals with more than two children.
Medical region. Medical regions, serving as the proxy for residential and geographic location, differed significantly in their impact on RA with P values ranging from 0.70 × 10 −2 to 0.01 across all models (Table 2; Supplementary Table S3). Effect sizes varied widely across Sweden and were generally small with large standard errors ( Supplementary Fig. S2).
Single child. We found that being a single child or having siblings made a significant difference regarding RA susceptibility with P values ranging from 0.01 to 0.02 across the models (Table 2; Supplementary Table S3). For the linear models, the mean estimated effect of being a single child was − 0.02 (± 9.00 × 10 −3 ), while an effect of − 0.15 (± 0.06) was observed under the threshold models.  ), and the P values (P) for all fixed effects on type 1 diabetes (T1D) and rheumatoid arthritis (RA), which were sex, birth year, social economic index (SEI), number of offspring (no. offspring), medical region, SEI of the mother (SEI mother ), years under observation (years obs ), and whether an individual was a single child or not (single child). Linear mixed models were used containing a gametic effect (g), a gametic effect as father (g s ), a gametic effect as mother (g d ), a maternal environmental effect (c), a maternal genetic effect (m), and a residual effect (e). www.nature.com/scientificreports/ Sex. The incidence of RA was considerably higher in women than in men with 11,442 female and 4,408 male cases, respectively. Sex effects were significantly different with P values ranging from 6.86 × 10 −149 to 6.61 × 10 −101 across models (Table 2; Supplementary Table S3). Interactions between sex and birth year were significant only under the linear models (P = 0.03; P = 1.00 in the threshold models). No clear trend was visible for interaction effects amongst male patients ( Supplementary Fig. S1). Significant interactions between sex and SEI (average P value of 4.56 × 10 −6 ) as well as between sex and the number of offspring were observed (average P value of 2.55 × 10 −15 ). The latter interaction was not significant using threshold models (average P value of 0.48; Supplementary Table S3).

Discussion
Parent-of-origin effects. Our finding that imprinting did not seem to affect T1D susceptibility supported previous findings by McCarthy et al. (1991), who analysed the importance of imprinting in an epidemiological study of clinical data from the Children's Hospital of Pittsburgh IDDM Registry in Pennsylvania, USA 29 . They rejected the imprinting hypothesis and suggested that other genetic and environmental factors may have caused disease occurrence 29 . In addition, Guo and Tuomilehto (2002) stated that the male preponderance in T1D prevalence, fecundity differences, misclassification of T1D and birth order effects could have led to a higher T1D-prevalence in children of T1D-affected fathers 30 . In contrast, a genome-wide association study of European T1D patients showed that an imprinted T1D-associated locus was located within the maternally expressed MEG3 gene 31 .
With regard to RA our findings are consistent with observations made by Zhou et al. (2007), who found that imprinting is unlikely to affect the susceptibility to RA 34 . In most tissues, the IGF2 (insulin-like growth factor   33 . This cell type forms the synovial intimal lining and contributes to cartilage destruction and synovial inflammation 49 . The authors reported that IGF2-linked "loss of imprinting" was responsible for the increased expression that contributed to the autonomous growth of RA fibroblast-like synoviocytes 33 . These findings demonstrated the effects of partially imprinted loci on RA susceptibility. In our study, the hypothesis that imprinted loci were the major cause influencing T1D and RA susceptibility was rejected. However, where a small number of fully or partially imprinted genes with small or moderate effect sizes does exist, the relative imprinting variance (ratio between the imprinting variance and the total additive genetic variance) is expected to be small. Although the imprinting model considers all kinds of imprinting 17,36,37 , our study may have been underpowered and therefore unable to obtain statistical significance. Power to detect significant imprinting variances depends on both the h 2 value of a trait and relative imprinting variance. In this study, heritability estimates for T1D and RA were comparatively small at 19.00% and 10.00%, respectively. However, h 2 estimates vary widely in literature; T1D estimates range between 40 and 92% 3 while RA estimates range between 12 and 60% [11][12][13] . Estimates also depend on the underlying data or the method used in the analysis. For example, twin studies result in broad sense heritability estimates where epigenetics are not taken into account. This increases the risk of inflated estimates 3,50 . Furthermore, sample size and pedigree information affect the ability to estimate genetic parameters and thus the power to detect significance. The depth of the pedigrees was theoretically sufficient to derive imprinting variances, but coancestry information between the maternal and paternal gametes of individuals with RA phenotypes would be a requirement for imprinting variance components analyses. In this context, the availability of genealogy databases in combination with genotypic data would increase the traceability of coancestries. Examples of such databases include the Utah Population database and a number of reliable Icelandic databases 51,52 . The availability of population-scale family trees 53 would further allow the determination of the parental origin of alleles 54 and generally enable large scale human population studies on epidemiological history 55 . Overall, the reliability of data and its impact on T1D and RA diagnoses must be discussed. According to Ludvigsson et al. (2011) the ratio of correct diagnoses for RA in the Swedish Hospital Discharge Register is 93.5% and 87.1% with and without lymphoma, respectively. In this study, no distinction was made between T1D and T2D from ICD-7 through ICD-9. T1D was therefore defined according the age at first hospitalisation as being not older than 20 years. For T1D and T2D 79% of cases were correctly diagnosed. The Swedish Hospital Discharge Register has provided complete national coverage since 1987 with more than 99% of all somatic hospital discharges currently registered 56 .
The replicability of the results of this study depends on the available data basis. Our results are based on data collected until 2007 (RA) and 2012 (T1D). The constant upgrade of the Swedish Hospital Discharge Register, supported by an increasing digitization of data collections, will further improve the size of good quality data and increase the pedigree depth so that higher power to detect significance can be expected when fitting our methods. Furthermore, while disease incidences in the data may not reflect the population incidences, the extension of data will lead to an accumulation of information on family affiliations as well as on the occurrence of the disorders across generations. This will improve the efficiency of estimating heritability in general and the separation of variance components in particular.
With regard to the impact of the mother, we found highly significant maternally derived environmental effects on T1D susceptibility. This finding is supported by Hirschhorn (2003) who stated that there are convincing data that non-genetic factors, such as environmental factors in early childhood, play a role in T1D susceptibility 7 , while Nisticò et al. (2012) found significant effects of the shared environment on T1D susceptibility in an Italian twin cohort study 57 . Maternal environmental effects are expected to especially affect an individual in utero, perinatally, or during early childhood (familial environment) 58 . Factors such as early exposure to cow milk and cereals or a shortened duration of breastfeeding have been mentioned in this context 58 . In addition, we found the mother's SEI to have a significant effect on T1D susceptibility in her offspring. It should however be considered that the SEI  (2017) concluded that T1D susceptibility is unlikely to be affected by autosomal dominant genes 59 . With regard to RA, the maternal genetic variance was found to be significant when maternal environmental effects were not part of the model. Moreover, results indicated an inflation of both, the additive genetic and the maternal environmental variance by maternal genetic effects if ignored in the models. The threshold model results underline these findings. Maternal genetic effects on RA susceptibility have not been reported before. However, as sex-linked effects such as X-chromosomal and mitochondrial effects were not considered in the analysis, an inflation of the maternal genetic variance cannot be excluded. The role of the X-chromosome for the development of RA has previously been discussed; however, these studies have focused on the impact of skewed X-chromosomal inactivation in the context of sex differences in RA susceptibility 60,61 . In addition, the existence of maternal environmental effects on RA susceptibility cannot be excluded. In utero effects, which include maternal smoking 62 or protective effects of maternal non-inherited HLA-antigens 63 , breast-feeding in perinatal life 64 , or hygiene standards during postnatal development 44 have been reported. This study found that RA development was affected by the mother's SEI. To conclude, while their relative importance could not clearly be quantified, both maternal genetic and environmental effects are indicated in RA susceptibility. Furthermore, while the ability to estimate covariances depends on the population structure 65 , the existence of covariances between additive and maternal genetic effects is nevertheless possible given that the h 2 value was reduced to zero when a maternal genetic effect was added to the various models.
Regardless, our study clearly indicated the significance of maternal effects on the development of T1D and RA. This was important as some of the contributory factors could be modified, possibly leading to the prevention of disease or treatment interventions 66 . Environmental and sex effects. The development of T1D and RA were found to be significantly different over the last century when considering the effects of the year of birth. With regard to T1D, Gale (2002) stated that an increase in the incidence thereof over the second half of the twentieth century within a genetically relatively stable population would imply that environmental factors play a role in its etiology 67 . Hence, birth year effects could be attributed to environmental factors which have changed during the last century. There is a reason to suspect that these factors are linked to adjustments in living conditions which are, among others, affected by the economic state of a country. The Swedish economy has been characterised by a steady acceleration in economic growth with decelerations observed during the 1970s and early 1980s 68 . Increased disposable incomes are usually associated with improved living conditions. That living conditions can have an effect on T1D susceptibility is supported by our finding that the SEI of the mother, and thus the environment she provides for her children, significantly affects their likelihood of being diagnosed with T1D. In earlier studies, the susceptibility to T1D has also been shown to be associated with increased, lifestyle-associated linear growth and obesity 69 . While T1D susceptibility varied between regions, medical and the presumed corresponding residential regions were found to significantly affect the likelihood of being diagnosed with this condition. This finding is supported by Tzaneva et al. (2001), who reported that the onset of T1D is strongly dependent on the area of residence 70 . Nevertheless, as observed in Fig. 3, the varying effects associated with the medical regions and year of birth might also be due to the periodic and regional variation in the construction of the Swedish Hospital Discharge Register and Outpatient Register. While the Swedish Hospital Discharge Register was founded in 1964 in six Swedish counties mainly located in the Uppsala region, its nationwide launch was only in 1986. Notably, the Swedish Hospital Discharge Register is now almost 100% with lower coverage of hospital-based outpatient care (approximately 80%) 56 .
As rheumatoid arthritis usually occurs later in life, it cannot accurately be determined which environmental factors have operated during a lifetime. There is evidence that an immune system can be permanently modified by environmental factors at an early age, with growth, nutrition and infectious exposure already having activated the immune system before disease onset 44 . This was supported by our novel finding that the SEI of the mother significantly affected RA susceptibility in children and that the SEI of the individual affected RA susceptibility later in life. The latter effect is also influenced by environmental factors such as the intake of oral contraceptives 71 . In addition, an inverse relationship between the total number of children and RA susceptibility was found. Pregnancy has been reported to typically ameliorate symptoms of RA 72 and breast feeding is known to decrease RA susceptibility 64 . With regard to men, causal effects remain to be discussed. It should be noted that while the effect of the number of children on RA susceptibility was not investigated separately in women and men, significant interactions between sex and the number of children underline existing sex-associated differences. The difficulty in determining when environmental triggers occurred in individual life times may be reflected in the challenging endeavor of estimating the effects of medical regions on disease susceptibility. No information was available for how long individuals have resided in the analysed regions, potentially explaining the small effects and large standard errors. However, as for T1D, the periodic and regional variation in the data collection might also be a reason for this 56 . Nevertheless, geographical variation in RA incidence has previously been reported. Indications were found for the lower incidence rates in South European countries compared to North American and north European countries 73,74 .
In the context of sex differences in T1D and RA susceptibility, effects of the X-chromosome, sex hormones, sex-specific behavior, and sex-specific differences in body composition and structure have been discussed 32,60,61,75,76 . Our study confirmed that there was a significant sex effect in the development of both disorders. Moreover, especially for T1D, significant interactions between sex and birth year were found. This Scientific RepoRtS | (2020) 10:11562 | https://doi.org/10.1038/s41598-020-68212-x www.nature.com/scientificreports/ raises the question whether environmental effects preferentially interact with either sex. Based on the hypothesis that the prevalence of T1D amongst males increases proportionally in relation to disease incidence when the underlying environmental causes preferentially affect males, Gale and Gillespie (2001) reviewed sex ratios in T1D-incidence in multiple populations. They found that, although disease incidence increased, the sex ratio does not change and the trend towards male cases is specific for some populations 75 . They concluded that environmental effects do not interact with males over females. In our study, the effect of sex on T1D was only visible from 1960 to the mid-1980s.
Linear versus threshold model. One difficulty in this study was that the affection status was measured as a categorical trait with binomial distributions. Therefore, the linear models are not appropriate to analyze this data type. Nevertheless, statistical significance of genetic parameters could only be tested using the REML log-likelihood of linear models via RLRT. When a threshold model is used, the ASReml-package employs an approximate likelihood (penalised quasi-likelihood) that cannot be used to test differences 77 . There are currently no alternatives to the ASReml-package for our specific imprinting analysis as it is the only package that allows setting an appropriate correlation between the two parental gametic effects. This function ensures equivalence between the Mendelian and imprinting models. Equivalence is needed to perform an RLRT. Regarding the utilisation of linear and threshold models, both models generated similar results in uncovering the underlying genetic variation for T1D and RA. High correlations between the estimated genetic values indicated a fairly good fit of both models to the data.

Conclusion
Not only was new knowledge gained on the environmental effects on T1D and RA development, the separate contributions of each POE was able to delineate the genetic and phenotypic variation in T1D and RA susceptibility for the first time. Results supported findings that imprinting was of minor importance, but confirmed the role maternal factors played in the occurrence of both diseases. The prospects of fitting complex genetic variance-covariance structures can be expected to further improve given the size of good quality epidemiological data and pedigree depth. Type 1 diabetes. The data selection for T1D was strict. T1D was not defined as an independent entity until the ICD-10 (International Classification of Diseases) classification system 80 . Therefore, T1D was defined through first hospitalisation until an age of 20 years in addition to the ICD codes (ICD-7 code 250, ICD-8 code 260, ICD-9 code 260 or ICD-10 code E10) to unambiguously delineate T1D from type 2 diabetes 81 . Furthermore, as AIs were recorded from 1964, it was assumed that the health status of individuals born prior to 1944 cannot be unambiguously defined. Therefore, only individuals born after 1943 were assigned a case/control status (Fig. 5). Furthermore, only individuals born in Sweden and with known maternal data were considered in the analyses. Overall 27,255 T1D patients (14,626   where Y is a vector of the response variable; b is a vector of fixed effects; g is the vector of random gametic effects; X and Z g are incidence matrices; and e is the vector of random residuals. They are assumed to be normally distributed with a mean of 0 and variances Gσ 2 g .

Model extensions.
Other POEs such as maternal genetic or maternal environmental effects have been reported to be potential nuisance factors, i.e. erroneously assumed to be maternal effects due to imprinting, if not considered within the model definitions 41 . To be able to distinguish maternal effects due to imprinting from other POEs as well as to investigate the importance of the latter, the following versions of Mendelian model 1 were applied: where c and m are vectors of random maternal environmental and maternal genetic effects, respectively. They are assumed to be normally distributed with a mean of 0 and variances I c σ 2 c and Aσ 2 m . Matrix I c is an identity matrix and A is the additive genetic relationship matrix; the incidence matrices Z c , and Z m relate observations and random effects. The significance of the additional random maternal environmental effect in Mendelian model 2 was tested by comparing the REML log-likelihood values of Mendelian model 1 and Mendelian model 2 using a one-sided RLRT. The RLRT is assumed to be χ 2 -distributed with a degree of freedom (DF) of one. The significance of the maternal genetic effect in Mendelian model 3 was tested by comparing the REML log-likelihood values of Mendelian model 2 and Mendelian model 3 using a one-sided RLRT (χ 2 -distributed with DF = 1).
As described for T1D, the maternal environmental variance turned out to be significant. Therefore, the following imprinting model was used to investigate the significance of the imprinting variance in T1D: This model is equivalent to Mendelian model 2 but assumes different contributions of maternal and paternal gametes to the susceptibility to T1D.  www.nature.com/scientificreports/ With regard to RA, the existence of maternal environmental and maternal genetic effects could not be excluded. Hence, the imprinting model was applied with an additional maternal environmental effect (see above) but also with an additional maternal genetic effect: In general, the significance of the imprinting variance was determined by comparing the REML log-likelihood value of the imprinting model with the REML log-likelihood outcome of the corresponding Mendelian models using a one-sided RLRT. The test statistic was assumed to be asymptotically distributed as a mixture of two χ 2 distributions with DF = 1 and DF = 2 36-40,82 . Calculation of population parameters. For the Mendelian models the direct heritabilities were calculated as h 2 Mendel = σ 2 a /σ 2 p , where σ 2 a = 2σ 2 g and σ 2 p = σ 2 g + σ 2 e in Mendelian model 1, σ 2 p = σ 2 g + σ 2 e + σ 2 c in Mendelian model 2, and σ 2 p = σ 2 g + σ 2 e + σ 2 c + σ 2 m in Mendelian model 3. The latter two expressions of σ 2 p were used to calculate c 2 Mendel ( c 2 Mendel = σ 2 c /σ 2 p ) for Mendelian model 2 and Mendelian model 3, respectively. To calculate m 2 Mendel ( m 2 Mendel = σ 2 m /σ 2 p ) for Mendelian model 3, σ 2 p was defined as σ 2 p = σ 2 g + σ 2 e + σ 2 c + σ 2 m . For the imprinting models the direct heritabilities were calculated as h 2 imp = σ 2 a /σ 2 p , where σ 2 a = σ 2 gs + σ 2 gd and σ 2 p = σ 2 gs + σ 2 gd + σ 2 e + σ 2 c or σ 2 p = σ 2 gs + σ 2 gd + σ 2 e + σ 2 m . The first expression of σ 2 p can be used to calculate c 2 imp as c 2 imp = σ 2 c /σ 2 p . The latter expression of σ 2 p can be used to calculate m 2 imp ( m 2 imp = σ 2 m /σ 2 p ).
Threshold model. As the case/control-status is denoted 1/0, the trait can be assumed binomially distributed. Hence, apart from the linear mixed models, which were needed for hypotheses testing, generalised linear mixed models (threshold models) were applied using a logit link and the pseudo-likelihood approach of Gilmour et al. (2015) 77 . The probability that an observation with index k belongs to class zero is: where -using Mendelian model 2 as an example-the linear predictor is: and x k , z g,k , and z c,k are the kth rows of the aforementioned incidence matrices X, Z g , and Z c and the vectors ß, g and c are defined as described for Mendelian model 2.
Fixed effects. Fixed effects included sex (2 levels; 1 = male, 2 = female), birth year (68 birth years for T1D; 57 birth years for RA), medical region (26 levels; individuals lived in 25 medical regions [26 = unknown]), and the SEI of the mother (6 levels) for which the following codes were used: 1 = unskilled/ semi-skilled workers 2 = skilled workers 3 = foremen in industrial production and assistant non-manual employees 4 = intermediate non-manual employees 5 = employed and self-employed professionals and higher civil servants and executives 6 = unknown Note that coding differed across the various versions of the population database. Therefore, they were adjusted and equalised according to the 2012 version. The time an individual was under observation was considered by including the Legendre polynomials from the years under observation up to order three (linear, quadratic and cubic). Additional effects included in the RA model were the SEI of the individual (6 levels), whether an individual was a single child or had any siblings (2 levels; 1 = single child, 2 = siblings), and the number of offspring (12 levels; range of 0-11 children). Furthermore, interactions between sex and birth years were fitted for both disorders. Interactions between sex and SEIs, as well as sex and the number of offspring were considered for RA.
Note that due to the underrepresentation of individuals in birth year levels, individuals were shifted to the next birth year level so that at least five individuals were obtained for each level.
The significance of all fixed effects was tested using an incremental Wald F test implemented within the statistical software package ASReml version 4.2 77 . Note that ASReml caters for linear dependencies in the model by setting singular effects to zero 83 . Variance-covariance components were also estimated via ASReml. The R-packages "kinship2" version 1.6.4 84 and "pedigree" version 1.4 85 in R version 3.5.2 86 were used to prepare the data.

Data availability
Individual-based data for research purposes is protected by strict confidentiality protections imposed by the Swedish Government and applied by Statistics Sweden, but can be made available for research after an ethical Y = Xb + Z s g s + Z d g d + Z m m + e. π(η k ) = exp (η k )/ 1 + exp (η k ) , η k = x k β + z g,k g + z c,k c