Genetic determinants of telomere length from 109,122 ancestrally diverse whole-genome sequences in TOPMed

Summary Genetic studies on telomere length are important for understanding age-related diseases. Prior GWASs for leukocyte TL have been limited to European and Asian populations. Here, we report the first sequencing-based association study for TL across ancestrally diverse individuals (European, African, Asian, and Hispanic/Latino) from the NHLBI Trans-Omics for Precision Medicine (TOPMed) program. We used whole-genome sequencing (WGS) of whole blood for variant genotype calling and the bioinformatic estimation of telomere length in n = 109,122 individuals. We identified 59 sentinel variants (p < 5 × 10−9) in 36 loci associated with telomere length, including 20 newly associated loci (13 were replicated in external datasets). There was little evidence of effect size heterogeneity across populations. Fine-mapping at OBFC1 indicated that the independent signals colocalized with cell-type-specific eQTLs for OBFC1 (STN1). Using a multi-variant gene-based approach, we identified two genes newly implicated in telomere length, DCLRE1B (SNM1B) and PARN. In PheWAS, we demonstrated that our TL polygenic trait scores (PTSs) were associated with an increased risk of cancer-related phenotypes.


Correspondence rmathias@jhmi.edu
In brief Taub et al. leverage TOPMed wholegenome sequencing (WGS) to estimate telomere length (TL) and report a transpopulation sequencing-based association study. They identify 36 loci associated with TL and correlate polygenetic trait scores with disease outcomes.

INTRODUCTION
Telomeres shorten in replicating somatic cells, and telomere length (TL) is associated with age-related diseases. 1,2 To date, 17 genome-wide association studies (GWASs) have identified National Heart, Lung, and Blood Institute (NHLBI) Trans-Omics for Precision Medicine (TOPMed) cohorts. Our analyses of TOPMed data offer the opportunity to address the limitations of prior TL GWASs with increased sample size, population diversity, and inclusion of rare variant analyses and fine-mapping of the OBFC1 locus.
In this study, we report a sequencing-based association analysis for telomere length in 109,122 ancestrally diverse individuals (European, African, Asian, and Hispanic/Latino) from the TOPMed program. We used WGS of whole blood for variant genotype calling. We used the TelSeq method for the bioinformatic estimation of telomere length from the WGS data and demonstrated that this approach has high phenotypic and genetic correlation with laboratory-based assays, providing a reliable measurement of TL. We identified 59 sentinel variants (p < 5 3 10 À9 ) in 36 loci associated with TL; 20 of these are newly associated, and 13 replicated in external datasets. We also identified new common and rare variant associations at previously reported TL loci. Using WGS data also allowed fine-mapping approaches for OBFC1. Finally, we conducted phenome-wide association studies (PheWAS) in BioVU and identified the association of our defined polygenic trait scores (PTSs) for TL with the increased risk of cancer-related phenotypes.

