Genetics of physiological dysregulation: findings from the long life family study using joint models

Empirical analyses and applications of joint models

Table 1 shows the characteristics of the LLFS sample (for the probands’ and offspring generations and the total sample) including information on variables used in fitting the joint models (see section Specification of joint models). Information on time-dependent variables is presented for each visit. Information on time-independent variables is given for each individual participating in the study (whether he/she was enrolled at baseline or at follow-up visit). See Notes under the table for the number of missing values for each variable.

Table 1 also includes information on prevalent (existing before the baseline visit) and incident (new cases reported after the baseline visit) cases of major aging-related diseases in the LLFS. It shows that participants from the probands’ generation had higher prevalence of major diseases compared to the younger (offspring) generation, as expected. However, for incidence the pattern is not uniform: the proportions of new cancer cases are almost the same in two generations and the proportions of new diabetes cases tended to be higher in the offspring generation. We note also that verified information on causes of deaths of LLFS participants was not available for this study. Therefore, the numbers and proportions of individuals dying from different causes were not determined in this sample.

We constructed the measure of physiological dysregulation (denoted as $D_M^1$) from a set of biomarkers moderately correlated with age as described in Methods and separate measures based on subsets of respective biomarkers negatively and positively correlated with age ($D_M^{1-}$ and $D_M^{1+}$), as well as the “age-dependent” D_M variants (denoted, accordingly, $D_{Ma}^1$, $D_{Ma}^{1-}$, and $D_{Ma}^{1+}$) that considered deviations of biomarker values from those typical of age peers rather than those of younger individuals as in the original measures ($D_M^1$, $D_M^{1-}$, and $D_M^{1+}$) (see Methods). Table 2 lists the biomarkers used in these computations along with respective numbers of observations in each LLFS visit and the total number of observations.

We compared the average values of D_M at visit 1 between (1) those who died between visit 1 and visit 2 and (2) those who survived beyond visit 2. The average values of D_M at visit 1 for those who died were significantly higher than those for the second group (see Table 3). Similar patterns were observed when stratified by sex, with larger and more significant differences for the “original” D_M variants. These analyses indicate that the situation depicted in Figure 1 is likely to be happening here, i.e., participants with higher values of D_M and/or higher rates of change tend to die earlier thus creating the paradigm for application of appropriate statistical techniques for joint analyses of longitudinal observations of D_M and time-to-event data on mortality.

Table 4 shows the results of application of the joint models to D_M variants ${\text{D}}_{\text{M}}^{1-}$, ${\text{D}}_{\text{M}}^{1+}$, and ${\text{D}}_{\text{M}}^1$ and mortality data in LLFS. The “Longitudinal” values represent estimates (coefficients, 95% confidence intervals [CI], and p-values for the null hypotheses of zero coefficients) from the fixed effects parts of joint models. Positive (and highly significant) values of the coefficient for age at visit 1 and follow-up period since visit 1 indicate that the trajectories of D_M increase with age for all variants, similar to other studies, e.g., [23]. The estimates for sex and country indicator varied in direction and significance showing that the levels of physiological dysregulation represented by the D_M variants constructed from different sets of biomarkers can be different in females and males and can have geographical/country-specific variation. Smoking and low education have consistent negative impacts on all D_M variants (i.e., smoking and low education are associated with an increase in the level of physiological dysregulation represented by the respective D_M variants). Use of medication for diabetes and hypertension, but not lipid-lowering medications also consistently showed negative impacts on the level of physiological dysregulation. The observations for anti-diabetic and anti-hypertensive drugs may reflect the destabilizing role of underlying disease in the physiological regulation since such drugs are typically prescribed after respective diagnoses. On the contrary, lipid-lowering drugs are commonly prescribed for prevention, often without clinical manifestation of a pathology (e.g., cardiovascular disease), which might contribute to their weaker association with physiological dysregulation. The fasting variable consistently reduced the values of D_M, but it reached significance only for ${\text{D}}_{\text{M}}^1$. The “Survival” values represent estimates (coefficients, hazard ratios (HR) along with their 95% CI and p-values) from the survival parts of joint models. The results indicate that D_M (current value estimated continuously through time using the joint model) is associated with mortality (HR per standard deviation (SD) ranging from 1.66 for ${\text{D}}_{\text{M}}^{1+}$ to 2.22 for ${\text{D}}_{\text{M}}^1$) even after controlling for age and the other covariates indicated in the table (see Methods), all of which increase the risk of death with varying magnitude and significance levels.

For comparison, the results of application of the Cox model to the same data (with the respective D_M variant included as a time-dependent covariate carried forward from visit 1 and updated if there is a visit 2) are presented in Table 5. This table shows that the HRs for D_M from joint models were 1.3 to 1.4 times larger than the estimates from the Cox model.

