Non-random length distribution of individual telomeres in immunodeficiency, centromeric instability and facial anomalies syndrome, type I

Abstract

Mutations in the de novo DNA methyltransferase DNMT3B lead to Immunodeficiency, Centromeric Instability and Facial anomalies (ICF) syndrome, type I. This syndrome is characterized, among other hypomethylated genomic loci, by severe subtelomeric hypomethylation that is associated with abnormally short telomere length. While it was demonstrated that the mean telomere length is significantly shorter in ICF type I cells, it is unknown whether all telomeres are equally vulnerable to shortening. To study this question we determined by combined telomere-FISH and spectral karyotyping the relative length of each individual telomere in lymphoblastoid cell lines (LCLs) generated from multiple ICF syndrome patients and control individuals. Here we confirm the short telomere lengths, and demonstrate that telomere length variance in the ICF patient group is much larger than in the control group, suggesting that not all telomeres shorten in a uniform manner. We identified a subgroup of telomeres whose relatively short lengths can distinguish with a high degree of certainty between a control and an ICF metaphase, proposing that in ICF syndrome cells, certain individual telomeres are consistently at greater risk to shorten than others. The majority of these telomeres display high sequence identity at the distal 2 kb of their subtelomeres, suggesting that the attenuation in DNMT3B methylation capacity affects individual telomeres to different degrees based, at least in part, on the adjacent subtelomeric sequence composition.

Introduction

Telomeres are DNA–RNA-protein complexes that play essential roles in maintaining the integrity and stability of chromosome ends [reviewed in (1–4)]. Telomeric DNA in mammals is composed of the hexameric repeat (TTAGGG)_n, and is complexed with shelterin and non-shelterin proteins that play various roles in telomere maintenance and function (5). In addition, telomeres contain the long non coding RNA, TERRA (telomeric repeat-containing RNA), transcribed from a subset of subtelomeres in human cells [(6,7), reviewed in (3)]. TERRA has been shown to bind telomeric regions and to be associated with proper telomere functioning through different mechanisms (8) [reviewed in (1,9,10)]. TERRA has also been shown to form DNA:RNA hybrids that promote recombination in telomerase negative cells, and can affect telomere stability of human cancer cells utilizing the alternative lengthening of telomeres (ALT) mechanism (11) [reviewed in (3)].

Telomere length homeostasis is the end result of opposing forces that lengthen and shorten telomeres (12). Telomerase, an RNA-protein enzymatic complex that elongates telomeres by reverse transcription, is the main elongating force in germ cells, embryonic and adult stem cells, and in the majority of human cancer cells (13,14). Telomerase preferentially elongates short telomeres (15,16), and in yeast TERRA has been shown to be involved in this preferential recruitment to short telomeres (17,18). The most characterized shortening force in normal dividing cells is the result of the ‘end replication problem’ occurring at the chromatid replicated by lagging strand synthesis, and this shortening affects all telomeres at a relatively constant rate [(19–22), reviewed in (12,23)]. Also associated with cell division is shortening due to telomere end processing at the chromatid replicated by leading strand synthesis, involving the Apollo nuclease and the 5’-3’ exonuclease, Exo1 (24). Stochastic shortening of telomeres is attributed to sporadic events of DNA damage due to environmental factors such as smoking and pollution (23), as well as intracellular factors such as oxidative stress, which can contribute to shortening in a non-uniform manner and to variable degrees [reviewed in (25)]. Recently, an additional shortening force in the form of telomere trimming has been characterized. Telomere trimming plays an important role in telomere length regulation by shortening overly long telomeres [reviewed in (26)] and has been shown to be carried out by XRCC3 and Nbs1 in embryonic stem cells and induced pluripotent stem cells (iPSCs) (27), and by TZAP/ZBTB48 in various other mammalian cells (28,29).

In normal human replicating cells lacking telomerase activity, when the shortening telomeres reach a critical length, they trigger entry into senescence, a cellular state in which the cell exits the cell cycle permanently (30). Numerous studies support the notion that a small subset of critically short telomeres is sufficient to drive cells into the senescent state (15,31–33).

