Measuring telomere length with advanced sequencing helps tell healthy aging apart from disease

Next, to verify that our method captures intact, full-length telomeres we set out to use DTM to observe de novo telomere addition by forcing a cell line with stable telomeres to overexpress the catalytic core of telomerase³³. To this end, we transiently transfected HEK293T cells with a plasmid expressing an hTR encoding a variant telomeric template sequence (TSQ) or wild-type hTR in addition to TERT, or with GFP alone. We harvested genomic DNA after three days and performed telomere capture sequencing using capture oligos designed against either the canonical human telomere sequence (5′-TTAGGG-3′) or the TSQ variant sequence (5′-TTGCGG-3′). Telomere capture using the WT oligo successfully generated libraries from both experimental groups, consistent with the presence of WT telomere sequence overhangs in each sample. In contrast, attempting telomere capture with variant-targeting oligos failed to generate a successful sequencing library using genomic DNA harvested from cells that did not express the TSQ variant sequence, indicating an absence of variant telomere repeats in these samples (Supplementary Fig. 1). Transient transfection with either hTR or TSQ and TERT substantially increased in vitro telomerase activity and resulted in a 600 bp increase in the mean telomere length after 3 days in culture (Fig. 2f–i). The detection of de novo telomere elongation by expressing either wild-type or variant hTR indicates that our method measures full-length telomeres (Fig. 2g).

DTM is capable of measuring telomeres with chromosome, allele specificity provided a phase-annotated telomere-to-telomere reference assembly is available. To demonstrate this and test the reproducibility of DTM, we performed telomere capture sequencing on genomic DNA isolated from HG002 cells sourced from the Coriell Institute alongside a second independent laboratory. The summary statistics for the HG002 telomere length distributions obtained from both laboratories are in remarkable agreement (Supplementary Data 5), and the availability of an annotated diploid genome for the HG002 cell line allowed us to measure telomere length at chromosome, allele-specific resolution (Fig. 2h). We observe significant variability from the bulk telomere length distribution in a chromosome, and often allele-specific manner. For example, the mean telomere length is 6.8 kb longer on the p-arm of the paternal copy of chromosome 1 relative to the maternal copy. This demonstrates DTM will be a useful tool to further investigate chromosome-specific telomere length maintenance, inheritance, or perturbation.

To determine if our method is capable of distinguishing healthy aging from TBDs, we sequenced 7 samples of PBLs from 6 individuals aged 6–50 either diagnosed with a TBD or identified as asymptomatic heterozygous carriers of a TBD-associated RTEL1 variant²¹. In addition, we analyzed genomic DNA from two paired bone marrow (TB32-M) and peripheral blood (TB32-B) samples from one patient at Stanford Hospital with a history of chronic myelomonocytic leukemia (CMML), pulmonary fibrosis, and a newly discovered and previously unreported frameshift mutation in TIN2 (p.Gln298fs) (Supplementary Data 1). We confirmed telomeres from these individuals were significantly shorter relative to healthy individuals in their age group (Fig. 3a, c, e) and observed a significant correlation (R² = 0.84) with mean telomere length measurements made by flow-FISH, the clinical standard for telomere measurement in peripheral blood (Fig. 1f). We additionally performed TRF with DNA from twelve of the fourteen healthy donors and the two samples obtained from the Stanford patient (Supplementary Fig. 2). TRF was found to systematically overestimate mean telomere length by one to three thousand base pairs in all samples for which both digital and analog telomere measurement were performed, most likely due to undigested subtelomeric sequence biasing the results of the Southern blot (Fig. 1f, g, Supplementary Figs. 2, and 4f). Furthermore, while it was not possible to distinguish the telomere lengths from TB32's blood and marrow samples by TRF, DTM revealed the bone marrow telomere length distribution contained longer telomeres than in the same individual's PBLs (Fig. 3a, Supplementary Fig. 2). An analysis of the telomere length fractions in our cohort of individuals with defective telomere maintenance revealed significant variation between individuals but altogether much higher proportions of the shortest fractions of telomeres compared to healthy donors (Fig. 3e). Strikingly, the abundance of telomeres measuring 0–2000 bp increased in PBLs from RTEL1 TBD patients (DC1-7) but not during healthy aging in PBLs from normal donors. The abundance of telomeres measuring 0–2000 bp also increased in TERT knockout hESCs with serial passage. Interestingly, engineered patient-derived TIN2 mutations in hESCs did not lead to an increase in the abundance of short telomeres, suggesting a unique mechanism for the action of TIN2 mutations in TBDs. In TERT−/− hESCs, the rate of shortening was similar for the mean, median, 25th percentile and 75th percentile telomeres (Fig. 2e). Similarly, the slope of the relationship between telomere length and age in TBD patients was comparable for all summary statistics (Fig. 3c, red lines). In contrast, during aging in PBLs from normal donors, short telomeres (25th percentile) decreased in length 3 times more gradually than long telomeres (75th percentile)(slope = −13 bp vs −35 bp, respectively)(Fig. 3c). These data suggest that the presence of intact telomerase in hematopoietic progenitor cells preferentially maintains short telomeres, blunting the rate of shortening in normal aging. Compromising telomerase through mutations or inactivation disrupts this homeostatic mechanism, resulting in the observed similarities in the slope of the short and long telomeres in TBD patients and in TERT−/− hESCs. These data from both cultured stem cells and patients with dysfunctional telomere maintenance by telomerase suggest that intact telomerase function is needed to deter the accumulation of short telomeres.

