Sex differences in genetic and environmental influences on frailty and its relation to body mass index and education

Sample characteristics

Table 1 shows the descriptive statistics of the sample, among the 42,994 twin individuals the mean age of which was 58.8 (standard deviation 10.7; range 41–99). The median FI, i.e., the proportion of endorsed frailty items out of 44 (items listed in Supplementary Table 1), was 0.108 (interquartile range 0.062–0.176). Women were on average older (mean age 59.2 vs 58.4 years) and had a higher median FI (0.119 vs 0.097) than men. Mean BMI and length of education were 25.0 kg/m² and 10.5 years respectively; both were similar across sex.

Phenotypic, intraclass and cross-twin cross-trait correlations

Age-adjusted correlations across zygosity are presented in Table 2 (unadjusted correlations in Supplementary Table 2). For FI, BMI and education, when correlating the same phenotype across twins in a pair (i.e., intraclass correlation), the MZ correlations were larger than DZ correlations and were less than 1, suggesting both genetic and individual-specific environmental influences on all three traits. We observed significantly greater intraclass correlations of FI in women than in men for same-sex twins, implying quantitative sex differences; nevertheless, intraclass correlation of FI in opposite-sex twins did not seem to be smaller than that in same-sex DZ twins, suggesting absence of qualitative sex differences.

FI had a positive phenotypic correlation with BMI (r = 0.13) and a negative correlation with education (r = -0.09). The higher MZ than DZ cross-twin cross-trait correlation between FI and BMI indicated that part of their covariation could be explained by genetic influences in common. In contrast, cross-twin cross-trait correlations between FI and education were comparable for MZ and DZ twins, indicating environmental factors common to twins in a pair, rather than genetic factors, contributing to their covariance.

Univariate twin modelling

Univariate sex-limitation models of FI were fitted as shown in Table 3, providing estimates of the variance components in men and women: additive genetic factors (A), dominance genetic factors (D), common environmental factors (C), and unique environmental factors (E). Overall, the best-fitting model was an ADE model (i.e., a model including A, D and E variance components) with only quantitative sex differences. Compared to the saturated model (i.e., a model that fully describes the observed data), the full ADE model did not provide a worse model fit, indicating that the model fit the data well, but both the ACE and AE models provided a poor fit to the data. Removing quantitative sex differences (i.e., different magnitude of heritability) from the ADE quantitative sex-limitation model resulted in a reduced model fit. Meanwhile, the genetic correlation for opposite-sex twins (r_fm), which is used for assessing qualitative sex difference (i.e., different sources of genetic variance), was close to, and statistically non-significantly different from unity. This suggests the presence of quantitative, but not qualitative sex differences in the genetic contributions to FI. We observed largely dissimilar proportions of additive and dominance genetic influences on FI for men and women. Alternatively, we estimated the broad-sense heritability, H, of FI (i.e., proportion of variance that can be explained by the total of additive and dominance genetic factors), which was estimated to be 45% [95% confidence interval (CI): 41–48%] in men, and at a statistically significantly higher proportion of 52% (95% CI: 50–55%) in women. The rest of the variation in FI was explained by unique environmental factors.

Bivariate twin modelling

Using bivariate twin models, we assessed the overlap of genetic and environmental variance components of FI with BMI and education and allowed for sex differences. The best-fitting bivariate model for FI and BMI was an ADE model (Figure 1A; model fit statistics and parameter estimates in Supplementary Tables 3, 4). Variance of BMI was predominantly explained by genetic factors (men: H = 66%; women: H = 69%). We observed modest genetic correlations between FI and BMI, which were 0.19 (95% CI: 0.14–0.23) and 0.26 (95% CI: 0.22–0.29) in men and women, respectively. The bivariate heritability for FI and BMI, i.e., the proportion of covariance explained by genetic factors, was 81% (95% CI: 65–97%) in men and 87% (95% CI: 78–95%) in women, indicating that a substantial part of the correlation between FI and BMI could be attributable to genetic factors in common to both traits (Figure 1C). By contrast, the best-fitting bivariate model for FI and education according to the Akaike’s Information Criterion (AIC) was an ACE model (Figure 1B), although, all bivariate models between FI and education provided a worse fit than the saturated model. A moderate proportion of the variation in education was due to common environmental factors, which was 25% (95% CI: 19–31%) in men and 29% (95% CI: 23–34%) in women. The C component of FI was small and non-significant in both sexes, however, they had perfect negative correlations with that of education (r_C = -1.00). Common environmental factors also accounted for considerable proportions of the phenotypic correlations between FI and education, which were 65% (95% CI: 23–107%) in men and 74% (95% CI: 22–126%) in women (Figure 1C). There was only weak overlap of genetic factors and unique environmental factors between FI and education.

