The association between telomere length and mortality in Bangladesh

Subject characteristics

Characteristics of cases and controls in each cohort are summarized in Table 1 (which includes all covariates included in the regression analyses). The distribution of TL in each cohort before and after adjustments is shown in Supplementary Figure 1.

The pooled sample is 68% male and the mean baseline age is 48.4 years. Mortality cases tend to have higher arsenic exposure based on arsenic measured in urine [16,17]. The BEST cohort has the highest mean exposure while HEALS-E has the lowest. As expected, in the combined case-control sample TL was inversely associated with age in HEALS (ß=-0.034, P<2x10^-16), HEALS-E (ß=-0.034, P=1.4x10^-8), and BEST (ß=-0.030, P=2.6x10^-11). TL was not significantly associated with sex in HEALS, HEALS-E, or BEST.

Association between TL and mortality

The associations between TL and mortality outcomes were estimated using conditional logistic regression, and are summarized in Table 2. In HEALS and BEST, we observed a significant association between shorter TL and increased overall mortality (OR=0.83, P=0.03 and OR=0.79, P=0.03). The inverse association was present in both HEALS and BEST for both mortality from chronic disease (OR=0.80, P=0.01 and OR=0.78, P=0.03) and mortality from circulatory disease (OR=0.77, P=0.03 and OR=0.73, P=0.04). In HEALS-E, TL was not significantly associated with any of the mortality outcomes.

In the pooled analysis, shorter TL was significantly associated with overall mortality (OR=0.86, P=0.008), mortality from chronic disease (OR=0.85, P=0.008), and mortality from circulatory disease (OR=0.81, P=0.008). The point estimates and confidence intervals of each outcome are shown in forest plots across all cohorts in Figure 1. For the three mortality outcomes showing an association with TL (P<0.05) in pooled analyses, the association between TL quartiles and mortality is shown in Figure 2. For each outcome, the shortest TL quartile shows the most prominent association with increased mortality.

TL did not have significant interactions with sex, age, BMI, smoking, or arsenic exposure in any pooled models.

Assessment of collider bias

Because all BEST participants have arsenical skin lesions at baseline (which is part of the eligibility criteria for BEST), we risked introducing collider bias into our analysis by effectively conditioning on skin lesion status (i.e., conducting our analysis among skin lesion cases only). Skin lesion status may act as a collider because skin lesion risk is affected by arsenic exposure and may also be affected by TL (i.e., our unpublished work suggests that individuals with short TL have higher skin lesion risk, but arsenic does not have a detectable impact on TL). Using very large simulated data sets, we were able to identify the expected direction of the collider bias introduced under the causal model described in Supplementary Figure 2.

Based on analyses of simulated data, the association between TL and mortality would be biased towards higher beta values in the SL subset as compared to a random sample of the population (Supplementary Table 1). We are only able to detect this bias when we simulate data with very large effect sizes and very large sample sizes. In our main analysis, we observe an inverse association between TL and mortality, and because this association is detected even in the presence of a potential bias towards the null, collider bias is not likely to be responsible for the association between TL and mortality observed in BEST.

Study participants

HEALS is a prospective cohort study of 11,746 married men and women from Araihazar, Bangladesh. Participants were recruited between October 2000 and May 2002. Demographic data, lifestyle data, and blood samples were collected at baseline interviews. HEALS was designed to evaluate the long and short-term effects of arsenic consumed in drinking water and has been described extensively elsewhere [16]. The HEALS-E cohort was establishing using the same methods as HEALS, and participants were recruited from the same area from 2006 to 2008, with 8,287 total participants. Informed consent was obtained from all participants. The study procedures were approved by the University of Chicago and Columbia University Institutional Review Boards and the Ethical Committee of the Bangladesh Medical Research Council.

BEST is a randomized chemoprevention trial of 7,000 participants from Araihazar, Matlab, and surrounding areas. All participants have skin lesions associated with arsenic exposure. The study was created to investigate the effects of vitamin E and selenium on non-melanoma skin cancer risk and oxidative stress [23]. Participant randomization was initiated in 2006. Demographic and lifestyle data and blood samples were collected at baseline. The study procedures were approved by the Institutional Review Boards of Columbia University and the University of Chicago and the Ethical Committees of the Bangladesh Medical Research Council and the International Center for Diarrhoeal Disease Research, Bangladesh (ICDDR,B).