Although it remains unclear what constitutes a critically short telomere, we know that the passage day 105 TERT−/− hESCs are 1 week away from reaching their replicative limit³². The abundance of short telomeres (between 0 and 2000 bp) in late passage TERT−/− hESCs is comparable to that observed in PBLs from TBD patients. This is evidence that bone marrow progenitors from these TBD patients are similarly approaching their replicative limit. In PBLs from normal donors, the abundance of short telomeres does not substantially increase with age (Fig. 3c). The absence of the shortest telomeres in aging PBLs suggests either that cells with very short telomeres are culled from the population or that intact telomerase actively prevents the development of these very short telomeres (0–2000 bp). These data suggest that information about the integrity of the underlying telomere maintenance machinery in a cell population can be inferred from the telomere length distribution and its change over time. Such information is not available from standard assays such as TRF Southern or Flow-FISH. Taken together, these data support a model of age-associated telomere attrition where limiting telomerase in the hematopoietic stem cell pool preferentially elongates short telomeres but is insufficient to prevent progressive shortening of the mean telomere length. This model predicts monotonic telomere shortening, on average, and more rapid shortening of longer relative to shorter telomeres with age in their cellular progeny and is consistent with our observations in our cross-sectional cohorts of human PBLs and genetically engineered hESCs.

TBDs are diagnosed by measuring PBL telomere length by flow-FISH, comparing the result to previously measured population statistics for telomere length for the patient's age, and subsequent detection of a genetic defect by targeted exon sequencing. Symptomatic patients often have observed mean telomere lengths at or below the first percentile for their age. TBDs also exhibit genetic anticipation, and the severity of clinical presentation is inversely proportional to telomere length. As a result, asymptomatic individuals carrying the same genetic defect as affected patients may still have telomere lengths within the lower quintile of the population for their age and subsequently develop symptoms later in life, leading to a delayed TBD diagnosis^19,21. DTM could also be used, therefore, as a tool in the diagnosis of TBDs. To explore this possibility, we performed a principal component analysis of age and the summary statistics of the telomere length distributions from our human samples and found this information was sufficient to distinguish healthy from affected individuals (Fig. 3f). Next, we trained a logistic regression binary classification model to distinguish the following groups by telomere lengths and age: symptomatic patients versus healthy donors; unaffected carriers versus healthy donors; or both symptomatic patients and unaffected carriers versus healthy donors. Telomere measurements and ages from our healthy donors, symptomatic patients, or unaffected carriers were randomly and iteratively assigned to serve as training or test data for the binary classification model with a 50:50 training-to-test ratio, respectively. In this small sample, our model successfully distinguished both symptomatic and asymptomatic carriers from healthy donors with high sensitivity (Supplementary Fig. 5, Supplementary Data 4). In the instance of healthy donors versus symptomatic patients, our binary classification test achieved an AUC of 0.95; when comparing healthy donors versus asymptomatic carriers, an AUC of 0.90; and when comparing healthy donors versus both symptomatic patients and asymptomatic carriers, an AUC of 0.91.