Figure 1. Parameter estimates from the best-fitting bivariate twin models. (A) ADE bivariate quantitative sex-limitation model between frailty index (FI) and body mass index (BMI), adjusted for age. Single headed-arrows represent the proportion of each traits explained by latent (circular) variance components; while double-headed arrows represent correlations between variance components. H indicates the sum of additive and dominance genetic factors; E indicates unique environmental factors; M and F are the estimates for men and women respectively. (B) ACE bivariate quantitative sex-limitation model between FI and education, adjusted for age. A indicates additive genetic factors; C indicates common environmental factors. (C) Phenotypic correlations of FI with BMI and education among men and women (with 95% confidence intervals), and the proportion of correlations explained by total genetic, common environmental and unique environmental factors. Note: Model-fitting results and parameter estimates can be found in Supplementary Tables 3, 4.

Moderation analysis

To determine whether genetic and environmental variance components of FI change with levels of BMI and education, we constructed a series of bivariate and univariate moderation models to examine whether there is moderation on the covariance between FI and the moderator (i.e., BMI and education), in addition to the variance unique to FI. For moderation by BMI, we observed significant moderating effects on the variance unique to FI as well as on the covariance between FI and BMI; in contrast, for moderation by education, we observed only significant moderating effects on the unique variance of FI (Supplementary Tables 5, 6). Therefore, a full bivariate moderation model and an extended univariate moderation model were fitted for moderation by BMI and education, respectively, as proposed by van der Sluis et al [23]. The H and E components of FI from the best-fitting moderation models are plotted in Figure 2 (see Supplementary Figure 1 for separate A and D components; and Supplementary Figure 2 for separate common and unique variance estimates of FI by BMI levels). Patterns of moderation were consistent in both men and women, although women had a higher absolute variance of FI than men. The total variance of FI was larger at low and high BMI levels and at fewer years of education. It appeared that genetic variance of FI increased at both ends of BMI levels, while unique environmental variance of FI was similar at all BMI values; hence, the heritability of FI was greater at low and high BMI. In contrast, both genetic and environmental variance of FI declined to a comparable extent with increased years of education; therefore, the relative proportions of genetic and environmental sources of FI variance did not seem to vary across education years.

Figure 2. Moderation analysis of frailty index (FI) by (A) body mass index (BMI) and (B) education, stratified by sex. First row shows the absolute variance of FI, while the second row shows the proportion of FI variance explained by total genetic (H, indicating sum of additive and dominance genetic factors) and unique environmental (E) factors, with changes in BMI and education. Variance estimates of moderation by BMI were obtained from the full ADE bivariate moderation model between FI and BMI; while variance estimates of moderation by education were obtained from the ADE extended univariate moderation model between FI and education. Quantitative sex-differences were allowed in the models to obtain estimates separately for men and women. Models were adjusted for age. Note: Model-fitting results can be found in Supplementary Tables 5, 6.

Sensitivity analysis

We additionally performed sensitivity analysis to examine if the skewed distribution of FI affected its heritability estimates. Variance component estimates of the square-root transformed FI were largely consistent with that of the non-transformed FI variable (Supplementary Table 7). Moreover, similar patterns of moderating effects by BMI and education on the transformed FI were observed in men, but there was less evidence of moderation by BMI in women (Supplementary Figure 3). We also fitted the univariate and bivariate models using an alternative “direct symmetric approach” that allows negative variance components [24]; all the conclusions were essentially the same as in the main analysis (Supplementary Tables 8–10).

Study population