The results for the “age-dependent” D_M variants ($D_{Ma}^{1-}$, $D_{Ma}^{1+}$, $D_{Ma}^1$) are shown in Tables 6 (joint models) and 7 (Cox model). In contrast to the “original” D_M variants ($D_M^{1-}$, $D_M^{1+}$, $D_M^1$), AgeV1, and TimeV1 have little influence on these measures (p=0.04 for $D_{Ma}^{1+}$, AgeV1; p=0.03 for $D_{Ma}^{1-}$, TimeV1; non-significant for the other cases; see Table 6, “Longitudinal”) reflecting the way such measures were constructed. Other variables in the “Longitudinal” part of Table 6 show varying significance and direction of impact on D_M. The “Survival” part of Table 6 reveals that the “age-dependent” variants still have strong associations with mortality after controlling for age and the other covariates indicated in the table (see Methods) which showed significance levels similar to the original D_M variants. However, the effect size for “age-dependent” D_M variants diminished compared to the original D_M variants (HR per SD ranging from 1.31 for $D_{Ma}^{1-}$ to 1.7 for $D_{Ma}^1$). Table 7 shows that the HRs from joint models are also larger (1.1 to 1.4 times) than the HRs from the Cox model, similar to the original D_M variants.

Genetic analyses of D_M-related traits