Ascertainment of deaths

In HEALS, mortality was evaluated at biennial in-person visits from 2000-2009. Cause of death was established using a verbal autopsy questionnaire validated by the ICDDR,B. Questionnaires were administered and reviewed by trained study physicians who determined a cause of death. Cause of death was coded using the World Health Organization’s tenth revision of the International Classification of Diseases (ICD-10) [25]. Deaths were ascertained in a similar way in BEST, but with semi-weekly household visits to alert study physicians to participant deaths [25].

Deaths from chronic disease excluded deaths from parasites and infectious disease (A00-B99), pregnancy and childbirth (O00-O99), and external or self-inflicted causes (S00-Y99). The majority of deaths from chronic disease were due to circulatory disease (I00-I99), respiratory disease (J00-J99), or neoplasms (C00-D49).

Nested case-control study sampling

To assess the association between TL and subsequent mortality, we conducted nested case-control studies within HEALS, HEALS-E, and BEST, frequency matching mortality cases in each cohort to an equal number of randomly-selected surviving controls based on sex and 5-year age intervals. Because age and sex are associated with TL, age and sex matching our case control groups avoids confounding of the association between TL and mortality by these variables [1]. Sampling was conducted in 2012 and only reflects deaths occurring prior to 2010.

In HEALS we observed 773 deaths during the study period. DNA was available for 426 of these cases, and 426 matched controls for TL measurement. After quality control (QC) exclusions (described below), we included data on 367 HEALS mortality cases and 381 HEALS controls. In HEALS-E we observed 149 deaths. DNA was available for 146 of these cases and 148 matched controls for TL measurement. After QC exclusions, we used 136 cases and 143 controls. From BEST we observed 276 deaths. DNA was available for 271 cases and 242 matched controls. Of these, TL data for 241 cases and 237 controls passed our QC filters.

Laboratory methods