Twin participants were from SALT, which is part of the population-based Swedish Twin Registry [45]. During 1998–2002, all twins born in 1958 or before were invited to attend a telephone interview with questions on common diseases, symptoms, medication use, demographics, and lifestyle factors. With a response rate of 65% for those born in 1886–1925 and 74% for those born in 1926–1958, a total of 44,919 twin individuals participated in the survey [46]. Zygosity was determined based on DNA, or questions on intra-pair similarities during childhood; the latter method was over 95% accurate when validated against DNA testing [47]. This study has been approved by the Swedish Ethical Review Authority in Stockholm. Informed consent was obtained from all participants prior to data collection.

For the present analysis, we excluded those with uncertain zygosity and those who had over 20% missing data across the 44 frailty items, leaving n = 42,994 consisting of 31,922 paired and 11,072 single respondents, aged from 41 to 99 years. Single twin individuals were retained as they contribute to the mean, variance and within-individual covariance estimates. Twins were grouped by zygosity: 4,788 MZ males (1,820 complete pairs), 5,997 MZ females (2,438 complete pairs), 7,640 DZ same-sex males (2,633 complete pairs), 8,808 DZ same-sex females (3,279 complete pairs), and 15,761 DZ opposite-sex individuals (5,791 complete pairs).

Measures

We constructed an FI based on Rockwood’s deficit accumulation model [3], using a total of 44 self-reported frailty items selected from a wide range of health-status related symptoms, signs, disabilities, and diseases in various biological systems (items listed in Supplementary Table 1). Each respondent’s frailty items were summed up and divided by the total number of deficits measured (i.e., 44 in this study), yielding a continuous FI score ranging from 0 to 1. For instance, an individual with five deficits would have an FI of 5/44 = 0.11. The FI in SALT has previously been validated for its ability to predict mortality [42].

Information on body mass index (BMI) and education were obtained from self-reported data. BMI was calculated as weight (kg) divided by height-squared (m²). We defined education as a continuous variable of the number of years of education completed.

Statistical analysis

Twin design

The classical twin design allows decomposition of phenotypic variance into genetic and environmental sources, based on the fact that MZ twins are genetically identical while DZ twins share on average half of their segregating alleles. Genetic sources of variance include additive genetic variance (A, representing the sum of allelic effects at different loci that influence the trait) and dominance genetic variance (D, representing interactions between alleles at the same locus). In MZ twins, both A and D correlate 100%; while in DZ twins, A and D are assumed to correlate 50% and 25% respectively. The sum of A and D influences is referred to as the broad-sense heritability (H), which is the proportion of phenotypic variance that can be explained by genetic effects. On the other hand, environmental sources of variance include common environmental factors (C) such as family environment and environments that are in common or correlated in adulthood, which correlate 100% in both MZ and DZ twins; and unique environmental factors (E) which are uncorrelated between twins and contribute to differences between them. The E component also contains measurement error. In addition to variance, phenotypic covariance can also be partitioned into the same sources of influence, indicating the extent to which genetic and environmental variance components overlap across two or more traits.

Phenotypic, intraclass and cross-twin cross-trait correlations

For initial inspection of the familial similarity for MZ and DZ twin pairs, we estimated intraclass correlations for FI, BMI and education, which show how much each trait correlates within twin pairs. A higher MZ correlation than DZ correlation indicates contribution of genetics to variance in the trait. If DZ correlation is less than half of the MZ correlation, the D component is inferred; on the contrary, the C component is implied if DZ correlation is greater than half of the MZ correlation. A MZ correlation less than 1 indicates presence of E.

Within-individual between-trait correlations, which we refer to as “phenotypic correlations”, were calculated for FI with BMI and education. To explore the genetic and environmental sources of covariances, we computed cross-twin cross-trait correlations, which represent the strength of relationship between one trait of the first twin and another trait of the second twin. They can be interpreted in a similar way as intraclass correlations. For example, a larger cross-twin cross-trait correlation in MZ than DZ twins suggests that phenotypic correlation between the two traits is partially owing to common genetic influences; and if MZ cross-twin cross-trait correlation is less than the phenotypic correlation then E is present.

Phenotypic, intraclass and cross-twin cross-trait correlations across zygosity were estimated from constrained saturated models in which means, variances and phenotypic correlations were equated across twin order. Age was regressed out of the means (as linear effect for FI and linear+quadratic effect for BMI and education, after initial investigation of associations between variables), so that twin correlations are not spuriously inflated due to the same age of twins in pairs.