Samples from human individuals used for the preparation of this manuscript were obtained with consent and all subsequent research was conducted with all relevant ethical regulations. Study protocols were approved by IRB #13942 for the Stanford Blood Center, approving the use of genomic DNA from healthy donors for sequencing assay development, and TBD patient samples were previously collected by Sharon Savage's laboratory as a part of an ongoing NIH clinical trial (NCT00027274). The National Cancer Institute's (NCI) Institutional Review Board-approved longitudinal cohort study of inherited BMF syndromes (IBMFS) (clinicaltrials.gov NCT00027274, https://marrowfailure.cancer.gov↗) is a retrospective and prospective observational study. All participants or their legal guardians provide informed consent in accordance with Health and Human Services regulation 45 CFR 46. Detailed questionnaires are completed by participants, including personal medical and family history, and medical records are obtained and reviewed by medical experts.

This study was conducted in accordance with the criteria set by the Declaration of Helsinki.

High-molecular-weight (HMW) DNA was extracted using the NEB Monarch HMW DNA Extraction kit for Cells and Blood (Catalog #T3050L) according to the kit manufacturer's instructions. Briefly, for cells maintained in tissue culture, cells were trypsinized until detached and then centrifuged at 1000 × g for three minutes before adding prep and lysis solution according to the manufacturer's instructions. Cells were incubated for 10 minutes at 1800 RPM in a thermomixer, and DNA from the lysed cells was precipitated onto glass beads, washed twice with 80% ethanol, and finally eluted in Monarch Elution Buffer II according to manufacturer's instructions. DNA was quantified using a Qubit 4 fluorometer and the Qubit BR dsDNA quantification reagents (Catalog #Q32850). For DNA extraction from peripheral blood, red blood cell lysis was first performed prior to DNA extraction according to the manufacturer's instructions. DNA quality was assessed by Nanodrop, and the average molecular weight was verified to be 60 kb or larger using an Agilent Tapestation (Catalog #5067-5365, 5067-5365).

Approximately 4 µg of genomic DNA was prepared in a 50 uL total volume restriction digest solution (1× Fast Digest Buffer, 3 µL HinfI, and 3 µL RsaI) and allowed to digest at 37 °C overnight. In the morning, 1 µL each of HinfI and RsaI was added to each digestion reaction and allowed to incubate at 37 °C for a further 3 hours. A 1% TAE agarose gel was prepared, and 3 µL per sample underwent gel electrophoresis (125 V, 70 minutes) to confirm restriction digest completed successfully (Fig. S2). A 0.8% TBE agarose gel was then prepared in a 20 × 27 cm casting tray after adding 15 µL ethidium bromide to the agarose solution. The entire volume of each restriction digest reaction was then loaded into each well with 1× NEB nucleic acid loading dye in addition to NEB 1 kb reference ladder, and gel electrophoresis was performed (85 V for 16 hours). In the morning, the gel was dried using a BioRad gel dryer (1 hour under vacuum then 1 hour under vacuum and 50 °C). A UV-translucent ruler was then overlayed on the dried gel before imaging with UV transillumination to establish reference distances for the ladder markers from the well positions. The dried gel was then incubated in denaturing buffer (1.5 M NaCl, 0.5 M NaOH) for one hour with gentle shaking. The denatured gel was washed with deionized water twice before a second one-hour incubation in neutralizing buffer (1.5 M NaCl, 1 M Tris-HCl, pH 7.4) with gentle shaking. The neutralizing buffer was decanted, and the neutralized gel was washed twice with deionized water. The gel was then rolled vertically into a glass hybridization tube (Thermo-Fischer Scientific) and incubated with pre-warmed hybridization buffer (Invitrogen #AM8670) at 42 °C for 30 minutes with rotation. 0.5 µM of γ-³²P labeled telomere probe was then added to the hybridization buffer tube and incubated at 42 °C overnight with rotation. The gel was washed once with 2× SSC buffer and twice more with 1× SSC buffer (0.15 M NaCl, 15 mM sodium citrate) before exposing onto a phosphor screen inside a lead exposure cassette for 24 hours. Following exposure, the phosphor screen was imaged on a Typhoon scanner.

Both the southern blot images and the ethidium bromide reference ladder were loaded onto ImageJ and aligned. The signal intensities at each position coordinate starting from the bottom of the well in the southern blot image were obtained with the ImageJ line and Measure tools after drawing a line from the bottom of the well to the bottom of the gel through each sample lane.

To measure telomerase activity, a two-step TRAP procedure was performed³⁶. Cell protein extracts (at 1× or 3× dilution with lysis buffer per transient transfection condition) were incubated with telomeric primers for 30 min at 30 °C in a PCR machine, followed by 5 min of inactivation at 72 °C (cold extension). 1 μl of the cold extension reaction was PCR amplified (24 cycle of 30 s at 94 °C, followed by 30 s at 59 °C) in the presence of ³²P end-labeled telomeric primers. The radiolabeled PCR reactions were resolved by 9% polyacrylamide gel electrophoresis at room temperature, and the gel was exposed to a phosphor-imager overnight and the phosphor screen was then scanned by a Typhoon scanner.

Nanopore library preparation for whole-genome sequencing was carried out according to Oxford Nanopore Technologies (ONT) protocol for native genomic DNA sequencing (LSK-110) with some modifications. Briefly, approximately 1 µg of DNA per sample was end-prepped using the FFPE DNA repair and Ultra II End-Prep enzyme mixes from the NEBNext companion module for ONT ligation sequencing (Catalog #E7180L). The end-prep reaction was incubated in a thermocycler at 20 °C for 30 minutes and then 65 °C for 30 minutes. End-prepped DNA was extracted from the reaction using Promega ProNex size selection beads (Catalog #NG2001) at a bead-to-reaction solution ratio of 1.6 and incubated on a Hula mixer at room temperature for 5 minutes prior to being pelleted on a magnet and then washed twice with 80% ethanol and then allowed to dry on the magnet for 3 minutes. DNA was eluted from the beads using ONT elution buffer at 37 °C for 15 minutes. For simplex experiments, sequencing adapters (ONT AMXF) were then ligated to the end-prepped DNA for one hour at room temperature using NEB Quick T4 DNA ligase (Catalog #E7180L). For multiplex experiments, ONT barcodes were ligated using NEB Blunt/TA ligase for one hour at room temperature, and barcoded DNA was extracted, pooled up to 1 µg of total DNA, and then ligated to sequencing adapters as described previously. Adapter-ligated DNA was extracted from the ligation reaction using Promega ProNex size selection beads at a bead-to-reaction solution ratio of 1.1 and incubated and eluted as described previously. 20–50 fmols of Adapter-ligated DNA was sequenced on R9.4.1 PromethION flow cells on a P2Solo for 24–72 hours. It is advisable to quantify DNA after every bead purification step using a Qubit fluorometer and the Qubit dsDNA BR DNA quantification assay.

Genetically modified hESCs were sequenced on R10.4 PromethION flow cells on a P2Solo, and therefore, the sequencing adapter used was changed to pair with the updated flow cell chemistry, per ONT's standard sequencing protocol (adapter NA instead of AMII for telomere capture).

Barcoded telomere capture oligos (Supplementary Data 3) were annealed to sequencing tether (seqTether) by mixing equimolar amounts of both oligos in low TE buffer, heating to 95 °C for 2 minutes and then being allowed to cool at room temperature for an hour. Approximately 3 µg of HMW genomic DNA was ligated to barcoded, freshly duplexed oligos in a 100 µL ligation reaction (10 µL 10× rCutSmart buffer, 5 µL 5 µM duplex capture oligos, 2 µL 2000U/µL T4 DNA ligase, 1 µL 10 mM ATP, 3 µg gDNA, nuclease-free H₂O up to 100 µL) overnight at 37 °C. The following day, the ligation reaction was heat-inactivated at 65 °C for 10 minutes. Potential gaps between the capture oligo and the double-strand/single-strand junction were then filled in using the same reaction tube by adding 2 µL (4U) Sulfolubus DNA Polymerase IV (NEB #M0327S), 12 µL 10× ThermoPol Buffer, 1 µL 20 mM dNTPs, 1 µL 10 mM ATP, and 4 µL of nuclease-free H₂O followed by incubation at 56 °C for 2 minutes and then 72 °C for 15 minutes with shaking (500 RPM). Promega ProNex size selection beads were then added to the ligation reaction at a bead-to-solution ratio of 1.6, and the solution was equilibrated on a rotating mixer for 5 minutes at room temperature. The bead solution was then pelleted on a magnet, washed twice with 80% ethanol, allowed to dry on the magnet for 3 minutes, and then eluted in 60 µL of ONT elution buffer for 15 minutes at 37 °C. The eluate contains the capture oligo ligated DNA, which will be referred to as the post-capture sequencing library. The complete volume of the post-capture library was then ligated onto the ONT sequencing adapter (AMII for R9, NA for R10 libraries) per ONT's barcode ligation sequencing protocol (5 µL sequencing adapter mix, 20 µL 5× Quick T4 DNA Ligase Buffer, 10 µL Quick T4 DNA Ligase, 5 µL nuclease-free H₂O) at room temperature for one hour. The rest of the library preparation and sequencing protocol is performed as for whole-genome sequencing above. Technical replicate comparisons of R9 and R10 flow cell chemistries and basecallers are available in the supplementary information (Supplementary Figs. 8 and 9).

To prepare multiplexed telomere capture sequencing samples, post-capture libraries containing up to 8 uniquely barcoded samples were pooled by quantifying each post-capture library with the Qubit BR dsDNA quantification assay and then combining equal molecular weights of each post-capture library up to a total of 1000 ng. Sequencing adapter ligation and subsequent preparation for flow cell loading is as previously described using the 1000 ng pool of post-capture libraries as input. The most current version of both our simplex and multiplex telomere capture methodology is available at https://github.com/santiago-es/Telometer↗.

Sequencing data was basecalled using Guppy 6.3.0 (ONT) high-accuracy basecalling and aligned to a human telomere-to-telomere reference genome (T2T-CHM13 + Stong 2014 subtelomere assemblies⁴¹) with minimap2. Aligned BAM files were then sorted and indexed with samtools (v.1.16.0), and telomere measurements were extracted from alignments using a custom script. Telometer is based on previously published work for telomere measurement from PacBio HiFi reads²⁵, which was adapted for telomere measurement from Oxford nanopore long reads. In brief, reads mapping to the terminal arms of each chromosome are searched for telomeric repeats by using regular expressions targeted at the sequence patterns corresponding to both the canonical human telomeric repeat (TTAGGG/CCCTAA) and the commonly miscalled motifs previously identified in the literature⁴². A check is then performed to ensure that identified telomeric sequences are terminal and telomeres are measured from the read terminus until the sub/telomeric boundary. Conceptually, the sub/telomeric boundary is encountered by moving from the telomere terminus inward until telomeric repeats give way to non-telomeric sequence. We formally define the sub/telomeric boundary as the genomic position corresponding to the starting position of the final two telomeric repeats before encountering non-telomeric motifs. The telomere length is then the length, in base pairs, spanned by consecutive telomeric sequence motifs at the terminus of chromosome arms. Figures were made, and statistical analyses were performed with R (v4.1.0).

For samples sequenced on R10.4 PromethION flow cells, raw pod5 files were basecalled using a custom bonito telomere calling model (HG002.k1) provided by ONTs. As the R10 nanopore chemistry differs significantly from R9, the raw signal from the flow cell is substantially changed, and telomere identification from basecalled reads with our custom scripts is not possible with data produced from the default guppy R10 basecaller; however, ONT's custom telomere basecalling model once again makes it possible for our algorithm to accurately identify telomeric reads from R10 data. Training ONT basecalling models to basecall telomeres more accurately has been previously demonstrated in the literature⁴². Since work for this manuscript was completed, ONT has incorporated improved telomere basecalling into the default dorado basecaller (v0.3.4, dna_r10.4.1_e8.2_400bps_sup@v4.2.0). Telomere capture experiments with all three basecalling models on both sequencing chemistries are compared in Figs. S7 and S8. Besides the supplementary basecaller comparison, all experiments in the same figure were performed on the same sequencing chemistry and basecalled with the same model. Otherwise, library preparation and data analysis were performed as described above. At the time of writing, we recommend the R10 flow cell chemistry and using the high-accuracy Dorado basecaller for optimal results.

In order to measure telomeres at the level of individual chromosomes and alleles (Fig. 2h), we aligned HG002 sequencing reads to the diploid, allele-annotated HG002 reference genome, which is available on the Telometer github page (https://github.com/santiago-es/Telometer↗). First, one should create maternal and paternal copies of the reference genome in addition to the combined diploid genome. Reads should be aligned to both allele-specific reference genomes and analyzed with Telometer. The read id can be used to identify reads appearing in both post-analysis datasets, and only the telomere measurement corresponding to the allele with the higher mapping quality score for that read id should be preserved prior to downstream analysis.

To develop a predictor of the disease or carrier status of patients, we used LogisticRegression from Scikit-learn (v1.2.2) to train three binary classifiers on the summary statistics and age of the telomere length distributions measured in 14 healthy individuals and 8 TBD patients: healthy versus carrier, healthy versus diseased (i.e., affected carrier), and healthy versus carrier or diseased (TBD patient phenotypes are annotated in Supplementary Data ). In training, all default values were used for LogisticRegression parameters except for max_iter (the maximum number of iterations before the solver converges), which was set to 1000. The models were trained and tested using a 50:50 trainig:test data split and a receiver-operator curve were generated with Scikit-learn, representing the variation of the true or false positive rates with the internal binary classifier threshold. 1

HEK293T cells were cultured in DMEM supplemented with 10% FBS and 1% penicillin/streptomycin. Transient transfection was carried out using a 3:1 polyethylenimine to plasmid ratio and 3 µg of DNA, equal parts pCDNA-2xStrep-3xFLAG-TERT and either pBS-U1-hTR or pBS-U1-TSQ containing plasmids. Following the addition of PEI and transfection plasmids, cells were maintained in culture for three days. Cells were then trypsinized and collected, and approximately one million cells were set aside for downstream interphase DNA FISH. High-molecular-weight DNA was collected from the remaining cells for each condition according to the NEB Monarch High-Molecular Weight DNA Extraction Kit for Cells and Blood (NEB #T3050L) protocol. Extracted DNA was quantified using the Qubit BR dsDNA assay (Invitrogen #Q32853) and a Qubit 4 fluorometer.

Frozen pellets of human embryonic stem cells for TERT knockout and TIN2 mutant hESCs were generously donated by Dirk Hockemeyer's laboratory and immediately processed for sequencing without restarting culture.

HEK293T cells were sourced from ATCC (ATCC CRL-3216). Parental wild-type subclones for hESCs used to generate PARN knock-out subclones were originally sourced from WiCell (WA01). HG002 cells were sourced from the Coriell Institute (GM27730).

Further information on research design is available in the linked to this article. Nature Portfolio Reporting Summary

Supplementary Information Peer Review File Description of Additional Supplementary Files Supplementary Data 1 Supplementary Data 2 Supplementary Data 3 Supplementary Data 4 Supplementary Data 5 Reporting Summary

Source Data