In the HEALS cohort, DNA was extracted from clotted blood using Flexigene DNA kit (Cat # 51204) from Qiagen, Valencia, USA. In the HEALS-E and BEST cohorts, DNA was extracted from whole blood using QIAamp 96 DNA Blood Kit (cat# 51161) from Qiagen. Samples’ quality was checked using a Nanodrop 1000.

For HEALS, quantitative polymerase chain reaction (qPCR) was used to measure the abundance of the telomere repeat sequence in relation to a two-copy gene (RPLP0). Lab technicians were blinded to demographic and most health information of the subjects, and case and control samples were randomized across plates and well positions. We used a qPCR protocol as described by Ehrlenbach et al [7]. Briefly, ninety-six well plates were loaded with six replicates of the subject sample: three for measuring telomere content and three for measuring the reference gene (RPLP0). A reference sample was included on each plate (six replicates per plate) to account for inter-plate variation.

The PCR was performed using 30 μL amplification reaction volume of 1× Qiagen Quantitect Sybr Green Master Mix (Qiagen) and 25 ng of template DNA. For the telomere sequences, the primer sequences used were 5′GGTTTTTGAGGGTGAGGGTGAGGGTGAGGGTGAGGGT3′ (forward) and 5′TCCCGACTATCCCTATCCCTATCCCTATCCCTATCCCTA3′ (reverse). For the reference sequence, the primer sequences used were 5’CAGCAAGTGGGAAGGTGTAATCC3’ (forward) and 5’CCCATTCTATCATCAACGGGTACAA3’ (reverse). Plates were run at 95° C for 15 minutes to denature the DNA. This was followed by 25 cycles of 95° C for 15 seconds, 54° C for 2 minutes, and 72° C for 1 minute. Plates were then left at 72° C for 4 minutes.

We extracted a C_t (cycle threshold) value from the qPCR data for each sample. We used the C_t to calculate a T/S ratio comparing the abundance of telomere sequence and the reference gene sequences on each plate. By including references on each plate, we standardized this measure across all samples. A sample’s T/S ratio is proportional to its relative TL [7]. T/S ratio was calculated using the following equation:

$\frac{T}{S}\mathrm{ }rat\mathrm{i}o=\mathrm{ }\frac{2^{{C_t}_{{tel\mathrm{o}mere}^{\left(ref\mathrm{e}rence\right)}}}\mathrm{ }/\mathrm{ }{2}^{{C_t}_{{tel\mathrm{o}mere}^{\left(sam\mathrm{p}le\right)}}}}{2^{{C_t}_{{RPL\mathrm{P}0}^{\left(refer\mathrm{e}nce\right)}}}\mathrm{ }/\mathrm{ }{2}^{{C_t}_{{RPL\mathrm{P}0}^{\left(sam\mathrm{p}le\right)}}}}$

where C_t is the cycle threshold.

For QC purposes, a coefficient of variation (CV) was calculated for each set of telomere and RPLP0 replicates using the following equation:

$CV=100\mathrm{ }\times \frac{\sqrt{\frac{{\sum {(x}_i-\stackrel{-}{x})}^2}{n}}}{\stackrel{-}{x}}$

For telomere and RPLP0 triplicates measured in subject samples, we excluded a single measurement if it caused the CV to rise above 3 percent or 1 percent respectively. If after this exclusion the CV was still above the thresholds, the subject was excluded from our analysis. For telomere and RPLP0 sextuplicates measured in reference samples, we excluded up to two measurements if they caused the CV to rise above the same 3 percent or 1 percent thresholds. If after this exclusion the CV was still above the thresholds, the whole plate was excluded from our analysis. Outliers, defined as measurements more than three standard deviations away from the mean of either the raw measurement (C_t) or T/S ratio, were excluded.

We verified reproducibility of our telomere measurements by re-running 37 samples on separate days. The overall pair-wise CV between the original and later measurements of these samples was 11.7 percent. The pairwise CVs were calculated using the following equation:

$CV=100\mathrm{ }\times \mathrm{ }\frac{\sqrt{\raisebox{1ex}{${\sum (x_1-x_2)}^2$}\!\left/ \!\raisebox{-1ex}{$2n$}\right.}}{\raisebox{1ex}{$\sum (x_1-x_2)$}\!\left/ \!\raisebox{-1ex}{$2n$}\right.}$

In the HEALS-E and BEST cohorts, TL was measured using Luminex QuantiGene Plex (QGP) assay [26]. This method is strongly correlated with the Southern blot TL measures and is an accurate and reproducible method for measuring TL [26,27].

The QGP assay measures telomere sequence abundance in relation to a two-copy reference gene (ALK). The QGP protocol and its accuracy and precision have been described previously [26–28]. In summary, fluorescent Luminex microbeads were equipped with capture probes, which cooperate with several additional probes (capture extenders, label extenders, and blockers) to bind the telomere repeat sequence or the reference gene sequence. A pre-amplifier probe then binds to the label extender, and allows for the binding of biotinylated amplifiers and fluorescent labels, whose abundance can be measured by a Luminex instrument. The fluorescence intensity is proportional to the number of target sequences. The Telomere Quantity Index (TQI) for each sample is the ratio of telomere DNA measured to reference DNA measured and accounts for inter-plate variation by running a reference sample on every individual plate [26].

In a prior study, Pierce et al. calculated the CVs for replicate QGP TL measures which ranged from 5.4 to 9.1 depending on the method of calculation [27]. We assume our measurements, which followed the same laboratory protocol, are similarly reproducible. Pierce et al. found that 14% of the variation in TL was attributable to age. Our study found a comparable 10% of variation attributable to age.

Our qPCR and QGP methods of telomere measurement both produced a ratio of telomere abundance to reference gene abundance. This ratio is expressed relative to a reference DNA sample. We standardized the TL measures within each cohort to a mean of 0 and a variance of 1.

Adjustments to TL measures

Because qPCR measures tend to vary according to plate and well position, we further attempted to remove this variation from our TL measurements using mixed linear models. We created a mixed effects model predicting TL with plate as the random effect (106 plates) and well position as the fixed effect (14 well positions). The residuals produced by this model serve as an estimate of TL with experimental variation removed. We standardized the residuals to a mean of 0 and standard deviation of 1 and then included them in the models as our TL measurement.

We did not adjust for well position in HEALS-E and BEST samples because, unlike qPCR, QGP does not show well position effects [28]. We removed any effect of plate on QGP TL using a linear regression; we predicted TL with categorical variables representing “plate” as predictors (52 plates) and then removed plate’s effect on TL.

Statistical methods

All analyses accounted for the following covariates: Age (years), sex (male or female), body mass index (BMI) quartiles (<17.5, 17.5-19.1, 19.1-21.5, >21.5), smoking status (never, former, current), arsenic concentration in urine quartiles (creatinine-adjusted urinary arsenic) (<87.87, 87.87-159.9, 160.0-298.2, >298.2), years of formal education (0, 1-5, 6-16), land ownership, and television ownership. BMI and arsenic exposure quartiles were based on the distribution among controls pooled across all cohorts.

We estimated the association between TL and mortality using conditional logistic regressions in order to account for the matched case-control design. Regressions were conditioned on the sex-specific, five-year age strata used to match controls to cases. In instances where within-strata sample sizes were too small to allow for matching, we combined adjacent age strata within sex categories.

Each cohort was analyzed separately for overall mortality and for mortality due to chronic disease, respiratory disease, circulatory disease, and neoplasm. In each cohort, subjects with TL measurement were discarded if missing data on covariates included in the model. We also conducted pooled analyses including all cases and controls from HEALS, HEALS-E, and BEST, adjusting for cohort as a categorical variable. In these pooled analyses we conducted interaction tests for TL with the following mortality-related covariates: age, sex, BMI, smoking, and arsenic exposure.

Simulation to assess potential collider bias

All BEST participants have arsenical skin lesions. Because skin lesion risk is effected by arsenic exposure and may be effected by TL, skin lesion status may act as a collider. We ran a simulation to measure any potential bias introduced by conditioning on skin lesion status in BEST. By simulating a very large dataset, we were able to determine the expected direction of the potential collider bias under the causal model described in SSupplementary Figure 2. We simulated 100 datasets of 1,000,000 observations and four variables: TL (W) and binary variables representing arsenic exposure (X), skin lesion status (Y) and mortality (Z). From each simulated dataset we additionally sampled 250 mortality cases and 250 mortality controls so that we could simultaneously simulate what we would expect to see in a dataset the size of our BEST data. We did this for both the overall dataset and the subset of SL cases. TL was generated as a random standard normal variable. Arsenic was generated as a random binary variable with a 0.3 probability of being high exposure (x=1). The skin lesion (Y) and mortality (Z) risk were influenced by W and X. Y and Z were generated as random binary variables using the following equations to determine their respective probabilities of a case:

$P(y=1)=\mathrm{ }\frac{1}{1+e^{-1\mathrm{*}({\mathrm{ß}}_{y0}+\left({\mathrm{ß}}_{yw}\mathrm{*}W\right)+\left({\mathrm{ß}}_{yx}\mathrm{*}X\right))}}$

$\mathrm{P}\left(z=1\right)=\mathrm{ }\frac{1}{1+e^{-1\mathrm{*}({\mathrm{ß}}_{z0}+\left({\mathrm{ß}}_{zx}\mathrm{*}X\right))}}$

We repeated the simulation three times with weak effect sizes similar to what we observed in our data, moderate effect sizes, and strong effect sizes. The effects of arsenic exposure on skin lesion and mortality (${\mathrm{ß}}_{YX}$ and ${\mathrm{ß}}_{ZX}$) were set to 0.405, 0.833, and 1.609 (OR=1.5, 2.3, and 5), respectively. ${\mathrm{ß}}_{YW}$ was set to -0.289, -0.693, and -0.1609 (OR=0.75, 0.5, and 0.2), and ${\mathrm{ß}}_{ZW}$ was set to 0 (OR=1) in all simulations.

We defined the baseline risk of skin lesions to be 5 percent (${\mathrm{ß}}_{Y0}$=-log((1-0.05)/0.05)) and baseline risk of mortality to be 3 percent (${\mathrm{ß}}_{Z0}$=-log((1-0.03)/0.03)). The association between TL and mortality was estimated using a simple logistic regression for the entire sample and just skin lesion cases.