Univariate twin modelling

Before model fitting, we checked and observed no violations in the assumptions of equal means and variances across twin order and zygosity. Univariate twin models were fitted to estimate variance components of FI, after adjusting for age. Since C and D are confounded in classical twin modelling, either an ACE or ADE model can be fitted at a time. Sex differences were allowed in the models to examine whether genetic influence on FI differs across sex quantitatively (i.e., different magnitude of heritability in men and women) and qualitatively (i.e., different genetic sources in men and women). Quantitative sex differences were modelled by allowing variance components to differ among men and women. Qualitative sex differences were estimated by multiplying a genetic correlation term (r_fm) to the expected genetic covariance of opposite-sex twin pairs, where a value of 1 would imply maximum genetic correlation (see Supplementary Methods for details). ACE or ADE sex-limitation models were then compared with their reduced models. We examined quantitative sex difference by constraining the broad-sense heritability to be equal in men and women (but allowing total variance to differ across sex) and tested if there was a significant reduction in model fit. Qualitative sex difference was assessed by testing whether r_fm was significantly less than 1. AE models were fitted to determine if the C and D parameters can be excluded from ACE and ADE models.

Bivariate twin modelling

During assumption testing, all means and variances could be equated across twin order and zygosity, except that equating means of education across zygosity led to a significant decrease in model fit, so they were estimated separately in subsequent bivariate models. A “correlated factor model” [48] was applied to estimate the etiological correlations between genetic (r_A and r_D) and environmental (r_C and r_E) variance components of FI with BMI and education, apart from estimating variance components of each variables. We define the correlation between A+D in trait X and A+D in trait Y as the broad-sense genetic correlation, r_H. Correlation of ±1 indicates complete overlap, while 0 indicates no overlap. Bivariate heritability was also calculated, which is referred to as the proportion of phenotypic correlation explained by genetic covariation. Since estimates of A and D components are highly dependent, as shown by non-significant A and D estimates, we estimated the broad-sense heritability and interpreted in ADE models to indicate the total genetic effects. Only same-sex twins were included in the bivariate quantitative sex-limitation models, and variance components and etiological correlations were estimated separately by sex, after adjusting for age.

Moderation analysis

In previous univariate and bivariate models, we assumed that variance components of FI were constant across the population (i.e., no gene-environment interaction). We relaxed this assumption and examined whether the genetic and environmental influences on FI are moderated by different levels of BMI and education. Moderation may occur on the variance that is unique to FI; for example, changes in genetic variance components of FI with years of education would suggest differential genetic sensitivity of FI to education levels. Furthermore, moderation may exist on the covariance between FI and the moderator, for example, the genetic correlation between FI and BMI may vary with different levels of BMI. To investigate these possible moderating effects, we adopted the full bivariate moderation models developed by Purcell [49], and the extended univariate moderation models proposed by van der Sluis et al [23] (Supplementary Figures 4, 5). By the bivariate moderation models, we first assessed if there were significant moderating effects on the covariance of FI with BMI and education on top of the variance unique to FI. In case of no statistically significant moderation on the covariance between FI and the moderator, the extended univariate moderation models were employed instead. All moderating effects were examined separately by sex in quantitative sex-limitation models using same-sex twins and controlled for age. Absolute genetic and environmental variance components of FI, as well as their proportions to the total variance, were plotted over BMI and education.

Sensitivity analysis

Since the FI variable is positively skewed (Supplementary Figure 6), we followed previous work and applied a square-root transformation to obtain an approximately normally-distributed FI [10, 11]. Univariate twin modelling and moderation analysis of the transformed FI were then performed to check for robustness of our results. We additionally compared our results of univariate and bivariate twin modelling with models fitted using the “direct symmetric approach” proposed by Verhulst et al. [24], which allows variance estimates to be negative and potentially reduces Type I error rates.

For all fitted models, we used the full information maximum-likelihood modelling in the R package OpenMx (version 2.17.3) for estimation of parameters that best fit the observed data. Changes in goodness of fit of models (distributed as χ²) were assessed by likelihood ratio tests, where significant values indicate worse fit of observed data. Akaike information criterion (AIC; a lower value is better) was used to select the most parsimonious and best-fitting models.