Several Mendelian disorders display abnormally short telomere length as part of their phenotype [reviewed in (23,34)]. These ‘telomere spectrum disorders’ include primary telomeropathies, such as Dyskeratosis Congenita, in which mutations affect genes fundamental for telomere homeostasis (23,35). The secondary telomeropathies involve genes which encode proteins frequently involved in DNA repair, and the abnormal telomeric phenotypes in these diseases may be rescued by ectopic expression of telomerase (23). In addition, in at least one Mendelian disease, Immunodeficiency, Centromeric instability and Facial anomalies (ICF) syndrome, type I (OMIM #242860), mutations in a gene encoding a protein involved in epigenetic modification of the genome leads to a telomeric phenotype (36,37).

The chromatin structure at chromosome ends constitutes an additional layer of protection of telomere integrity (38). Telomeres and subtelomeres, the regions immediately proximal to telomeres, are packaged together as heterochromatin, characterized by typical histone modifications and high levels of subtelomeric-DNA methylation (38,39). The majority of human subtelomeres contain CpG-rich repetitive sequences that undergo de novo methylation by DNA methyltransferase 3B (DNMT3B) during early embryonic development (40–42). In the rare autosomal recessive ICF syndrome type I, bi-allelic mutations in DNMT3B result in insufficient residual activity of DNMT3B (43,44). This syndrome is often fatal in childhood due to markedly reduced serum immunoglobulin levels [reviewed in (45,46)], however the exact mechanism/s by which compromised methylation leads to this immunodeficiency is as yet unclear. ICF patient DNA exhibits hypomethylation at various genomic regions, among them repetitive regions such as satellite 2 and 3 repeats that consequently fail to condense properly in mitosis (47).

DNMT3B-mutated ICF syndrome cells display, among other repetitive regions, severe hypomethylation at subtelomeres. The subtelomeric hypomethylation is associated with accelerated telomere shortening, significantly elevated levels of TERRA (36,37) and enhanced levels of subtelomeric/telomeric DNA:RNA hybrids (48). In accordance with abnormally short telomere length, ICF fibroblasts enter premature replicative senescence at a low population doubling (PD), a phenomenon that can be rescued by ectopic expression of the catalytic subunit of human telomerase (hTERT) (49). Despite endogenous expression of telomerase, telomeres are also abnormally short in lymphoblastoid cell lines (LCLs) generated from ICF type I patients (37), indicating that the level of telomerase in these cells is insufficient to compensate for the accelerated shortening occurring in this syndrome.

Telomere-FISH analysis of ICF syndrome metaphase cells reveals absence of hybridization signals from either one or both sister chromatids at an abnormally high frequency (37,50). However, only a fraction of telomeres in each cell demonstrates this abnormality, suggesting that not all telomeres shorten evenly. While the short mean telomere length has been demonstrated clearly by terminal restriction fragment (TRF) analysis in several patients (37), it is unknown whether all telomeres are equally prone to shortening in this disease. A scenario in which compromised methylation at subtelomeric regions or elsewhere in the genome may place certain telomeres at a higher or lower risk to shorten would result in a non-random distribution of telomere lengths across several ICF patients.

Here, we study whether all telomeres are equally vulnerable to shortening in DNMT3B-mutated ICF syndrome cells by analyzing individual telomere lengths in ICF versus control LCLs. We find significantly enhanced telomere-length heterogeneity in the ICF LCLs, suggesting that telomere shortening in this disease does not occur at a constant rate for all telomeres per each cell division. We find that this unequal shortening among the telomeres is not random. Using various statistical approaches to analyze the telomere length data acquired in this study, we identified a subgroup of telomeres whose frequently short lengths can distinguish between a control and ICF metaphase. This suggests that in ICF syndrome cells, certain individual telomeres are at greater risk than others to shorten.

Results

Methylation analysis of chromosome end-specific subtelomeres in LCLs

To confirm that indeed all ICF type I LCLs subjected to analysis in this study display an abnormal telomeric phenotype, we first examined their average telomere length and subtelomeric methylation status (all further reference to ICF patients or cells in this study relates to type I solely). Each of the four ICF LCLs, P1-P4 (Table 1), carries different mutations in DNMT3B (51). Out of these four ICF LCLs, abnormally short telomere length and subtelomeric hypomethylation had been demonstrated previously by TRF analysis for the three patient LCLs P1, P2 and P3 (37). To demonstrate that ICF patient P4 LCL shows a similar telomere length phenotype, we performed a TRF analysis on these cells (Supplementary Material, Fig. S1) which indeed confirmed that telomeres in ICF P4 LCL are abnormally short. In addition, we could demonstrate by this approach that the telomere lengths of ICF patient P3 are strikingly low in comparison to his parents (carriers C1 and C2) (Supplementary Material, Fig. S1). Similarly, we had previously demonstrated that the telomere lengths of ICF P1 LCL (pCor) are significantly shorter than those of her parents, C3 and W1 (37) [the mutation transferred by W1 is a de novo germline mutation (52)].

Table 1.

Lymphoblastoid cell lines (LCLs) utilized in this study

LCL Sample	Additional Name and/or Catalogue ID	DNMT3B Genotype	Cell Source	Additional Information	Reference
P1	GM08714 (pCor)	ICF	Coriell Institute, USA	Daughter of C3 and W1	52
P2	pG	ICF	Claire Francastel and Guillaume Velasco, Paris Diderot University		51
P3	pY	ICF	,,	Son of C1 and C2	51
P4	pT	ICF	,,		51
C1	mY	Carrier	,,	Mother of P3	48
C2	fY	Carrier	,,	Father of P3	48
C3	GM08728 (mCor)	Carrier	Coriell Institute, USA	Mother of P1	52
W1	GM08729 (fCor)	WT	,,	Father of P1	52
W2	GM119116c	WT	,,		48
W3	3125	WT	Nadine Cohen, Technion, Haifa, Israel		79

We extended the methylation analysis of subtelomeric regions by bisulfite sequencing of several subtelomeres for all ICF LCLs, and for three carrier LCLs (C1-C3) and three WT LCLs that served as controls in this study (W1-W3). As shown previously for ICF fibroblasts and iPSCs (50), the seven amplicons that represent 12 subtelomeric regions demonstrate strikingly lower methylation levels in ICF LCLs compared to all control samples (Fig. 1A). Some variation in subtelomeric methylation can be observed between the ICF LCL group members, with LCL P3 displaying slightly higher methylation levels (Fig. 1B). This variation may be influenced by the level of residual DNMT3B activity, and/or by the effect of modifier loci. Additionally, these findings reveal that the subtelomeric methylation levels of two (C1 and C2) out of three carrier LCLs are highly similar to those of WT LCLs (W1-W3). The subtelomeric methylation levels of the LCL derived from carrier C3, the mother of ICF patient P1, are slightly lower in comparison to the two other carriers (C1 and C2, the parents of P3) and all WT LCLs.

Figure 1.

Subtelomeric methylation in LCLs derived from ICF, carrier and WT individuals. (A) Methylation analysis of 12 subtelomeric regions by high-throughput sequencing of bisulfite-converted DNA. The heat map shows the methylation status across 7 amplicons that represent the 12 indicated subtelomeric regions, in WT LCLs (W1-W3), three DNMT3B-mutation carriers (C1-C3) and three ICF LCLs (P1-P4). In this heat map, each subtelomeric region consists of several lines, each line representing a specific CpG site analyzed in that subtelomeric region, in the 5′ to 3′ direction of the sequence. (B) Box plot presenting the distribution of the methylation percentage of all subtelomeric CpG sites for each of the ICF, carrier and WT samples as determined in (A).

Generation of chromosome end-specific telomere length data

We next prepared chromosomes from the ICF, carrier and WT LCLs described above. As demonstrated previously for ICF LCLs (53), we observed decondensation of pericentromeric regions, mainly in chromosome 1. We also detected dicentric chromosomes arising from telomere fusions at a rate of 4.3% (two fusions in 46 ICF metaphases), similar to the rate observed in this previous study (6% - 89 out of 1500 metaphases, described as telomere associations in the former study). However, surprisingly, we did not observe multiradial chromosomes involving chromosomes 1 and 16 that are typical of ICF cells. We then determined the relative telomere lengths of individual chromosome ends in 9–14 metaphase cells from each LCL, using telomere-FISH followed by spectral karyotyping (Fig. 2A–C). Analysis of the telomere-FISH images by the TFL-Telo software yielded arbitrary intensity values of Telomere Fluorescent Units (TFU) for both sister chromatids of each chromosome end, and these were assigned to the appropriate chromosome ends based on the spectral karyotyping. After validating a strong and significant correlation between sister chromatids (see ‘Methods and Materials’), we continued our analysis on the average TFU value of both chromatids of each chromosome end. Values of the homologous chromosome ends were not averaged but were analyzed independently (Supplementary Material, Table S1).

Figure 2.

Combined Telomere-FISH and Spectral Karyotyping. (A) A metaphase spread prepared from ICF P1 LCL hybridized to a telomeric PNA probe. (B) The same metaphase spread after spectral karyotyping. (C) Chromosomes from the above metaphase arranged in a karyotype table following analysis with the HiSKY software.

ICF telomere lengths are abnormally short and highly variant

Examination of telomere lengths by TRF analysis in LCLs generated from ICF patients [(37) and Supplementary Material, Fig. S1] demonstrated short telomere length in comparison to LCLs generated from their parents (possible to examine only in the case of P1 and P3), as well as to other WT samples. TRF analysis reflects the telomere length status of all telomeres in a large number of cells. However, the variation between different cells from a specific sample and between different telomeres within a specific cell cannot be determined from such an analysis. To this end, we analyzed the TFU dataset generated as described above, and determined the TFU median and distribution of all telomere values among the metaphases of each sample (Fig. 3A).

Figure 3.

TFU values distribute differently in ICF versus carrier and WT LCLs. (A) The TFU values of each telomere across all metaphases of each LCL were averaged and are presented as box-plots (the middle line depicts the median). All ICF LCLs show a lower median than the lowest WT sample, but a broader distribution of values in comparison to the other LCLs. (B) TFU values of all samples were compared by Anova (P<0.0001) followed by the Tukey‘s HSD test. The values and color of the cells reflect the extent of the mean difference between each pair of samples. The significant comparisons following adjustment for multiple comparisons (adjusted P-value<0.05) are marked by bold text. (C) The TFU values of all telomeres per LCL group [ICF, carriers (excluding C3) and WT] are presented as violin plots that combine the boxplot and the probability density plot. The ICF group exhibits both the lowest median value and the highest distribution of values, reflecting the high variance in this group. The carrier and WT groups show similar medians and distributions. Anova test and Tukey‘s comparisons between the three groups exhibit a significant difference between ICF and the other groups (P-value<0.0001) but no significant difference between the carrier and WT groups. (D) Each metaphase was plotted based on the mean TFU value (X axis) and the metaphase variance divided to the metaphase mean TFU (Y axis in log10 scale). The red, dark blue and light blue dots represent ICF, WT and carrier metaphases, respectively. The shapes represent the different samples, as depicted to the right.

As expected from the TRF results, we found that the ICF samples demonstrate notably low median TFU values. Interestingly, we found that the distribution of these values was very broad in each of the ICF LCL in comparison to the WT and carrier LCLs (Fig. 3A). Despite the low median TFU values, all ICF LCLs carried subgroups of telomeres with lengths comparable to those of the carrier and WT groups, and some of the ICF TFU values were even greater than those documented for the WT and carrier LCLs. This finding implies that not all ICF telomeres are shorter than WT telomeres, but rather there is a significant subgroup of telomeres that are extremely short, and this leads to the low median telomere lengths in ICF LCLs. We conclude from this data that not all telomeres shorten evenly in ICF syndrome cells.

Among the group of carrier LCLs, C1 and C2 exhibited medians and TFU value distributions very similar to those of the WT LCLs. In contrast, C3 even though showing low variation, revealed a very low median TFU (Fig. 3A). C3 LCL was generated from the mother of P1 who was reported by the Coriell Institute for Medical Research (GM08728) to be healthy. Despite this, C3 LCL shows a molecular phenotype that is skewed towards the ICF phenotype by several parameters. These include subtelomeric methylation levels (Fig. 1), TERRA levels (48) and profiles of RNA-seq and Chip-seq for H3K4me3, H3K9me3 and H3K27me3 (Maria Matarazzo and Miriam Gagliardi-personal communication). Intermediate molecular phenotypes in carriers of other recessively inherited diseases have been identified previously (54,55), and these may be dependent on the severity of the carrier‘s mutation, as the two other carriers (C1 and C2) presented normal subtelomeric methylation and telomere lengths (Fig. 3A and Supplementary Material, Fig. S1). In order to statistically evaluate the difference between the means of each sample, we used an Anova test (P-value < 0.0001) followed by a Tukey’s HSD test for multiple comparisons between samples with an adjusted P-value of <0.05 considered as significant (Fig. 3B). This analysis confirmed that telomere lengths in the C3 LCL differ significantly from those of the other two carriers and three WT LCLs. Based on this analysis, and on the molecular data mentioned above, we excluded sample C3 from the following analyses.

We next assembled the TFU values of all telomeres of all LCLs belonging to a certain group (ICF, Carriers, WT) (Fig. 3C). Comparison of the group means by an Anova test followed by a Tukey‘s HSD test demonstrated no significant difference between WT and carrier groups (P-value = 0.95), while both these groups differed significantly from the ICF group (P-value < 0.0001). In addition, examining the degree of similarity in the TFU variance values by the Levene’s test revealed that the ICF group differs significantly from both WT and carrier groups (P-value < 0.0001), however the WT and carrier groups were not significantly different from each other (P-value = 0.26).

The previous analyses indicated that the TFU values of all chromosome ends could segregate the ICF group from the carrier and WT groups. However, they did not indicate whether all metaphases in a particular group show similar TFU medians and distributions or whether certain ICF metaphases carried only short telomeres while others carried mostly long telomeres. In order to determine whether the metaphases within each group cluster together we studied the variability within each metaphase, compared to the mean TFU value of the metaphase (Fig. 3D). This analysis demonstrated that, remarkably, although several metaphases belonging to the ICF group had a mean TFU similar to metaphases from within the WT and carrier groups (Fig. 3A), they differed substantially by the variance within the metaphase, segregating the ICF group clearly from the carrier and WT groups. In addition, the metaphases of all ICF samples did not cluster according to the specific ICF LCLs that they belonged to, but rather as a distinct group (Fig. 3D). With regard to the groups of the carriers and WT LCLs, the metaphases belonging to these two groups clustered together.

ICF LCLs segregate from non-ICF LCLs based on ranking of telomere-specific lengths in individual metaphases

We next wanted to determine whether a classification model utilizing the telomere length data could infer which telomeres contribute most substantially to group separation. The previous analyses indicated that the control and carrier groups display similar data distribution (Fig. 3). Therefore, we first compared the classification accuracy by the random forest method of supervised classification and regression (56), when classifying the individual metaphases to either two (ICF/non-ICF) or three (ICF/carrier/WT) groups (Fig. 4A). We found that under separation to two groups the classification accuracy was 89%, while separating into three groups yielded only 72.4% accuracy. Accordingly, we combined the carrier and the WT groups and designated this new group as the ‘control group’.

Figure 4.

Classifying metaphases by random forest analysis. (A) Confusion matrices of metaphase classification. The metaphases were classified either by a two-group classification (ICF/non-ICF) or a three-group classification (WT/carrier/ICF). The values in the matrices are the number of analyzed metaphases. (B) Utilizing the two-group classification random forest model we extracted the variable importance list by the Gini importance index. The line was set at the point showing a shoulder in the Gini values. A high value indicates the high importance of a variable for an accurate classification.

We next proceeded, based on the Gini importance index, to utilize the random forest model to extract the variable importance list for classification. The Gini importance index indicates the discriminative value of a specific variable in obtaining an accurate classification. This was done by randomizing the TFU values of each individual telomere and measuring the degree by which this randomization affected the classification accuracy. Thus, this analysis ranks the telomeres based on the degree by which their relative length contributes to the differentiation between the ICF and control groups (Fig. 4B). This analysis indicated that telomeres 12p, 22q, 20q, 9p, 19q, 19p and 17q ranked among the seven most important telomeres whose lengths influenced the classification of a metaphase as belonging to the ICF vs. the control group.

Specific telomeres are prone to shortening in ICF

To further validate that a subgroup of telomeres in ICF syndrome exhibits preferentially short lengths, we compared the mean difference of each telomere between the ICF and control groups by the student t-test. For this analysis, we determined the mean TFU value for each telomere in all metaphases of a specific LCL sample and examined the mean of the four ICF LCLs compared to the five control LCLs (comprising the three WTs and the two carriers). After adjustment to multiple comparisons by the Benjamini & Hochberg false discovery rate method, 12 telomeres (19p, 19q, 17q, 12p, 15q, 22q, 3q, 2q, 20p, 1p, Xq/Yq and 17p) were found to exhibit significant differences between the two groups, all showing a lower mean value in the ICF group (Fig. 5, Supplementary Material, Table S2). Six out of the seven most important telomeres ranked by the random forest analysis (12p, 22q, 20q, 19q, 19p and 17q) (Fig. 4B) were among the telomeres that demonstrated a significant difference between both groups. Interestingly, nine telomeres had similar or higher mean values in the ICF versus the control group (15p, Xp/Yp, 5p, 14p, 2p, 4q, 21p, 3p and 7q) (Fig. 5), suggesting that certain telomeres may be protected from shortening in ICF syndrome cells. The ‘difference between means’ analysis also illustrated clearly the broad distribution of the mean value of the different telomeres in ICF LCLs in comparison to the relatively narrow distribution in the control LCLs.

Figure 5.

Telomere-specific differences between ICF and control groups. The mean difference between the ICF and control TFU values is presented for each telomere. The ordering of the telomeres from left to right is based on the degree of mean difference of telomere lengths calculated by mean length of the ICF group subtracted from the mean length of the control group, from the highest difference to the lowest difference. Dots and circles represent the mean TFU values of the ICF and control groups, respectively. The smoothing lines were fit by local polynomial regression fitting (ICF - full line, control – dashed line). Telomeres that show a significant difference following multiple comparison adjustments (Benjamini & Hochberg adjusted P-value<0.05) are marked by an asterisk. Full values appear in Supplementary Material, Table S2. The data for Xp and Xq includes data for Yp and Yq, respectively in the case of male-origin samples.

The majority of the most influential telomeres display subtelomeric sequence similarities

We next searched for common characteristics in the group of telomeres whose lengths contribute most to the segregation between ICF versus WT groups. All human telomeres end with the same hexameric repeat, however subtelomeres, while sharing many characteristics, differ from each other to varying degrees (57). We postulated that the variation in subtelomeric-CpG dinucleotide frequency might not be the sole factor that influences the propensity of a nearby telomere to shorten in DNMT3B-mutated cells. Therefore, we focused our analysis on the entire distal 2000 bp sequence of each chromosome end, with exception of the acrocentric p-arms for which no sequence is yet available. Sequence similarity between the available 41 subtelomeres was explored by the MAFFT sequence analysis tool (58), after which we organized the telomeres in a similarity matrix (Fig. 6) based on their order of importance in segregating the groups, determined by the random forest analysis (Fig. 4B). Interestingly, we found that seven out of the eight most influential telomeres that separate the ICF from the control group by the Gini index (12p, 22q, 20q, 9p, 19q, 19p and 2q) show a very high degree of sequence similarity based on the MAFFT tool analysis. This finding suggests that common sequence motifs in these subtelomeres are influenced by the DNMT3B loss-of-function in a manner that contributes to preferential telomere shortening in ICF syndrome.

Figure 6.

Sequence similarity matrix of the most distal 2 kb of 41 human subtelomeric regions. The sequence similarity matrix was created by the MAFFT algorithm and telomeres were ordered based on their importance in segregating the ICF group vs. the control group by the random forest analysis. Color and value in each cell indicate the percent of similarity as calculated by the MAFFT algorithm.

Because the human acrocentric p arm telomeres were excluded from this analysis, and since it has been shown that they are composed of typical subtelomeric repeat elements and are highly similar to each other (57), we asked whether these telomeres share any feature related to their telomere length. We focused on the five autosomal acrocentric chromosomes 13, 14, 15, 21 and 22, and excluded chromosome Y since our data analysis did not distinguish between the telomeres of chromosomes X and Y. Interestingly, the acrocentric p arm telomeres display significantly increased telomere lengths in comparison to all other autosomal telomeres in both groups. This difference appeared greater in the ICF group compared to the WT group (Supplementary Material, Fig. S2), however given our sample size, significance of this difference was limited (P = 0.14). These findings reinforce the notion that telomeres sharing similar subtelomeric sequences are likely to be similar in their relative lengths. In the case of the acrocentric p arm telomeres, the factors that influence telomere lengths may be influenced to some degree by the hypomethylated state of the ICF genome, either at subtelomeres or elsewhere, resulting in even longer telomeres in the patient cells.

Discussion

In this study, we describe the unique characteristics of telomere lengths in ICF syndrome cells mutated in the de novo DNA methyltransferase DNMT3B gene. Our analysis was carried out on LCLs that had grown for many passages in culture, thus reflecting the telomere length status following many generations of replication. Our current analysis, carried out by combined telomere-FISH and spectral karyotyping, confirms previous findings based on TRF analyses of significantly short telomere lengths in this syndrome (37,49,50). However, the ability to examine the relative length of each individual telomere by a molecular cytogenetic approach enabled us to reach deeper insights regarding the nature of telomere shortening in this disease.

In cells that express telomerase, telomere lengths reflect a balance point produced by opposing forces that lengthen and shorten telomeres (12). In primary fibroblasts of ICF patients that lack telomerase activity, the short telomere phenotype leads to premature senescence (49,50). However, in ICF iPSCs and in ICF fibroblasts expressing ectopic hTERT, the high levels of telomerase expression eclipse the short telomere phenotype (49,50). The significant shortening found in ICF LCLs suggests that the shortening forces, both the natural ones and those that stem from the molecular pathology of this disease, overcome the lengthening potential of the endogenous low levels of telomerase (37). In line with this, all ICF LCLs grown in culture arrested their growth following continuous culturing for several weeks in contrast to WT and carrier LCLs that continued to proliferate with no apparent impact on their replicative potential (results not shown).

The clear short telomere phenotype in ICF syndrome patient cells strongly suggests that this syndrome represents a new member in the category of secondary telomeropathies (23). Telomeropathies are a wide spectrum of Mendelian diseases in which the maintenance of telomeres is defective, telomeres are abnormally short and cells enter premature senescence. In primary telomeropathies, mutations occur in genes involved in telomere maintenance such as components of the Shelterin or Telomerase complexes. Patients with these disorders may start life with short telomeres, and consequently, senescence is triggered at an abnormally young age that varies substantially depending on the genetic defect. In this category of telomeropathies, the expectation is that the entire population of telomeres will demonstrate linear shortening at an accelerated pace (23).

The secondary telomeropathies are an emerging category of Mendelian diseases in which mutations occur mainly in genes responsible for DNA repair and replication, however the mechanism behind the telomeric phenotype in several of the secondary telomeropathies is yet unclear (35). It is well established that the telomeric chromatin is crucial for proper telomere function and maintenance (38). In this regard, subtelomeric methylation plays an important role in telomere function, and a mechanism by which hypomethylated subtelomeric DNA leads to DNA damage at telomeric regions in ICF syndrome was recently demonstrated (48). In secondary telomeropathies, telomere loss is predicted to occur stochastically such that not all telomeres are affected at a constant rate (23). Nevertheless, whether or not all telomeres are at an equal risk of shortening in this group of telomeropathies is yet unknown.

To this end, we asked whether all telomeres will be similarly affected in ICF syndrome, as has been shown during normal aging (59–61), or whether specific telomeres will demonstrate enhanced vulnerability to shorten due to the hypomethylated state of their adjacent subtelomeric region or other genomic regions affected by the mutated DNMT3B. The impact of subtelomeric methylation on telomere length homeostasis has been demonstrated in additional organisms such as mouse (62) and Arabidopsis (63), however the effect of methylation of a specific subtelomere on its neighboring telomere has not been studied previously. The combined telomere-FISH and spectral karyotyping approach enabled us to measure the relative lengths of each individual telomere in ICF and control LCLs. Two clear findings emerged from this analysis: a. telomere lengths in ICF cells are abnormally heterogeneous, and b. the shortening of telomeres in ICF syndrome LCLs is not random.

The development of peptide nucleic acid (PNA) probes for telomeric repeats greatly advanced the study of individual telomere lengths (64). Using such probes, telomere lengths of sister chromatids were found to be highly comparable (61,64). However, lengths of homologous telomeres were found to differ to a great degree (60,61,65) and the lengths of all telomeres of an individual revealed high intracellular heterogeneity (64,65). Indeed, we documented length variation in the control LCLs (Fig. 3A), however, variation in the ICF LCLs was significantly higher than in controls (Fig. 3A) and this variation supports the scenario that not all telomeres shortened at a constant rate throughout the replicative history of these cells. In the majority of the ICF metaphases a subset of telomeres were very long, extending even beyond the longest telomeres in control LCLs (Fig. 3A). This condition prevailed alongside a group of extremely short telomeres in the same metaphases. While this length variation was indeed abnormally wide, it was not as extensive as the variation documented in ALT cancer cells that undergo telomere sister chromatid exchange (T-SCE) [(66–68), reviewed in (69)]. In addition, we previously showed that in ICF LCLs, -fibroblasts and -iPSCs there is no evidence for elevated T-SCE (37,50). We therefore conclude that telomere loss occurred in ICF LCLs over many passages in a non-linear manner, with a subset of telomeres experiencing relatively large step-wise losses at each round of replication. The erosion of only a fraction of telomeres at any time point is also in concordance with our previous findings of DNA damage signals appearing only at a small group of chromosome ends in each metaphase spread of ICF LCLs (48). We therefore next asked whether or not telomere loss occurs preferentially at specific telomeres.

Analysis of the relative lengths of individual telomeres revealed that telomeres varied by the degree of mean length difference between the ICF and control groups (Fig. 5). Telomeres could be ranked according to this difference and certain telomeres could be designated as being more prone to erosion in ICF cells. Moreover, random forest analysis demonstrated that based on the relative length of a subgroup of telomeres, a metaphase could be assigned at very high confidence to either the ICF or the control group. This group of ‘influential telomeres’ only partially overlapped with different telomeres that were found in previous studies to be frequently the shortest telomeres in human cells (22,33,59,70), suggesting that disease-specific factors contribute to the typical pattern of telomere shortening in this syndrome.

The finding that not all telomeres are at similar risk to shorten in ICF syndrome raises the question regarding the mechanistic basis for the preferential shortening of a subset of telomeres in this disease. Most likely many factors affect in concert the length of a specific telomere both in health and disease states (12). The varying sequence characteristics documented at human subtelomeres, such as the distribution of repetitive sequences, have been proposed early on to play a role in telomere-length regulation in cis (71). The influence of subtelomeric sequence may be expected to stem from regions in close vicinity to the hexameric repeats, however long range interactions of telomeres with more internal chromosomal regions cannot be dismissed (57). Clearly, the study of telomere lengths in this syndrome calls to scrutinize subtelomere-specific characteristics relevant to DNA methylation. Indeed, a previous analysis of the distal 2 kb of the majority of human subtelomeres identified variations in the CpG density among chromosome ends (48). In addition, these subtelomeric regions vary with respect to characteristics such as GC density and GC-skewness (48). Moreover, the degree of GC skew played a role in DNA:RNA hybrid formation at telomeres. This in turn was implicated in generating DNA damage at telomeric regions that could lead to telomere loss (48). With respect to this last study, several but not all of the telomeres that demonstrated high DNA:RNA hybrids were amongst the shortest telomeres in the ICF LCLs. We predict therefore that telomere-hybrid formation in ICF cells contributes to the telomere-specific shortening documented in the present study, although not as a solitary factor, but rather in concert with other, yet to be identified, chromosome end-specific attributes that regulate telomere length in the scenario of compromised DNA methylation. The proposition that additional sequence characteristics of the most distal 2 kb of chromosome ends are instrumental in telomere length regulation is supported by our similarity matrix analysis (Fig. 6) that reveals a strong degree of similarity among seven out of the eight most ‘important telomeres’, based on the random forest analysis. We predict that subtelomeric sequence characteristics influence in addition to hybrid formation, many additional functions such as chromatin configuration, replication timing and nuclear positioning (57,72–74). Such characteristics may in turn affect both the vulnerability of telomeres to undergo damage, and the capacity to repair this damage once produced.

Since our study was carried out on LCLs expressing telomerase, we cannot exclude that in addition to the possible mechanisms described above, the pattern of relative telomere lengths in ICF LCLs is a consequence of differential recruitment of telomerase to various telomeres. While telomere-specific differential recruitment of telomerase has not been reported, it has been shown in several organisms that telomerase is recruited preferentially to the shortest telomeres (15,16). In budding and fission yeast, this occurs via the TERRA transcribed off these telomeres (17,18). Since short telomeres were found to transcribe higher levels of TERRA, this could explain the preferential elongation of short telomeres (17,18). Conversely, TERRA was recently shown to inhibit telomerase in mammalian cells (8). Taking in consideration TERRA’s involvement in telomerase recruitment, and the fact that TERRA expression is enhanced unequally at various telomeres in ICF LCLs, the capability of telomerase in ICF cells to elongate telomeres could vary at different chromosome ends, such that highly-transcribing TERRA telomeres could recruit telomerase at higher efficiency. On the other hand, high TERRA levels result in DNA:RNA hybrids at telomeres which are postulated to lead to DNA shortening through DNA damage (48). Analysis of telomere lengths in ICF cells that do not express telomerase could enhance our understanding of the mechanism/s responsible for setting the characteristic pattern of telomere lengths in ICF cells and lead to a better understanding of the contribution of TERRA-telomerase interactions in telomere length homeostasis.

The analysis of the various LCLs in this study indicates that DNMT3B mutation-carriers might exhibit a partial molecular phenotype, as seen in the case of carrier C3. This example joins at least two additional examples of recessive diseases in which carriers have been found to exhibit partial phenotypes. Thus, this calls for awareness to the possibility that carriers of recessive diseases may exhibit an abnormal phenotype later in life, albeit at less severity in comparison to the patients, due to a partial loss of function of the mutated protein. In the case of DNMT3B mutations, this might place the carriers at an elevated risk to conditions such as earlier onset of aging or propensity to develop malignancies.

The finding that a common pattern of relative telomere length is shared among the various ICF patient LCLs, each carrying different mutations in DNMT3B, supports a significant contribution of subtelomeric DNA methylation to the regulation of the individual telomere length. The demonstration that an epigenetic modification can play a role in telomere-specific length regulation calls for future investigations regarding the contribution of other telomere-specific epigenetic modifications to length regulation of individual telomeres. In addition, our findings are the first to show that a telomeropathy may display an uneven but consistent pattern of shortening at different telomeres. We believe that as such, similar analyses carried out on additional primary and secondary telomeropathies will contribute to a better understanding of the various factors that regulate telomere-specific length homeostasis.

Materials and Methods

Cell lines and cell culture

Four ICF LCLs designated as P1, P2, P3 and P4 (details of all LCLs appear in Table 1) were grown in RPMI supplemented with 20% FCS, 2 mM glutamine, 100 U/ml penicillin and 100 μg/ml streptomycin. Control WT LCLs designated W1, W2 and W3 and DNMT3B mutation carriers designated C1, C2 and C3 (Table 1) were grown in similar media as above, supplemented with 15% FCS.

Bisulfite sequencing and analysis

Bisulfite-converted DNA was sequenced using the Ion Torrent PGM high-throughput technology (75). Genomic DNA (1 μg) was bisulfite-converted with the Methylamp DNA Modification kit (EPIGENTEK, NY, USA) according to the manufacturer‘s instructions. After bisulfite conversion, DNA was amplified using FastStart Taq DNA Polymerase (Roche) as follows: 95 °C for 5 min; 4 cycles of 95 °C for 1 min, 53 °C for 3 min, 72 °C for 3 min; 2 cycles of 95 °C for 30 s, 55 °C for 45 s, 72 °C for 45 s; 40 cycles of 95 °C for 30 s, 72 °C for 1.5 min; 72 °C for 10 min. Barcode and adaptor containing-primers for amplification were described previously (50). PCR products were cleaned with the QIAquick PCR purification kit (#28104, Qiagen) and all amplicons were quantified, diluted to 10pM, pooled into one tube and subjected to emulsion PCR (E-PCR) on Ion Sphere Particles (ISPs; Life Technologies) using the Ion PGM™ Sequencing 400 Kit on the OneTouch system (Life Technologies). Reaction efficiency was estimated by calculating the percentage of DNA-containing ISPs compared to the background of blank ISPs using the Qubit Ion Sphere quality control kit (Life Technologies). Enriched ISPs were subjected to sequencing on the Ion 314 Chip using the Ion PGM™ Hi-Q™ Sequencing Kit. The multiple reads generated for each amplicon were aligned to the target sequences, after which FASTQ files were generated and the percent of methylation for each CpG site was calculated, all using Bismark software (Babraham Institute website) (76).

Telomere FISH

Cultured cells were treated with 0.2 μg/ml colcemid for 30 minutes. Cells were then swollen with 67 mM KCl for 7 min at room temperature followed by methanol/acetic acid (3:1) fixation. The fixative was replaced three times. Fixed cells were stored in 5 ml of fixative at −20 °C. Cells were dropped on SuperFrost Plus slides (Thermo Science) and slides were air-dried overnight. Hybridization was performed with a PNA oligo probe of (CCCTAA)₃ conjugated to Cy3 according to the Q-FISH protocol (77). Hybridized slides were stored in a light-protected storage box at 4 °C until images were captured. Slides were visualized on a BX51 Olympus microscope and images were captured using dedicated acquisition software for capturing FISH signals (GenASIs FISHView) from Applied Spectral Imaging (ASI). Metaphases were captured per each LCL within one week following hybridization with the telomere FISH probe. For each captured metaphase the slide coordinates were recorded. Analysis of the hybridization signal intensities following telomere-FISH was carried out using the TFL-Telo software (78). This analysis produces arbitrary fluorescence intensity values for each chromosome end (see below – ‘Data Transformation’).

Spectral karyotyping (SKY)

Following image capture as described above, slides were washed in 2XSSC and processed for a second round of hybridization using the Spectral Karyotyping (SKYPaint) kit from Applied Spectral Imaging (ASI). Hybridization was carried out according to the manufacturer's protocol with slight modifications, as described (15).

The metaphases captured after the Telomere-FISH, were recaptured following SKY hybridization using the GenASIs Hyperspectral system equipped with the HiSKY software (Applied Spectral Imaging, Yokneam, Israel). This system enables measurement of the full visible light spectrum at each pixel of the image by using a Sagnac interferometer. A classification algorithm was used to differentiate between different spectra in the image and to assign pseudo colors to all the pixels that share similar spectral characteristics. The DAPI images were captured separately and inverted to give a G-banding like pattern. The chromosomes were then sorted automatically into a karyotype table. A trained cytogeneticist verified the karyotyping, as well as the correct designation of the p and q arms in the case of metacentric chromosomes.

Data transformation

The telomere intensity values of each metaphase were assigned to the correct chromosome arm by integrating manually the data of the telomere-FISH with the SKY-karyotype. Two individuals verified independently the assignment of the telomere-FISH signals to the specific chromosome arms.

The TFL-Telo software produces four fluorescent intensity values for each metaphase; P1 and P2 represent the two p-arm sister chromatids, and Q1 and Q2 represent the two q-arm sister chromatids. Similar to previous reports (61,64) the values of sister chromatids were found to be highly correlated (by Pearson correlation coefficient values), therefore we averaged the sister chromatid data for each chromosome end. In the rare instances of a missing telomere fluorescence value due to technical reasons, it was assigned a value that corresponded to the average of the corresponding telomere in all other metaphases of the same sample. This ensured that such a telomere would not be recognized erroneously as one lacking a signal. To note, it was necessary to include these data points to preserve the data matrix dimensions and facilitate error-detection. Chromosomes X and Y in male samples (W1, C1, P2 and P3) were treated as homologous and assigned to the features Xp/Yp and Xq/Yq. After the aforementioned substitutions, homologous telomeres were randomly designated identities of ‘1’ or ‘2’. Each row of data contains the following features: data key (WT, Carrier or Patient, sample-ID, metaphase-ID, homolog-ID, sex of sample) followed by the arbitrary fluorescence values, as determined by the TFL-Telo software, for the 46 telomeres (p and q arm values for 22 somatic chromosomes and one X/Y chromosome).

Statistical analysis

All statistical analyses were done using R software packages ′reshape2′, ′car′, ′randomForest′, ‘plyr’ and ′MASS′. For graphics ′ggplot2′, ′cowplot′ and ′RColorBrewer′ packages were used. For each parametric statistical test, we ascertained the test assumptions. For the normality assumption we used the Shapiro-Wilks test and histograms of each group. For the independency of observation assumption, all statistical analyses were done on the average TFU of each telomere across all metaphases of a specific sample, thus for each sample for each telomere end there is only a single observation. The t-test was conducted under the assumption of non-equal variances.

Supplementary Material

Supplementary Material is available at HMG online.

Acknowledgements

We are very grateful to Claire Francastel and Guillaume Velasco for providing us with the majority of patient and carrier LCLs. We thank Eyal Bergmann and Omer Schwartzman for help with data analysis and Maria Matarazzo and Miriam Gagliardi for sharing with us unpublished data. We thank Daniel Kornitzer and Aya Tzur-Gilat for comments on the manuscript. S. Sagie is grateful to the Azrieli Foundation for the award of an Azrieli Fellowship.

Conflict of Interest statement.None declared.

Funding

The Legacy Heritage Bio-Medical Program of the Israel Science Foundation [657/15 to S. Selig].

Author notes

^†Shira Sagie, http://orcid.org/0000-0001-7234-0256

^‡Shir Toubiana, http://orcid.org/0000-0003-3001-5281