The R-package JM provides estimates of two types of random effects (random intercept, D_M-RI, and random slope, D_M-RS) for each individual in the analytic sample; these random effects are the traits used in the genetic analyses (see Methods). Figure 2 presents the results of the GWAS of $D_M$-RS for the "age-dependent" $D_M$ variants. We found two genome-wide significant variants on chromosome 5 in the TRIO gene (rs12652543, p=7.2×10^-9, and rs16903264, p=1.2×10^-8, which are in linkage disequilibrium (LD) [r² ~ 0.95] in whites, according to the NIH-supported online tool LDlink,https://ldlink.nci.nih.gov) that are associated with $D_M$-RS for $D_{Ma}^{1-}$. Biological interpretations of these findings are provided in Discussion. Several SNPs from this region and some others showed suggestive associations (p < 10^-5) with $D_M$-RI for $D_{Ma}^{1-}$ (Supplementary Figure 1). Analyses of other D_M variants (Figure 2: D_M-RS for $D_{Ma}^{1+}$ and$D_{Ma}^1$; Supplementary Figure 2: D_M-RS for $D_M^{1-}$, $D_M^{1+}$, and $D_M^1$; Supplementary Figure 3: D_M-RI for $D_M^{1-}$, $D_M^{1+}$, and $D_M^1$) did not yield any genome-wide significant signals, but there were several suggestive signals on different chromosomes.

The D_M variants were constructed using several biomarkers (see Methods). To check whether the signals observed in the analyses of $D_M$-RS for $D_{Ma}^{1-}$ are due to associations of the genetic variants with some particular biomarker (or a subset of biomarkers) used in the construction of $D_{Ma}^{1-}$, we performed GWAS of $D_M$ constructed from individual biomarkers following the same procedures as for $D_M$ based on several biomarkers. The results are shown in Supplementary Figure 4. As one can see from this figure, there were no signals on chromosome 5 for any of the biomarkers constituting $D_{Ma}^{1-}$. However, the GWAS of $D_M$-RS for grip strength revealed a genome-wide significant signal (see description of sensitivity analyses below) on chromosome 7 (rs10231286, p=4.6×10^-8; this SNP is located near GRM3, a glutamate receptor gene involved in brain functioning). Supplementary Table 1 shows top SNPs from this analysis.

We also performed sensitivity analyses running the models with different numbers of principal components (PCs: 1, 2, 5, 10). The top two signals for $D_M$-RS of $D_{Ma}^{1-}$ shown in Figure 2 remained genome-wide significant in all cases (p-values for the top SNP rs12652543 varied from p = 2.1 ×10^-8 to p=9.7×10^-9; p-values for the second-ranked SNP rs16903264 varied from p=3.3×10^-8 to p=1.4×10^-8). However, the genome-wide significant signal for grip strength found in the analyses with 20 PCs Hazard ratios (HR) of all-cause mortality risk for D_M variables are per standard deviation: $D_{Ma}^{1-}$: 0.474, $D_{Ma}^{1+}$: 0.385, $D_{Ma}^1$: 0.237; HR for other variables are per a unit change.

(Supplementary Figure 4) did not reach the genome-wide significant level in some models with different numbers of PCs (e.g., p=8.8×10^-8 for the model with 2 PCs). In addition, we performed analyses of D_M variants constructed using different sets of biomarkers selected according to other thresholds for correlation with age (absolute value of the correlation ≥ 0.1 and ≥ 0.2) which did not produce any genome-wide significant signals. We also ran the model with time-interaction terms for all covariates (except age) in the fixed effects part of the longitudinal sub-model in JM. All interactions were non-significant except that for country (estimate: -0.008; p = 0.03). The GWAS results for $D_M$-RS of $D_{Ma}^{1-}$ constructed from this model were similar to those shown in Figure 2 (the top two SNPs remained genome-wide significant: p = 4.3×10^-9 for rs12652543 and p = 7.4×10^-9 for rs16903264).

Applications of joint models to composite measures of physiological dysregulation (D_M) and genetic analyses of individual characteristics of D_M, in the context of research on aging

In this work, we constructed the statistical (Mahalanobis) distance measure (D_M) using multiple biomarkers collected at two visits in the LLFS using the original approach from [17] and its “age-dependent” modification that considers deviations of biomarker values from those typical of age peers rather than those of younger individuals. Analyses of longitudinal trajectories of such measures present methodological challenges (see Introduction) that require applying appropriate statistical methodology for correct inference. Here we applied one such method, joint models, for joint analyses of longitudinal observations of D_M and follow-up data on mortality for LLFS participants. Applications confirmed that, as in other studies [17, 18, 20, 22, 23], the association of D_M with mortality in LLFS is significant and effect sizes are substantial (with larger HRs for the “original” D_M variants). Comparisons of joint models with the Cox regression model with D_M considered as a time-dependent covariate indicated that the values of the association parameter for D_M in the hazard are underestimated in the Cox model, as expected [9, 10, 15]. Even though both models reveal the same direction of the influence and the result remains highly significant so that one may argue that the interpretation of the results in this case is the same in these two approaches, application of an inappropriate approach can still have substantial consequences if, for example, one needs to build a predictive model based on these results. As shown in many studies [4], including our applications to LLFS [24] and other data [23], inclusion of composite measures improves the predictive accuracy of the models for mortality and other health-related outcomes. Therefore, such predictive models should be based on an appropriate statistical approach such as joint models, which effectively account for informative missingness (death), and thus are able to correct for that type of bias.

Joint models is an active area of research in statistics with numerous extensions of the basic model (analyzed in this paper) suggested in the literature that cover a wide range of research applications such as latent classes, competing risks, multivariate models, non-linear models, dynamic predictions, stochastic processes, etc. (see books [15, 16] and recent review papers and tutorials [25–33]). Such extended models can be applied to analyze dynamic characteristics of composite measures such as D_M with various outcomes in more comprehensive ways.

One particular approach for joint analyses of longitudinal and time-to-event outcomes, the stochastic process model, SPM (see non-technical introduction in [4] and technical reviews in [34, 35]), is especially relevant in the context of research on aging as the model has components that permit clear biological interpretation in terms of fundamental features of aging-related changes in an organism. Our recent applications [22, 23] of SPM to D_M confirmed its significant association with mortality and proxy measures of physiological robustness and resilience, and revealed significant relationships of physiological dysregulation with other hidden aging-related characteristics, such as decline in stress resistance and adaptive capacity which typically are not observed in the data and thus can be analyzed only indirectly through such an analytic approach. The availability of genetic data in longitudinal studies makes it possible to explore genetic determinants of biological aging of individuals based on the dynamics of such composite measures as D_M. However, this is still a largely unexplored area. The “genetic” SPM [36–38] allows investigation of genetic determinants of such aging-related characteristics in applications to longitudinal observations of composite biomarkers (such as D_M). In particular, recent developments in the SPM methodology [39] have considerably enhanced the computational speed (which was a critical barrier in implementing this approach to large-scale genetic analyses) and opened new avenues for applying this model to GWAS, with far reaching implications for significantly improving our understanding of the genetic underpinnings of complex aging-related traits.

Interpretation of genetic associations with D_M: Biological and health effects of genes associated with the increase in physiological dysregulation with age

Table 8 shows the five genes (TRIO, FNBP1, PLXNA4, CADM1, and UBE2E2) corresponding to the top SNPs found in the association analysis of the random slope of D_M (for the “age-dependent” D_M constructed from biomarkers that decline in late life), including the SNP rs12652543 that reached genome-wide level of significance (other three genes shown in Table 8 will be discussed later). Note that multiple SNPs in these five genes were associated with slopes of D_M, many of which were in LD with each other (according to LDlink,https://ldlink.nci.nih.gov). We eliminated the redundant SNPs (r²>0.8), so the SNPs shown in Table 8 represent not only themselves but also LD blocks (r²>0.8) of other (not shown) SNPs associated with D_M slopes. We then performed in depth review of scientific literature and information provided by the NCBI Gene (https://www.ncbi.nlm.nih.gov/gene/), and found that four of the above five genes have been implemented in cancer, especially in its progression/prognosis. Three of these genes (TRIO, PLXNA4 [Plexin A4] and CADM1 [TSLC1, SynCAM1]) are also involved in axon guidance and growth (Table 8, last column) [40–46].

For a broader functional analysis, we selected 36 genes corresponding to the top 100 SNPs (all with p-value < 10^-4) from GWAS of the D_M-RS for the "age-dependent" D_M (Supplementary Table 2). We performed the pathway/process enrichment analyses for these 36 genes using several online tools available through Enrichr (https://amp.pharm.mssm.edu/Enrichr/ [47]) and MetaScape (http://metascape.org/gp/index.html#/main/step1) portals that exploit traditional ontologies and pathway sources, such as Gene Ontology [GO] processes, KEGG, BioCarta and Reactome pathway collections, among other. We also run the enrichment analysis using a commercial MetaCore platform for the functional analyses, by Clarivate Analytics (https://clarivate.com/products/metacore/), which uses custom-made manually curated libraries of pathways and processes, along with open-access ontologies, such as the GO and other [48]. We used several tools rather than just one since we wanted to feature pathways/processes that consistently show up among the top results of the enrichment analyses using the different tools.

We found that in most cases axonal guidance was among the top biological processes enriched for the 36 genes associated with the rate of increase in physiological dysregulation with age (D_M-RS) (see examples in Supplementary Figure 5). These results were further supported by the information provided by NCBI Gene (https://www.ncbi.nlm.nih.gov/gene/) about biological effects of these 36 genes, and relevant research publications (e.g., [49–51]).

Then we used the Pathway Map Creator tool, a part of the MetaCore platform [48], to create a custom map showing only products of those genes (of the 36) that participate in functionally related biological processes (Supplementary Figure 6). This figure, again, pointed to a common involvement of TRIO, PLXNA4 (Plexin A4) and CADM1 (TSLC1) in axon guidance and growth, and in cell-cell adhesion, which plays a role in both the axon guidance and cancer, and also featured the products of three more genes among those 36 (ALCAM (CD166), CNTN6 and RTN4 (Reticulon 4)) as involved in the axon guidance and nerve growth. We added these three genes to Table 8, to show their biological effects in the context of the top significant genes.

In summary, our analysis of the biological effects of the top 36 genes from GWAS of the D_M-RS, based on (i) the up-to-date scientific literature and the NCBI Gene resource, (ii) commercial (MetaCore) and open online pathway/process enrichment tools, and (iii) a custom pathway map creator (a part of the MetaCore platform), pointed to a common biological process that shows up across all these analyses, namely axon guidance. Although axon growth is mainly observed during early development, the axon guidance genes can be functional in adults and impact the maintenance of neural circuits, synaptic function and plasticity, neuroinflammatory responses, and as result neurological disorders [44, 52, 53]. Also, a recent study of the changes in human proteome across the lifespan revealed that proteins corresponding to genes involved in axon guidance and synaptic function are significantly over-represented among the clusters of proteins whose plasma levels show the strongest correlation with increasing age [54], thus supporting the role of respective biological processes in human aging.

Our results thus indicate that the decline in axons ability to form and maintain complex neuroregulatory networks may potentially play an important role in the increase in physiological dysregulation during aging.

In our recent paper we showed that the level of physiological dysregulation (estimated through D_M) can be a useful aggregate indicator of the whole-body resilience and robustness [23]. In this context, our current results of the functional analysis of genetic associations with D_M may also indicate that the declining ability to form and maintain complex neuroregulatory networks could contribute to the decline in physical resilience with age, which is the key universal feature of aging [55]. This potential connection deserves further investigation.

One should note that our results do not imply that aging can be explained by a single biological process, such as the decline in axons ability to maintain complex networks that may lead to the increase in physiological dysregulation, in turn resulting in the decline in resilience and the increase in mortality risk with age. Aging is heterogeneous, and the increase in physiological dysregulation per se is one of potentially many processes contributing to its heterogeneity. Respectively, genes associated with the decline in physiological dysregulation are not the only genes involved in aging; however, they may significantly contribute to the genetic heterogeneity of aging.

Concluding remarks

The “geroscience” hypothesis posits that interventions aimed at slowing biological aging could prevent or delay many different diseases simultaneously thus prolonging healthy lifespan and total lifespan [56]. Recent projections showed that the economic value of delayed aging (with a moderate increase in life expectancy by about 2.2 years, most of which would be spent in good health) is estimated to be $7.1 trillion over fifty years and, “in contrast, addressing major diseases such as heart disease and cancer separately would yield diminishing improvements in health and longevity by 2060 –- mainly due to competing risks” [57] providing additional arguments on the importance of identifying systemic factors that can underlie increased vulnerability to multiple diseases (rather than a specific pathology) in aging organisms. One of the critical barriers in developing interventions to slow or delay aging is that aging in humans is a gradual and slow process spanning years and decades which is not feasible to investigate in the timeframe of clinical trials. Thus, developing “proxy” measures quantifying the process of biological aging and investigation of the effects of different genetic and non-genetic factors on such measures is of paramount importance for moving research on aging forward.

Different measures to quantify biological aging (including D_M) have been recently suggested in the literature and, as the recent comparative study of such measures reveals [58], they may not measure the same aspects of the aging process thus calling for further evaluation and refinement of such measures in additional studies. This, in particular, requires rich data containing relevant information on human aging and appropriate statistical methodology that would help utilize the full potential of such data. Our previous results [23] using Framingham Heart Study and Cardiovascular Health Study data suggested that multiple deviations of biomarkers from their baseline physiological states (reflected in higher physiological dysregulation levels summarized by D_M) could be promising indicators of declining robustness and resilience during aging, and may precede clinical manifestation of not just one but many diseases (thus supporting a “geroscience” concept), even though deviations can be small and not clearly abnormal for individual biomarkers. The current paper is, to the best of our knowledge, the first study which revealed the significant genetic underpinnings of such composite measures of physiological dysregulation (D_M) in the framework of the statistical approach relevant for joint analyses of longitudinal and time-to-event outcomes (joint models).

Results of GWAS of dynamic characteristics of D_M constructed from the output of joint models yielded genes (Table 8) broadly involved in the axon guidance, synaptic function, neuroinflammatory responses, cognitive disorders and cancer, which points out to a potentially important role of the decline in neurons ability to maintain complex regulatory networks in the increase in physiological dysregulation and related mortality risk during aging.

These encouraging findings call for further exploration of the genetic mechanisms of the change in physiological dysregulation with age, and its role in the heterogeneity of human aging. They also call for future replication in independent large cohorts that collect repeated measurements of biomarkers similar to those used in the construction of composite measures of physiological dysregulation in the LLFS data.

Data

The Long Life Family Study (LLFS) is a family-based, longitudinal study of healthy aging and longevity that enrolled more than 4,900 participants from 583 families selected for exceptional familial longevity [59]. Participants were recruited at three U.S. (Boston, New York, Pittsburgh) and one European (Denmark) field centers during 2006–2009 based on age, capacity to understand the study, and their Family Longevity Selection Score (FLoSS) [59]. This score was developed specifically to select the families for the LLFS and it takes into account both the exceptionality of family members’ survival and the presence of very old living family members. The FLoSS was later validated in an independent large-scale genealogically-based resource (the Utah Population Database [60]) as a selection criterion for family longevity studies [61]. Sibships were eligible for the LLFS if their FLoSS was greater than 7 (this threshold was chosen because it was determined that such families are rare but are still detectable with sufficient frequency [59]) and they had at least one living sibling and at least one offspring willing to be enrolled in the study. Written informed consent was obtained from all subjects following protocols approved by the respective field center’s Institutional Review Boards. In this paper, we performed secondary analyses of LLFS data collected at all field centers. The data used in this study were provided by the LLFS Data Management and Coordinating Center (Washington University, St. Louis). The LLFS data are also available in the database of Genotypes and Phenotypes (https://www.ncbi.nlm.nih.gov/gap; Study Accession: phs000397.v1.p1).

Socio-demographic variables, data on past medical history and current medical conditions, medications use, physical and cognitive functioning, and blood samples were collected via in-person visits and phone questionnaires for all subjects at the time of enrollment, as described elsewhere [62]. Participants are followed-up annually to track their vital and health status. The analyses reported in this paper used the April, 20, 2018 release of LLFS data with the latest recorded date of death on January, 24, 2018. Ages at death for those participants who died within the follow-up period were computed from available dates of birth and death. Ages at censoring for those who did not die within the follow-up period were determined from dates of birth and last follow-up. The ages of the oldest participants were validated against external data [63]. Surviving participants underwent a second in-person evaluation in 2015–2018. Blood assays were centrally processed at a Laboratory Core (University of Minnesota) and protocols were standardized, monitored and coordinated through the LLFS Data Management and Coordinating Center. Genotyping was performed by the Center for Inherited Disease Research using Illumina Human Omni 2.5 v1 BeadChip array (see details on genotyping and quality control (QC) procedures in [64]).

We also reported disease-related characteristics of the analyzed sample that include the disease status at the baseline (prevalence) and new cases reported during the follow-up (incidence) for four major aging-related diseases available in the study: Alzheimer’s disease/dementia (AD), cancer, cardiovascular diseases (CVD), and diabetes. Information on diseases and health conditions was collected during the interviews either from the participants or proxies (if the participant was unable to respond). Using responses to questions about specific diseases (AD or dementia: Alzheimer’s Disease or Dementia; cancer: All cancer cites; CVD: Myocardial Infarction, Heart Attack, Coronary Angioplasty, Coronary Artery, Bypass Grafting, (Congestive) Heart Failure, Stroke, Cerebrovascular Accident, Transient Ischemic Attack, or Mini-Stroke; diabetes: Diabetes) from the baseline and the follow-up interviews, we computed the numbers of prevalent cases at the baseline and the numbers of new cases reported since the baseline.

Construction of the measure of physiological dysregulation (D_M)

The measure of physiological dysregulation (D_M) is a recently developed approach for constructing a composite measure from multiple biomarkers [17, 18, 21]. It is a continuous measure which is essentially the statistical (Mahalanobis) distance [65] from “optimality” constructed for the joint distribution of multiple biomarkers and it uses the correlation structure of the biomarkers to measure how “aberrant each individual’s profile is with respect to the overall average (centroid) of the reference population” [19]. The “reference” centroid is assumed to represent the optimal physiological state. The “reference” population can be either a subsample of the same study population or it can come from some other study. For a set of biomarkers represented by a column vector x measured in an individual at age t, x(t), D_M is defined as [17]:

$D_M(t)=\sqrt{{(x(t)-\overline{x})}^T{S}^{-1}(x(t)-\overline{x})},$(1)

where $\overline{x}$ is a vector of means and S is the variance-covariance matrix for the respective biomarkers calculated from the “reference” population, and superscript T denotes transposition.

Information on biomarkers measured in the LLFS (number of measurements at each visit, correlations with age and pairwise correlations between biomarkers, p-values for testing the null hypothesis of a zero correlation, and number of observations used for computation of correlations) is given in Supplementary Table 3. For the purpose of this paper, we initially selected a set of biomarkers collected at both visits, for a total of 30 out of 47 biomarkers available in the study (see Supplementary Table 3). We then further reduced the list of biomarkers including only those moderately correlated with age (absolute value of the correlation ≥ 0.15; see description of sensitivity analyses for different correlation thresholds in Results) to consider for the computations of the statistical distance D_M, following the ideas from previous work [17, 18]. Further, for the groups of related biomarkers (such as systolic/diastolic/pulse pressure; forced expiratory volume/forced vital capacity; red blood cell count/hematocrit/hemoglobin; total/low-density lipoprotein cholesterol; and white blood cell count/absolute neutrophil count), we randomly selected one biomarker for inclusion in D_M. We constructed the D_M variants from the resulting set of biomarkers, separating those negatively and positively correlated with age, i.e., declining versus increasing at old ages (~65+) in most people [66]. Resulting variants were denoted $D_M^{1-}$ and $D_M^{1+}$, respectively. $D_M^{1-}$ includes the following biomarkers: diastolic blood pressure, forced vital capacity, grip strength, hematocrit, and total cholesterol, and $D_M^{1+}$ includes: absolute monocyte count, creatinine, glycosylated hemoglobin, mean corpuscular volume, pulse pressure, red cell distribution width, and white blood cell count. $D_M^1$ is computed using the combined list of biomarkers from $D_M^{1-}$ and $D_M^{1+}$.

We first constructed the D_M variants ($D_M^{1-}$, $D_M^{1+}$, $D_M^1$) using the conventional approach suggested in [17]. Specifically, the observed values of each biomarker were transformed using the Box-Cox transformation and standardized so that the transformed biomarkers were all on the same scale (with a zero mean and a unit variance). These standardized and transformed values were used in calculations of D_M as in Eq. 1. We used individuals aged <60 years at the LLFS visit 1 as a “reference population.” This cutoff gave a reasonable number of participants in the reference population: 1,407 (834 females, 573 males), 1,437 (847, 590), 1,389 (821, 568), for $D_M^{1-}$, $D_M^{1+}$, $D_M^1$, respectively. Computations of the means and variance-covariance matrix in the “reference” population were performed separately for females and males using the observed values of the biomarkers included in the definition of the respective D_M. Note that in the LLFS there are at most two observations of biomarkers per individual. Therefore, we did not impute missing values of biomarkers as we did in our previous studies [22, 23] in applications to other datasets with longer series of measurements. Rather, the appropriate approach to handle missing not at random (MNAR) data is used here to jointly model the dynamics of D_M and the risk of death (see section Specification of joint models). For $D_M^{1-}$, we computed 6,244 values (4,056 at Visit 1, 2,188 at Visit 2) in 4,365 individuals. For $D_M^{1+}$, we calculated 6,557 (4,344, 2,213) values in 4,598 participants, and for $D_M^1$, the observations of 5,921 (4,010, 1,911) values in 4,290 persons were available. The resulting D_M variables were also transformed using the Box-Cox transformation. We performed sensitivity analyses (for the joint models analyses and genetic associations, see below) using the D_M in the original and the transformed scales which showed similar results so that only those for the Box-Cox transformed values are reported in the paper.

Next, we constructed the “age-dependent” D_M variants (denoted, accordingly, $D_{Ma}^{1-}$, $D_{Ma}^{1+}$, $D_{Ma}^1$) using a conceptually different approach in constructing the measure of physiological dysregulation. In the original specification [17], it was assumed that there is a fixed “normal” physiological state represented by the reference population (typically, a relatively younger sample) so that the D_M measures deviations from that state at respective ages. Here, in the definition of the “age-dependent” D_M we considered deviations of the biomarker values from those typical of age peers rather than those of younger individuals. That is, the means and variance-covariance matrix in Eq. 1 were defined in the respective age groups. Ideally, such values should be calculated for each age, however, due to small sizes of one-year age groups in our data, we computed these quantities for five-year age groups <50, 50-54, …, 90-94, 95+. For example, for a female aged 83 years at visit 1, the means and variance-covariance matrix were calculated based on a sample of females aged 80-84 years at visit 1 from which the D_M value at visit 1 was computed using the observations of biomarkers at visit 1 according to Eq. 1. Altogether, we computed 6,244 values (4,056 at Visit 1, 2,188 at Visit 2) for 4,365 individuals for $D_{Ma}^{1-}$, 6,557 (4,344, 2,213) values for 4,598 participants for $D_{Ma}^{1+}$, and 5,921 (4,010, 1,911) values for 4,290 persons for $D_{Ma}^1$.

Specification of joint models

We used joint models [15, 16] as a tool to jointly estimate the longitudinal (D_M) and time-to-event (mortality) outcomes (see Introduction). The R-package JM [67] version 1.4-8 was used to estimate the parameters of joint models. We applied the standard version of joint models as described below using the notations from [67]. The survival part of joint models quantifies the association between the longitudinal outcome and the risk of an event:

$h_i(t|M_i(t),w_i)=h_0(t)\exp \{{\gamma }^T{w}_i+\alpha m_i(t)\},$(2)

where $h_i(t|\cdot )$ is the hazard rate (mortality rate in our applications) for i-th individual at time point t, $M_i(t)=\{m_i(u),0\le u<t\}$ denotes the history of the “true” (i.e., unobserved) longitudinal outcome ($m_i(\cdot )$, see below) up to t, $h_0(\cdot )$is the baseline hazard, $w_i$ is a vector of baseline covariates, $\gamma$ is a corresponding vector of regression coefficients, and $\alpha$ is the parameter quantifying the effect of the longitudinal outcome on the risk of an event (this is usually the main parameter of interest in applications of joint models). A linear mixed effects model describes the dynamics of the longitudinal outcome (D_M in our case):

$y_i(t)=m_i(t)+{\epsilon }_i(t)=x_i^T(t)\beta +z_i^T(t)b_i+{\epsilon }_i(t),$(3)

where $y_i(t)$ is the observed longitudinal process (with available observations $y_{ij}=\{y_i(t_{ij}),j=1\dots n_i\}$ for i-th individual), $\beta$ is a vector of fixed effects parameters, $b_i$ is a vector of random (normally distributed) effects, $x_i(t)$ and $z_i(t)$ are corresponding design vectors for fixed and random effects, and ${\epsilon }_i(t)$ denotes a normally distributed error term (assumed to be independent of $b_i$) with zero mean and variance ${\sigma }^2$. The quantity $m_i(t)=x_i^T(t)\beta +z_i^T(t)b_i$ (that is, the observed value minus the error term) represents the “true” value of the longitudinal outcome which is included in the formula for the hazard rate (Eq. 2). This distinguishes joint models from the Cox model with time-dependent covariates which includes the observed values of the longitudinal process in the hazard rate. Ignoring measurement errors (or natural biological variation) in such variables in the latter model can result in underestimation of the strength of association between the hazard and the underlying longitudinal process [9, 10] which can be correctly inferred from joint models in such cases.

In our applications, the longitudinal trajectories of D_M (the longitudinal sub-model in joint models) were specified as a linear mixed effects model with linear (random intercept and random slope) random effects and time since visit 1 as a time variable (as implemented in the R-package JM). The list of covariates in the fixed effects part included: sex (1 – male, 0 – female), age at visit 1, country (1 – Denmark, 0 – USA), education (1 – below high school, 0 – otherwise), smoking (smoked >100 cigarettes in lifetime: yes [1]/no [0]), medication use (anti-diabetic, lipid-lowering, anti-hypertensive) (1 – used, 0 – did not use), and fasting (1 – ≥8 hours, 0 – otherwise). The groups of medications indicated above were constructed by the LLFS investigators in a separate study from original medications records using the corresponding Anatomical Therapeutic Chemical Classification System codes. We note that this list of medications does not include all possible groups of medications that might be relevant for this analysis (e.g., osteoporosis related medications). The time-to-event outcome (the survival sub-model in joint models) was modeled as the standard relative risk form [68] with the “true” or unobserved value of D_M (i.e., the estimate from the longitudinal sub-model [15]) included in the hazard (the “value” parameterization in the R-package JM, as in Eq. 2) along with the baseline covariates (the same list as above except medication use and fasting which are time-dependent covariates). The baseline hazard was represented as a piecewise constant function and the pseudo-adaptive Gauss-Hermite quadrature rule [69] was chosen to approximate the required integrals in the estimation algorithm (the "piecewise-PH-aGH" method in the R-package JM). We kept the default values for the number of internal knots (6 knots) in the baseline hazard and for the number of Gauss-Hermite quadrature points (3 points) used to approximate the integrals over the random effects. In some cases, the estimation algorithm in the R-package JM did not converge for the default values. In such situations, we varied the numbers (10 knots and/or 6 points) to achieve convergence. Sensitivity analyses confirmed that, in the cases of convergence, using models with different values for knots and points had little effect on the estimates of the parameter of interest (association parameter for D_M in the survival sub-model). See also sensitivity analyses with other specifications of JM described in Results.

Genome-wide association study (GWAS) of D_M-related traits

We performed GWAS of two D_M-related traits constructed using the output of the joint models estimation procedure for the respective D_M variants specified above: random intercept (D_M-RI) and random slope (D_M-RS). These characteristics were computed for each study participant and define how the baseline value of D_M (D_M-RI) and the age dynamics of D_M (D_M-RS) in the particular individual differed from the average values in the study sample (adjusted for the respective covariates, as estimated by the fixed effects part of the longitudinal sub-model in joint models).

The QC procedure was performed according to the protocol described in [70] before running the association analyses. The original data contained 4,693 genotyped individuals (2,581 females, 2,112 males) and 2,225,478 SNPs. The sample QC check removed individuals with call rate below 95% and/or heterozygosity rate beyond ± 3 standard deviations (SD) from the mean as well as individuals of divergent ancestry (those for whom the first two principal components (PC) scores were beyond 8 SD from the respective mean scores for the HapMap Phase III European reference populations). The SNP’s QC check removed duplicated SNPs, variants with missing allele code information, indels, SNPs with call rate below 95%, minor allele frequency (MAF) below 1%, and those with a significant deviation from Hardy-Weinberg equilibrium (HWE) with p-value < 10^-5. The resulting sample after QC contained 4,608 individuals of European ancestry (2,536 females, 2,072 males) and 1,464,300 autosomal SNPs which passed the QC procedure.

The R-package GENESIS (Bioconductor) [71] was used for the association testing and for computation of PCs using the PC-AiR method [72] to take into account the relatedness among individuals in the LLFS sample. The KING-robust kinship coefficient estimator [73] was used as the measure of ancestry divergence to identify a mutually unrelated and ancestry representative subset of individuals, as implemented in the PC-AiR algorithm. The mixed model was used for genetic association testing which included 20 PC-AiR PCs as fixed effects covariates to adjust for population stratification in the baseline scenario and a genetic relationship matrix (GRM) to account for genetic similarity among sample individuals. The kinship coefficient estimates from PC-Relate [74] (implemented in GENESIS) were used to construct the GRM. Genomic control (GC) [75, 76] was also applied as an additional tool to control for population stratification. In all analyses the GC lambdas were close to 1 indicating that the PCs were sufficient for this purpose.

We also performed sensitivity analyses to test the sensitivity of our results to various aspects of the analytic procedures. This included running GWAS for D_M variants computed from biomarkers selected using different thresholds for correlation with age, GWAS for separate biomarkers constituting the specific D_M variant, and analyses with different number of PCs, as described in Results.