RESULTS
We selected TelSeq 20 to bioinformatically estimate TL due to its computational efficiency and high correlation with Southern blot 21 and flowFISH 22 measurements (Figures S1A-S1C; STAR Methods). We developed a principal components-based approach to remove technical artifacts arising from the sequencing process that affected TL estimation, which further improved accuracy (Figures S1D and S1E; STAR Methods). We found high phenotypic correlation of TelSeq-derived TL with TL measured by Southern blot 21 in the 2,398 TOPMed samples from the Jackson Heart Study (JHS) and with TL measured by flowFISH 22 in a set of 19 TOPMed GeneSTAR samples (r = 0.68 and 0.80, respectively; Figures S1C-S1E; STAR Methods). In addition, we observed high genetic correlation between Tel-Seq and Southern blot assays of TL in the subset of 1,083 family-based JHS samples (r g = 0.8069, SE = 0.05, estimated using SOLAR 23 ). Together, the phenotypic and genetic correlations with the more traditionally used Southern blot or flowFISH assays suggest the suitability of TelSeq-based TL as a potential TL measure for large-scale genetic epidemiologic study.
Pooled trans-population association analysis was performed with n = 109,122 individuals from TOPMed (including 51,654 of European ancestry, 29,260 of African ancestry, 18,019 Hispanic/Latinos, 5,683 of Asian ancestry, and 4,506 of other, mixed, or uncertain ancestries, as determined by harmonized ancestry and race/ethnicity [HARE] 24 ; Figure 1A; STAR Methods); 44% were male and ages ranged from <1 to 98 years old (Table S1). Genome-wide tests for association were performed across 163 million variants. Using a series of single-variant tests for association (primary to identify loci, iterative conditional by chromosome to identify additional independent variants, and joint tests, including all independent variants to summarize effect sizes; see STAR Methods), we identified 59 independently associated variants mapping to 36 loci, meeting the significance threshold of p < 5 3 10 À9 ( Figure 1B; Tables 1 and S2); 16 of these were previously reported and 20 were newly associated loci, as described further below.
We examined 25 previously reported loci for TL identified through GWASs, using qPCR or Southern blot assays to directly measure TL, for evidence of replication in our study. For 16 loci (PARP1, ACYP2, TERC, NAF1, TERT, POT1, TERF1, OBFC1,  ATM, TINF2, DCAF4, TERF2, RFWD3, MPHOSPH6, ZNF208/ ZNF257/ZNF676, and RTEL1), there was at least 1 variant with p <5 3 10 À9 in our trans-population TL analysis that was in linkage disequilibrium (LD) (r 2 R 0.7) with a published genome-wide significant (p <5 3 10 À8 ) variant from a previous study (Tables 1  and S3). Directionally consistent and nominal evidence for replication was noted for CTC1 (rs3027234, p = 7.97 3 10 À5 ) and SENP7 (rs55749605, p = 0.023). A signal previously attributed to PRRC2A is located <200 kb from our signal for HSPA1A but may be distinct given the low LD (r 2 = 0.26). We found no evidence of replication (all variants with p > 0.05) for the remaining previously reported TL loci (CXCR4, PXK, MOB1B, DKK2/ PAPSS1, CARMIL1, and CSNK2A2; Table S3). Our comprehensive conditional analyses identified R1 independent sentinel variants at 9 of the 16 previously reported loci ( Figure 2A; Table  1). The resolution possible with our trans-population WGS data identified a sentinel variant different from the one previously reported by tagging-based GWASs for 11 of the 16 known loci. At known loci RTEL1, RFWD3, POT1, ACYP2, and PARP1, our WGS-based sentinels included a coding missense variant in genes RTEL1, RFWD3, POT1, TSPYL6, and PARP1, respectively. For the remaining known TL loci, many of the non-coding sentinel variants are annotated as having regulatory evidence (RegulomeDB score < 7; Table 1), as illustrated further for OBFC1 below.
A total of 22 independent sentinel variants were located at the 20 newly associated loci (Table 1). We examined 19 of these sentinel variants for evidence of association in 2 previously qPCR-based TL GWASs with non-overlapping subjects 18,19 ( Figure S2A). Variants at 10 of these loci (BCL2L15, CXXC5, HSPA1A, NOC3L, NKX2-3, ATP8B4, CLEC18C, TYMS, SAMHD1, and TYMP) had a Bonferroni-corrected p < 0.05/ 19 = 0.0026, and an additional 3 had variants with p < 0.05 (TNP03, KBTBD7, and BANP), as did a second variant at TYMS. The variant at SAMHD1 was previously reported at an false discovery rate (FDR) < 0.05 (p = 1.41 3 10 À7 ), 19 but here has genome-wide significance (p = 1.58 3 10 À19 ). Proteins encoded by two of these genes have strong biological connections to TL: CXXC5, which physically interacts with ATM and transcriptionally regulates p53 levels 25 ; two proteins implicated in telomere length regulation; and BANP (also known as SMAR1), which forms a complex with p53 and functions as a tumor suppressor. 26 There is high consistency in the effect sizes at the 59 sentinel variants observed in our TelSeq-based TL GWAS compared to prior GWASs using qPCR assays of TL ( Figure S2B). The Pearson correlation was 0.92 (p = 2.1 3 10 À15 , for 37 overlapping variants) for our study compared to Dorajoo et al. 18 (n = 23,096 Singaporean Chinese; Figure S2B, upper panel) and was 0.86 (p = 1.2 3 10 À13 , for 43 overlapping variants) for our study compared to Li et al. 19 (n = 78,592 European; Figure S2B, lower panel). However, qPCR is a relative measure of TL (relative to a single copy gene, see Note S1); it has different units from our TelSeq measurement, which is in base pairs. A direct comparison between effect sizes in base pairs from Southern blot and TelSeq in the 2,398 JHS samples (STAR Methods) confirms (B) Trans-population genome-wide tests for association using 163 million sequence-identified variants on n = 109,122 samples with sequence-generated telomere length from TOPMed. All loci had a peak p < 5 3 10 À9 in the pooled trans-population analysis. Previously reported loci for TL are indicated in red, and loci newly associated in the present study are indicated in blue. Note the shift in scale above the y axis break; no peak variants had a p value within the y axis break. Loci are labeled as known if the sentinel variants in the locus were in LD (r 2 R 0.7) with previously reported GWAS association for telomere length. There are 5 variants marked with an asterisk where the primary analysis did not meet our threshold of p < 5 3 10 À9 ; however, they reached significance after conditioning on significant variants mapping to the chromosome (detailed in Table  S2). Variants marked with a double asterisk are direct matches to prior reported sentinel variants. Percentage of trait variation explained by each variant is provided from single-variant association tests. p values and effect sizes (in base pairs) are reported from a joint model including all variants. p values for effect heterogeneity across population groups were generated using Cochran's Q statistic. MAC is the minor allele count from the full combined sample. For all exonic variants, detailed annotation is provided, while for all non-coding variants, the RegulomeDB score is given. See also Tables S2 and S4.  Table S2.
Cell Genomics 2, 100084, January 12, 2022 9 Short article ll OPEN ACCESS using GWAS summary statistics from the European ancestry group (n = 51,564) in our TOPMed WGS-based analysis and the Li et al. 19 analysis on n = 78,592 European ancestry individuals with qPCR-measured TL. The genetic correlation was 0.8066 (SE = 0.09, p = 1.8 3 10 À17 ), indicating a high degree of shared genetic determinants of TL from the 2 different measurement technologies. Each of the 59 sentinel variants individually accounted for a small percentage of phenotypic variation (Table 1), consistent with prior GWASs of TL but cumulatively accounted for 4.35% of TL variance, compared to 2%-3% from prior GWASs. 3 The 37 variants mapping to 16 known loci explained 3.38% of TL variability, and an additional 0.96% was explained by the 22 variants mapping to our 20 newly associated loci, representing a sizable gain in explained variability for TL from the present study. Prior GWASs using Southern blot and qPCR report allelic effects ranging from $49 to 120 bp. 3,4,11,13 In the TOPMed data, our estimated effect sizes for common variants (minor allele frequency, MAF R 5%) ranged from 2 to 59 bp per allele. In comparison, the effect sizes were larger for rare and low-frequency variants (MAF < 5%) in the TOPMed data (40-1,063 bp per allele).
Stratified association analyses were performed in population groups with at least 5,000 samples to evaluate effect heterogeneity of the 59 sentinel variants (Table S4). Reduced sample sizes, coupled with variation in allele frequency, often limited our power to detect population-specific associations at GWAS thresholds in the individual strata (Table S4); no additional loci were identified. A major advantage of our analysis was the ability to rely on the individual-level WGS data for the iterative conditional approach to identify the final set of independent sentinel variants at each locus. The identified sentinel variants show little evidence for heterogeneity across populations (Table 1). All Cochran's Q 29 p values (Table 1) were above a Bonferroni correction threshold (p > 0.001), and the 5 with nominal significance (0.001 < p < 0.05) appear to be primarily driven by differences in the (smallest) Asian stratum. An interesting illustration of a locus with strong allele frequency differences between groups is TINF2; the evidence at the peak variant (rs28372734) in the trans-population analysis was driven by the smaller Hispanic/Latino and Asian groups (group-specific p 4.6 3 10 À9 and 7.3 3 10 À10 , respectively), and the secondary peak (rs8016076) was driven by the African group (group-specific p 1.7 3 10 À10 ; Figure 2B; Table 1). No association is noted in the European group, where these variants are nearly monomorphic ( Figure 2C).
Gene-based tests in the combined sample of all 109,000 individuals identified 8 protein-coding genes with deleterious rare and low-frequency (MAF < 1%, including singletons) variants associated with TL (p < 1.8 3 10 À6 , see Figure S3; STAR Methods). Six of these genes were previously identified GWAS loci (POT1, TERT, RTEL1, CTC1, SAMHD1, and ATM), now adding support for rare variant associations in these genes. Both DCLRE1B and PARN have been implicated in short telomere syndrome (STS) patients. [30][31][32] DCLRE1B protein localizes to the telomere via interaction with the protein of another previously implicated GWAS gene, TERF2, and contributes to telomere protection from DNA repair pathways. 33,34 Notably, two PARN loss-of-function variants included in our gene-based test were previously identified in STS patients. 30 Both rs878853260 and rs876661305 produce frameshift mutations; rs876661305 produces an early termination codon, truncating most of the nuclease domain. 35 For each of these 8 genes, a leave-one-out approach iterating over each variant included in the aggregate test showed there were no detectable main driver variants and indicated that these gene-based association signals arise from cumulative signals across multiple rare deleterious variants (Figure S3), with the possible exception of ATM. When conditioned on the 59 sentinel variants, all of the genes, except POT1, maintained or increased statistical significance ( Figure S3). For POT1, while the removal of the single variant identified in Table 1 (rs202187871) and conditioning on all 59 sentinels resulted in a decrease in significance from 1.52 3 10 À24 to 5.53 3 10 À18 , it nonetheless remained strongly significant, meeting Bonferroni thresholds.
The identification of multiple independent sentinel variants for several loci offers the unique opportunity to evaluate the potential for distinct regulatory mechanisms (Figures 2A and S4). OBFC1 is part of a complex that binds single-stranded telomeric DNA 36 and is expressed across multiple tissues in GTEx 37 and in whole-blood studies meta-analyzed in eQTLGen. 38 All four signals at the OBFC1 locus are in the promoter and early introns of OBFC1 ( Figure 3A and 3B). Evidence for expression quantitative trait loci (eQTL) colocalization was detected at the primary, tertiary, and quaternary signals in various tissues (STAR Methods). While all 3 signals colocalized with OBFC1 eQTLs, the strongest colocalization evidence in each case was in a distinct tissue: sun-exposed skin from the lower leg (posterior probability of shared signal, PPH4 = 98.0%) for the primary, skeletal muscle (PPH4 = 84.4%) for the tertiary, and whole blood (GTEx PPH4 = 75.5%, eQTLGen PPH4 = 75.5%) for the quaternary signal (Figures 3C-3E and S5E; Table S5). Data from the Roadmap Epigenomics Consortium 39 indicate that all 4 signals are consistent with promoter or enhancer regions across blood cells and skeletal muscle tissue ( Figure 3B). We were unable to perform colocalization analysis on the secondary signal with data from either GTEx or eQTLGen as it is driven by rare variants only in the Hispanic/Latino and Asian individuals (rs111447985 ;  Table S4).
Using individual-level data within the Vanderbilt University biobank BioVU, we performed a PheWAS (Table S6) using 49 available sentinel variants individually in addition to a TL polygenic trait score (PTS). The PTS was generated separately for European and African individuals in BioVU as a simple linear combination of the effect sizes from the stratified joint analysis in European or African individuals, respectively (Table 1, Effect sizes from joint model). PTS values were significantly higher in BioVU African Americans (AAs, mean = À217 bp, SD = 96 bp) compared to European Americans (EAs, mean = À279 bp, SD = 96 bp, p < 0.05, Welch's 2-sample t test; Figure S6A), offering evidence that previously observed differences in TL by ancestry (longer TL in individuals of African ancestry 1 ) may be explained in part by the genetic contribution to TL. The largest cumulative effect of the sentinel variants, as evidenced from the PTS, is for the category of neoplasms in the EAs, with higher TL PTS associated with increased risk to the individual cancer phenotypes (11 of the 14 significant results after Bonferroni  Table S6); associations were only nominal in the BioVU AAs, likely due to lower power from the smaller sample size. Single variant PheWAS (Table S6) in the BioVU EAs are largely replicated within the UK Biobank (UKBB ; Table S7), again showing strong associations with neoplasms, and in general, demonstrating the alleles that increased TL also increased risk for these cancer-related phenotypes. In addition, analyses of both the UKBB and BioVU data identified an association between the HSPA1A locus (rs1008438) and type 1 diabetes-related endocrine/metabolism phenotypes (BioVU p 1.2 3 10 À8 to 4.1 3 10 À28 , UKBB p 6.7 3 10 À6 to 2.6 3 10 À27 for a range of phecodes grouped under 250.1); here, the allele decreasing TL increased risk for these phenotypes (BioVU odds ratios 1.4-2.1, UKBB odds ratios 1.4-2.2). This agrees with prior associations between shorter TL and increased risk of type 1 diabetes, 40 and between the protein product of HSPA1A (Hsp72) and diabetic ketoacidosis. 41

DISCUSSION
Leveraging WGS available through the NHLBI TOPMed program, we have illustrated the value of a large, trans-population WGS study for a harmonized phenotype of broad interest, bioinformatically estimated TL, to identify new loci associated with TL. The well-powered study enabled identification of rare (B) Roadmap Epigenomics Consortium data in hg19 coordinates for skeletal muscle tissue, Primary T CD4 + memory cells from peripheral blood, and primary T CD8 + naive cells from peripheral blood (Roadmap samples E108, E037, and E047, respectively; data were not available for sun-exposed skin). The ChromHMM state model is shown for the 18-state auxiliary model. The state model suggests the primary (rs9420907), secondary (rs111447985), and tertiary (rs112163720) signals are in the promoter region, while the quaternary signal (rs10883948) is in an enhancer region in all Roadmap blood cell types but is transcriptional for peripheral blood monocytes and CD19 + B cells.
(C-E) GWAS and eQTL results for the primary (C), tertiary (D), and quaternary (E) signals. The top panels are the GWAS summary statistics from the primary, and iterative conditional analyses that were used to perform colocalization analysis (secondary signal was rare and not available for colocalization). Bottom panels are eQTLs for OBFC1 in the indicated tissue from GTEx. The GTEx eQTLs for these tissues do not colocalize with one another (PPH4 < 4.4 3 10 À7 ), and each signal did not significantly colocalize in the other tissues. LD was calculated from the pooled trans-population samples with respect to the sentinel (black diamond). See also Figures S4 and S5 and Table S5.
Cell Genomics 2, 100084, January 12, 2022 11 Short article ll deleterious variants with estimated effect sizes larger than those of common variants. Using WGS allowed us the unique opportunity to hone in on causal variants using fine-mapping approaches for one locus, OBFC1, and begin to characterize tissue-specific genetic effects for this locus. We were also able to establish that for most population groups, effects are highly consistent at sentinel variants, despite differences in association strength at loci such as TINF2 and OBFC1, in which allele frequencies varied among populations. One of the main limitations to the interpretation of human genetic studies of TL pertains to the heterogeneity and lack of standardization of various TL assays, and therefore comparability of results (including the genetic effect sizes) between studies 42,43 (Note S1). To date, there is a paucity of data directly comparing WGS TL estimates with laboratory-based TL measurements in large-scale genetic epidemiologic studies. 44 In the present study, we have performed an in-depth analysis of the robustness of our TelSeq-derived TL and the resulting GWAS statistics and report the following: (1) we confirm previous observations that TelSeq estimates are consistently shorter than Southern blot (mTRF), but that the 2 values are highly correlated 44 ; (2) we demonstrate a high degree of shared heritability (i.e., genetic correlation) between TelSeq-derived and Southern blot-derived TL using phenotype-on-phenotype measures of heritability in the same subjects; (3) we see similarly high genetic correlation using GWAS summary statistic measures between qPCR-and TelSeq-derived TL GWASs in Europeans; (4) we show high correlation of effect sizes at sentinel variants between TelSeq-and qPCR-derived GWASs; and (5) we show that effect sizes from TelSeq were consistently $50% lower compared to Southern blot at the same variants measured on the same subjects mirroring the correlation patterns noted in the original phenotypes themselves.

Limitations of the study
The limitations of this study include the lack of datasets that include TL measurements in diverse populations for replication of the newly associated loci and genes; external studies are limited to largely European and Asian ancestry. For example, the identified effects at the ZMYM4 and P3H2 loci may be larger in Hispanic/Latino and African ancestry populations, respectively, than European. In both, the strength of the association, the effect sizes, and percent variation explained in the context of relative sample size in our data are larger in these non-European groups. Our lack of replication described here may be overcome with an ability to evaluate these loci in additional studies with greater population representation.
In addition, in this study, we evaluate statistical significance for association using p values that could result in anunequal ability to define significance thresholds across allele frequencies (lower allele frequencies need higher effect sizes, for example). Alternative approaches that consider effect sizes as a prioritization scheme could be applied in the future.
Finally, our fine-mapping approach based on tissue expression is limited for many of the associated loci, due to their lack of expression in GTEx tissues. Here, we followed up on the OBFC1 locus because the gene of interest is expressed in multiple adult tissues readily available in eQTL resources such as GTEx. In contrast, while TERT and TERC are important components of telomerase, they have low to undetectable expression in most GTEx tissue samples. As discussed by others, to adequately fine-map these loci, data on stem cell and/or developmental tissues will be important. 45 While our TelSeq-based TL measurements and the resulting genetic effect sizes appear to be robust based on our comparison to laboratory-based assays (summarized above), caveats described (Note S1) necessitate attention when interpreting and comparing the results between large-scale TL genetic studies, especially from the perspective of clinical risk quantification. Nonetheless, the ability to implement sequence-based TL phenotype estimation in a large, transpopulation WGS dataset creates opportunities to meaningfully expand our ability to evaluate the role of genes influencing TL in human health and disease, to dissect the genetic basis to TL differences across populations, and to set in place a model to leverage preexisting resources of WGS to bioinformatically quantify TL.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

ACKNOWLEDGMENTS
We thank Chen Li, Claudia Langernberg, and Veryan Codd for providing summary statistics from their TL GWAS for replication analysis. The whole-genome sequencing (WGS) for the Trans-Omics in Precision Medicine (TOPMed) program was supported by the National Heart, Lung, and Blood Institute (NHLBI). Specific funding sources for each study and genomic center are given in the supplemental information. Centralized read mapping and genotype calling, along with variant quality metrics and filtering were provided by the TOPMed Informatics Research Center (3R01HL-117626-02S1; contract HHSN26820 1800002I). Phenotype harmonization, data management, sample-identity QC, and general study coordination, were provided by the TOPMed Data Coordinating Center (3R01HL-120393-02S1; contract HHSN268201800001I). We gratefully acknowledge the studies and participants who provided biological samples and data for TOPMed. The full study specific acknowledgments, as well as individual acknowledgments, are detailed in the supplemental information. The BioVU projects at Vanderbilt University Medical Center are supported by numerous sources: institutional funding, private agencies, and federal grants. These include the NIH-funded Shared Instrumentation Grant S10OD017985 and S10RR025141 and CTSA grants UL1TR002243, UL1TR000445, and UL1RR024975 from the National Center for Advancing Translational Sciences. Its contents are solely the responsibility of the authors and do not necessarily represent official views of the National Center for Advancing Translational Sciences or the National Institutes of Health. Genomic data are also supported by investigator-led projects that include U01HG004798, R01NS032830, RC2GM092618, P50GM115305, U01HG0 06378, U19HL065962, R01HD074711. and additional funding sources listed at https://victr.vumc.org/biovu-funding/. Support for this work was provided by the National Institutes of Health, National Heart, Lung, and Blood Institute,

RESOURCE AVAILABILITY
Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Rasika Mathias (rmathias@jhmi.edu).

Materials availability
This study did not generate new unique reagents.

Data and code availability
TOPMed genomic data and pre-existing parent study phenotypic data are made available to the scientific community in study-specific accessions in the database of Genotypes and Phenotypes (dbGaP) (https://www.ncbi.nlm.nih.gov/gap/?term=TOPMed). Telomere length calls were derived from the raw sequence data as described in the Method details, and the phenotype covariates of age, sex, and ancestry are available through the study-specific dbGaP accession IDs as listed in the supplementary information, Table S8. This TOPMed work includes multiple studies, some of which are based on sensitive populations, precluding the unrestricted sharing of GWAS summary statistics. The TOPMed dbGaP accession (dbGaP: phs001974.v3.p1) provides a mechanism for sharing sensitive results, using the controlled-access mechanism of dbGaP to provide protections for sensitive populations. The full results from this GWAS for the full group, European, African, Asian and Hispanic/Latino subgroup analyses have all been deposited at phs001974.v3.p1 and are available publicly through dbGaP access. All original code has been deposited at Zenodo (Zenodo: https://doi.org/10.5281/zenodo.5360775) and is publicly available. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. currently consists of more than 80 participating studies, 51 with a range of study designs as described in Taliun et al. 46 Participants are mainly U.S. residents with diverse ancestries (self-reported European, African, Hispanic/Latino, Asian, and Other). Smaller representation comes from non-US populations including Samoan, Brazilian, and Asian studies. Details on the specific samples included for telomere length analysis are outlined below, details on the population groupings using HARE (harmonized ancestry and race/ ethnicity) are described in detail below, and final categories summarized in Table S1; additional information is also described by TOPMed. 51 Counts of subjects by sex are also included in Table S1. While sex is included as a covariate in all relevant models in our analysis, we do not specifically investigate the effect of sex on telomere length in this work.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
TOPMed whole-genome sequencing (WGS) WGS was performed to an average depth of 38X using DNA isolated from blood, PCR-free library construction, and Illumina HiSeq X technology. Details for variant calling and quality control are described in Taliun et al. 46 Briefly, variant discovery and genotype calling was performed jointly, across all the available TOPMed Freeze 8 studies, using the GotCloud 52 pipeline resulting in a single multistudy genotype call set.

METHOD DETAILS
Estimating telomere length for WGS samples A variety of computational tools exist that leverage WGS data to generate an estimate of telomere length. 44 Here, we performed a thorough comparison of two leading methods for estimating telomere length from WGS data to choose the preferred scalable method for performing the estimation on all available samples from TOPMed. The first method, TelSeq, 20 calculates an estimate of individual telomere length using counts of sequencing reads containing a fixed number of repeats of the telomeric nucleotide motif TTAGGG. Given that 98% of our data was sequenced using read lengths of 151 or 152 (as confirmed from the SEQ field in the analyzed CRAM files), we chose to use a repeat number of 12. These read counts are then normalized according to the number of reads in the individual WGS dataset with between 48% and 52% GC content to adjust for potential technical artifacts related to GC content. The second method, Computel 53 uses an alignment-based method to realign all sequenced reads from an individual to a ''telomeric reference sequence.'' Reads aligning to this reference sequence are considered to be telomeric and are included in the estimate of telomere length. Because Computel performs a complete realignment, additional computational steps are involved compared to those needed for TelSeq. To compare the results and scalability from these two methods, we first directly compared estimates obtained from TelSeq and Computel on 2,398 samples from the Jackson Heart Study (JHS) and found them to be highly correlated with one another (Pearson correlation r = 0.98, Figure S1A). We also compared computational time to generate the telomere length estimates on these samples and show that Computel is around ten times more time-consuming ( Figure S1B). This is in part due to the fact that Computel requires CRAM-formatted files (as the WGS data are currently stored) to first be converted back to Fastq format (while TelSeq requires a CRAM to BAM conversion), but also due to the computationally expensive step of realignment to the telomeric reference genome that the Computel algorithm employs.
TelSeq generates an estimate of TL in bp similar to laboratory assays such as Southern blot 21 and flowFISH; 22 in contrast qPCR approaches are represented as T/S ratios. 54,55 As a further comparison to orthogonally measured telomere length values, we used data on the same 2,398 samples from JHS with Southern blot 21 telomere length estimates. 56 For these samples, the Southern blot assay was performed on the same source DNA sample that was used to generate the WGS in TOPMed. The Pearson correlation values between the TelSeq and Computel estimates and the Southern blot estimates did not differ (r = 0.58 and 0.56 for TelSeq and Computel, respectively, Figure S1C). Based on our observation that both Computel and TelSeq showed similar correlation to the Southern blot estimates and high correlation with each other, and that TelSeq was an order of magnitude more computationally efficient, we chose to use TelSeq to perform telomere length estimation on our data. Final telomere length estimation was performed on a set of 128,901 samples whose CRAM-files were available for analysis at the TOPMed IRC at the time of analysis.

Batch adjustment to correct for confounders
To account for technical sources of variability in our telomere length estimates, both within a study (see, for example, colors in Figures S1A and S1C which indicate grouping by shared 96-well plate for shipment to the sequencing center) and across studies, we developed a method to estimate components of technical variability in our samples. We estimated these covariates using the sequencing data itself, similar to methods developed for other multivariate genomics data types (SVA or PEER factors 57,58 ), using aligned sequencing reads and relying on the fact that genomic coverage patterns of aligned reads can reflect technical variation.
We computed average sequencing depth for every 1,000 bp genomic region (''bin'') genome-wide using mosdepth. 59 We removed bins known to be problematic: those containing repetitive DNA sequence with difficulty mapping (mappability < 1.0 using 50bp kmers in GEMTools v1.759 60 ) or that overlap the list of known problematic SVs 61 or overlap known CNVs in the Database of Genomic Variants. To avoid overcorrecting for sex, bins were limited to autosomes. After normalizing the approximately 150,000 remaining bin counts within sample, we performed Randomized Singular Value Decomposition 62 (rSVD), a scalable alternative to principal components analysis, to generate batch principal components (bPCs). We included increasing numbers of bPCs in a linear regression model predicting TelSeq TL, and computed the correlation of the resulting residuals with external data measurements, including Southern e2 Cell Genomics 2, 100084, January 12, 2022 Short article ll OPEN ACCESS blot measurements for JHS (n = 2,398) and the Women's Health Initiative (WHI; n = 596) and age at blood draw (JHS n = 3,294; WHI n = 10,708). Based on the observed correlation, the final decision was to include the first 200 bPCs across all samples. Using the n = 2,398 JHS samples described above, we compared TL estimates before and after batch correction. The percent of variance in TL explained by sequencing plate reduced from 21.9% (baseline) to 10.5% (200 bPCs), and the variance explained by age increased from 8.0% (baseline) to 10.3% (200 bPCs), evidence that the signal-to-noise ratio was improved. Overall, the correlation between the bPC corrected TL and Southern blot data improved from r = 0.58 to 0.68 ( Figure S1D) in the JHS data and from r = 0.54 to 0.72 for the WHI data. Further, we compared TelSeq estimates of 19 samples within a single sequencing batch from the GeneSTAR study to the clinical gold standard of flowFISH 22 ( Figure S1E) and observed a correlation of 0.80 in both granulocytes and lymphocytes. Therefore, our data show that we are able to reduce the sequencing artifacts stemming from batch variability to attain correlation of TelSeq to Southern blot similar to the correlation of TelSeq to flowFISH.

Samples included in genetic analysis
All samples with telomere length estimated from the WGS data from TOPMed Freeze 8 were considered for inclusion, provided they had consent that allowed for genetic analysis of telomere length. Only samples with sequencing read lengths of 151 or 152 base pairs and having age at blood draw data available were included. For the set of samples that were part of a duplicate pair/group (either part of the intended duplicates designed by TOPMed, or a duplicate identified across the studies through sample QC) only one sample from each duplicated pair/group was retained. The final counts and demographic summary statistics for subjects grouped by TOPMed study for all 54 studies included in our analysis are shown in Table S1.
While self-reported race (Asian, Black and White) and Hispanic ethnicity group (Central American, Costa Rican, Cuban, Dominican, Mexican, Puerto Rican, South American) data are available in TOPMed, these data have limitations for analysis that include individuals with missing information or non-specific responses (e.g., 'other' or 'multiple') and high variability in genetically inferred measures of ancestry among individuals with the same reported race/ethnicity. To overcome these limitations, we used a computational method called HARE (harmonized ancestry and race/ethnicity), a newly developed machine learning approach for jointly leveraging reported and genetic data in the definition of population strata for GWAS. 24 HARE uses provided race/ethnicity labels and genetic ancestry principal component (PC) values to compute probability estimates for each individual's membership in each race/ethnicity stratum. For our HARE analysis, we used provided race (Asian, Black, White) or Hispanic ethnicity group (Central American, Costa Rican, Cuban, Dominican, Mexican, Puerto Rican, South American) as input labels to define population strata, and we used 11 PCs computed with PC-AiR 63 using 638,486 LD-pruned (r 2 < 0.1) autosomal variants with minor allele frequency > 1% to represent genetic ancestry. Genetic outliers for population strata were identified as individuals for whom their maximum stratum probability was more than 5 times greater than their reported stratum probability. Stratum values for genetic outliers and individuals with missing or nonspecific race/ethnicity were imputed as the stratum for which they had the highest membership probability.
Our primary analysis allowed for heterogeneous residual variance (see Primary single variant tests for association for details) among groups defined jointly by study and HARE-based population stratum assignment, with minor study-specific modifications to account for small strata. We required at least 30 individuals within a study-HARE grouping and collapsed individuals into merged HARE groups within a study as necessary to retain everyone for analysis. For our population-specific analyses, we used HARE assignment to stratify individuals into the following population groups: African (corresponding to the Black HARE stratum), Asian (Asian), European (White), and Hispanic/Latino (Central American, Costa Rican, Cuban, Dominican, Mexican, Puerto Rican, and South American). To better preserve genetic ancestry similarity among individuals in population-specific stratified analyses, we restricted to individuals for whom their HARE population stratum membership probability was at least 0.7; the population stratum counts in Table S1 reflect the counts in the stratified analyses, where individuals not meeting this criterion are labeled as ''Other/ Uncertain.'' Samoan individuals from the Samoan Adiposity Study and Brazilian individuals from the Reds-III Brazil study were excluded from the HARE analyses due to their unique ancestry in the TOPMed dataset; these studies were treated as their own population groups for analyses.

QUANTIFICATION AND STATISTICAL ANALYSIS
Primary single variant tests for association Genome-wide tests for association were performed using the R Bioconductor package GENESIS. 47 The primary analysis included all available trans-population TOPMed samples (n = 109,122). A secondary analysis was performed for all population groups with n > 5,000, which included European (n = 51,654), African (n = 29,260), Hispanic/Latino (n = 18,019) and Asian (n = 5,683) groups as defined above using HARE. Prior to genetic modeling, we generated residuals from a linear regression model on all 109,122 samples with 200 batch principal components (bPCs), as described above; for clarity we call these residuals TL bPC below. For the pooled trans-population analysis, we used a fully adjusted two-stage model, as described in the next two bullets. 64 For each population-specific analysis, the same approach was used, limited to samples within that population group. proportional to a sparse empirical kinship matrix computed with PC-Relate 65 to account for genetic relatedness among samples; and allowing for heteroskedasticity of residual variance across study-HARE groups as defined above. The marginal residuals from this Stage 1 model were then inverse-normalized and rescaled by their original standard deviation. This rescaling restores values to the original trait scale, providing more meaningful effect size estimates from subsequent association tests. 66 Stage 2: We fit a second LMM on all n = 109,122 samples, using the inverse-normalized and rescaled residuals from Stage 1 as the outcome; all other aspects of the model including fixed effects adjustment, random effects, and residual variance structure were identical to the model in Stage 1. This two-stage covariate adjustment has been shown to be most effective at controlling for false-positives and increasing statistical power in this setting. 64 The output of this Stage 2 model was then used to perform both single variant and gene-based tests for association.

Single variant tests for association
We used the output of the two-stage LMM to perform score tests of association for each variant with minor allele count (MAC) R 5 that passed TOPMed Informatics Research Center (IRC) at the University of Michigan quality filters 46 and which had < 10% of samples with read depth < 10. Genotype effect size estimates and percent of variability explained (PVE) were approximated from the score test results. 67 Significance, conditional analysis, locus definitions A p value cutoff of 5x10 À9 was used to determine genome-wide significance in the primary trans-ethnic analysis. We identified our set of independent significant variants (as reported in Table 1) through an iterative conditioning process within each chromosome. For a given chromosome, if at least one variant from the primary analysis crossed the genome-wide significance cutoff, this peak variant was included as an additional fixed-effect covariate in a new two-stage LMM (see Stages 1 and 2 described above), and score test results were examined to see if any remaining variants crossed the 5x10 À9 threshold. If so, we performed a second round of conditioning, including both the original peak variant and the new conditional peak variant as fixed-effect covariates in another two-stage LMM; and so on, adding conditional peak variants for additional rounds (Table S2). For each chromosome, the conditioning procedure was completed when no additional variants crossed the genome-wide threshold (p < 5x10 À9 ) on that chromosome. At each step, all variants passing the p < 5x10 À9 threshold were examined in BRAVO 68 to assess quality, and 334 variants were filtered out due to variant call quality issues. In the case where a current peak variant was flagged for quality, the next most significant variant, provided its p value was below the 5x10 À9 cutoff, was considered the peak variant instead. Variants were grouped into loci based on physical distance and an examination of linkage disequilibrium (LD) patterns, and locus names were determined using a combination of previous literature, known telomere biology, and physical location.

Cumulative percent of variability explained (PVE)
Through the iterative conditional approach, we identified a total of 59 variants ( Table 1) that met our genome-wide significance threshold of p < 5x10 À9 . The cumulative PVE values for this full set of 59 variants (4.35%), the set of 37 variants mapping to known loci (3.38%), and the set of 22 variants mapping to previously un-identified loci (0.96%, see Overlap with prior published GWAS below for definition of previously un-identified variants) were each estimated jointly using approximations from multi-parameter score tests. This joint PVE approximation is similar to the single variant PVE approximation described above, except that the set of variants is tested jointly, accounting for covariance among the estimated variant effect sizes. This approach avoids inadvertently double counting any partially shared signal among the set of identified variants.

Joint tests for association, cross-population heterogeneity
We then performed joint association analyses for the full multi-ethnic sample (n = 109,122), as well as each of the four population groups with n > 5000, to determine effect sizes and p values when all 59 variants were considered together. Using the inverse-normalized and rescaled residuals from the primary analysis Stage 1 LMM as the outcome, we fit a new Stage 2 LMM that was the same as described above, except with the additional inclusion of the genotypes for these 59 variants as additive genetic fixed effects. Given this joint modeling framework, the variant effect size estimates are all adjusted for one another. These estimates were used as input for calculation of a polygenic trait score used for the PheWAS described below. Finally, we tested for heterogeneity of effect sizes from these analyses among the population groups by adapting Cochran's Q statistic and its p value, 69 commonly used to test for effect heterogeneity in meta-analysis (Table 1). For each variant, the effect size estimates and standard errors from each population group analysis were used to calculate Q, and a Bonferroni threshold of 0.001 (0.05/59) was used to assess significance.
Overlap with prior published GWAS For each of the 59 variants identified, we examined the linkage disequilibrium (LD) with previously reported sentinel variants from 17 published GWAS. Only sentinel variants with p < 5x10 À8 in their published study were considered, which included a total of 56 variants (Table S3). If one of our variants had LD R 0.7 with a published variant, it was labeled as a known variant/part of a known locus in Table 1. Within a locus, we then compared each independent variant to the prior GWAS reported sentinel variant. If they were identical, the variant was labeled as a known sentinel variant in Table 1. Additionally, locus names for the final set of independent variants were selected based on (i) prior GWAS study definition for known loci, and (ii) the specific gene annotation for each variant mapping directly to a gene for previously un-identified loci.