High-throughput telomere length measurement at nucleotide resolution using the PacBio high fidelity sequencing platform

Abstract

Telomeres are specialized nucleoprotein structures at the ends of linear chromosomes. The progressive shortening of steady-state telomere length in normal human somatic cells is a promising biomarker for age-associated diseases. However, there remain substantial challenges in quantifying telomere length due to the lack of high-throughput method with nucleotide resolution for individual telomere. Here, we describe a workflow to capture telomeres using newly designed telobaits in human culture cell lines as well as clinical patient samples and measure their length accurately at nucleotide resolution using single-molecule real-time (SMRT) sequencing. Our results also reveal the extreme heterogeneity of telomeric variant sequences (TVSs) that are dispersed throughout the telomere repeat region. The presence of TVSs disrupts the continuity of the canonical (5'-TTAGGG-3')n telomere repeats, which affects the binding of shelterin complexes at the chromosomal ends and telomere protection. These findings may have profound implications in human aging and diseases.

Fig. 1. Enrichment of telomere-containing genomic DNA fragments using telobaits for PacBio HiFi sequencing.

a Schematic illustration of the key steps in the enrichment of telomere repeat-containing genomic DNA fragments for PacBio HiFi sequencing. b Teloblot results showing that the telomere-containing genomic DNA fragments enriched following telobaits ligation (without trimming using RsaI and HinfI restriction endonucleases) and eluted after EcoRI digestion have similar length distributions as genomic DNA fragments digested directly using EcoRI restriction endonuclease. Representative teloblot results from one of the two independent experiments are shown (Source Data 1).

Fig. 2. Telomere length measurement in culture cells using PacBio HiFi sequencing.

a Schematic illustration for the estimation of raw telomere length or telomere length with or without the inclusion of heterogenous TVSs in sequencing reads, respectively. b, c Histograms and violin plots of telomere length distribution in HCT116, 293T, T24, IMR90 (PD47), and WI38 (PD40) cells obtained from two biological replicates (Rep.) in a single sequencing run. PD population doubling. Within each violin, the white center circle denotes the median value, the bounds of box represent the 25th to 75th percentile values, the whiskers represent adjacent values within 1.5 interquartile range, the ends of the whiskers depict the minimum and maximum values within the range, and the violin shape reflect the kernel density plot of the entire dataset. Source data are provided in Source Data 2. d. Table showing the characteristics of telomere length distribution in HCT116, 293T, T24, IMR90 (PD47) and WI38 (PD40) cells obtained from two biological replicates as shown in the violin plot in c. e Beeswarm plot showing the over-representation of extremely short telomeres (<1000 bp) in sequencing reads, specifically from T24 cells. Source data are provided in Source Data 2. f Bar graph showing the over-representation of extremely short telomeres (<1000 bp) in sequencing reads, specifically in T24 cells after normalization. g Teloblot of HCT116, 293T, T24, IMR90 (PD47), and WI38 (PD40) genomic DNA digested with the same HinfI and RsaI restriction endonucleases that were used to trim down the genomic DNA fragment during telomere enrichment for PacBio HiFi sequencing. Representative results from one of the two independent experiments are shown (Source Data 1).

Fig. 3. High-throughput telomere length measurement in patient PBL samples.

a Mean telomere length distribution of patient PBL samples shows age-dependent shortening as expected. Pearson correlation coefficient (two-sided) was used to measure the correlation between mean telomere length and age (r = −0.252; p-value = 0.023). Source data are provided in Source Data 2. b Mean raw telomere length distribution of patient PBL samples shows age-dependent shortening as expected. Pearson correlation coefficient (two-sided) was used to measure the correlation between mean raw telomere length and age (r = −0.228; p-value = 0.041). Source data are provided in Source Data 2. c Bar graph showing the percentage of different G strand end sequences across all 81 patient PBL samples. The most common G strand end observed is 5′-GGTTAG-3′. Source data are provided in Source Data 2.

Fig. 4. The unique length distribution and sequence of telomeric variant sequences (TVSs) in patient PBL samples.

a Boxplots showing the length distribution of TVSs within the telomere repeat region that is unique to each PBL sample. Please refer to the Source Data 2 file for the n value of each sample. Within each box, the red center line denotes the median value, the bounds of box represent the 25th to 75th percentile values, the whiskers represent adjacent values within 1.5 interquartile range, the ends of the whiskers depict the minimum and maximum values within the range, and gray diamonds denote outliers beyond the range. Source data are provided in Source Data 2. b The unique sequence preference of TVSs that are 1 to 6 nucleotides in length, compiled from all the telomere repeat-containing HiFi sequencing reads obtained from the 81 patient PBL samples. Only TVSs with proportion >5% are displayed. Source data are provided in Source Data 2.

Fig. 5. Presence of TVSs may uniquely affect chromosome end protection for each individual.

a Histograms showing the unique TVS size abundance profile of P01 and P02. The Y-axes are presented in logarithmic scale. Only TVSs with sizes ≤2 kb are shown. b TVS signature plots of samples P01 and P02 for telomeres at Chr. 7q. Longer and more abundant TVSs are present near the sub-telomeric ends of telomere repeat regions in the patient PBL samples as shown by an increase in TVS percentage. c The TVS signature correlation of Chr. 7q for P01 to P12. TVS signatures were pairwise compared by Euclidean distance. d Boxplots showing the density distribution of the TRF1/TRF2 binding motif on the bottom 5% of telomeres with the fewest binding sites in each sample (P01-P81). Only binding motifs occurring within the first 1 kb of canonical telomere repeats next to the sub-telomeric region were counted. Source data are provided in Source Data 2. Please refer to the Source Data file for the n value of each sample. Within each box, the red center line denotes the median value, the bounds of box represent the 25th to 75th percentile values, the whiskers represent adjacent values within 1.5 interquartile range, the ends of the whiskers depict the minimum and maximum values within the range, and gray diamonds denote outliers beyond the range. e An example of differential TRF1/TRF2 binding densities on telomeric reads at Chr. 7q from P14 and P15 PBL samples. The black lines represent telomeres starting from the 5′ end (next to the sub-telomeric region) until the end of the first 1 kb of canonical telomere repeats. The red markers represent the TRF1/TRF2 binding motifs. Source Data 2.