Genome-wide identification of expression quantitative trait loci for human telomerase

Supplemental Digital Content is available in the text


Introduction
Telomeres play key roles in human genome stability through replenishing their G-rich sequences by telomerases. [1] A variety of human cancer cells show high expression levels of telomerase with short telomeres whereas most human somatic cells have low expression levels of telomerase. [2] The telomerases consist of human telomerase reverse transcriptases (hTERT) that regulate telomerase activity and telomerase RNA components (TERC) that are used as template RNAs for lengthening telomeres. [2] It has been reported that mRNA expression level of hTERT is highly correlated with telomerase activity in tumor cells. [3] TERC are targeted by imetelstat sodium (GRN163L), an inhibitor of telomerase activity in cancer therapy. [4] Nevertheless, genetic factors that regulate the variability of gene expression are not well understood. In particular, knowledge on genetic factors is limited to genetic variability in or near the genes encoding hTERT and TERC. Some nucleotide variants (rs2736108 upstream of hTERT, rs7705526 in intron 2 of hTERT, and rs12696304 downstream of TERC) were associated with telomere length in leukocytes. [5,6] Some intronic variants (rs10069690, rs2242652, and rs7725218) of hTERT have been reported to have association with hTERT expression in prostate cancers. [7] This study aimed to identify nucleotide variants associated with mRNA expression of hTERT and TERC through genome-wide analysis.

Subjects
Expression data of the genes encoding hTERT and TERC in lymphoblastoid cell lines generated from the Geuvadis RNAsequencing project [8] were used to identify expression quantitative trait loci (eQTLs). Cell lines were derived from 373 Europeans of the following 4 populations: Utah residents with northern and western ancestry (n = 91), Finns (n = 95), British (n = 94), and Toscani (n = 93). [8] We excluded Yoruba population from the project to avoid false positive associations produced by heterogeneous genetic background. Gene expression was calculated as the sum of reads per kilobase per million mapped reads (RPKM) for all transcripts of each gene in each individual. [8] Their corresponding genotypic data were obtained from the 1000 Genomes Project (http://www.1000ge nomes.org/). Genotypes with minor allele frequency <5% or with missing rate of >5% were removed. After the quality control, genotypes of 5,851,914 SNPs were used for final analysis. Ethical approval was not necessary because we dealt with publically available data.  Chromosome: chromosomal position (bp) from GRCh38.p2. † The "upstream" in the parentheses indicates that the variant is located within 2 kb upstream of gene. The "ncRNA" indicates that the variant is located in a noncoding RNA gene. The "intronic variant" indicates that the variant is located in introns of LOC105370092, but the intronic number is uncertain because the gene is uncharacterized. ‡ Minor allele/major allele.
x SNP in bold indicates representative variant with the most significance at each eQTL.

Statistical methods
Linear regression analysis was performed to discover eQTLs of hTERT and TERC. Multiple testing was employed with significance threshold value of 5 Â 10 -8 . All the statistical analyses were conducted using PLINK. [9] Linkage disequilibrium (LD) blocks at association signals were constructed using HaploView. [10] The identified eQTLs were analyzed to determine if they were located in transcription factor binding sites using ChIP-seq data from the regulomeDB. [11] 2.3. Transcriptome-wide association analysis of eQTLs for human telomerase reverse transcriptase Further associations between 10,518 genes and SNPs identified to be associated with the expression of hTERT were examined. For transcriptome-wide association analysis, the significance threshold value was set at 4.75 Â 10 -6 (= 0.05 divided by the total number of 10,518 genes). Functional enrichment analysis was conducted with the identified genes to examine their functional relevance using the DAVID functional annotation tool. [12] 3. Results Genome-wide association analysis revealed 37 eQTLs associated with mRNA expression of hTERT (P < 5 Â 10 -8 ; Fig. 1, Table 1). However, no significant eQTL was identified to be associated with mRNA expression of TERC (P > 5 Â 10 -8 ; Fig. 1). Some of the identified eQTLs located in chromosomes 6 and 12 turned out to be in strong linkage, and linkage disequilibrium blocks were constructed in Fig. 2. As a result, we found 6 association signals (rs224514, rs112953754, rs17755753, rs35070061, rs2636908, and rs187444335) from the eQTL analysis (Table 1). Among them, the intronic nucleotide variant rs224514 of the gene encoding sedoheptulokinase (SHPK) showed the most significant association signal (P = 1.50 Â 10 -10 ). Another intragenic association signal was located in intron 1 of R-spondin-3 (RSPO3) gene (rs17755753, P = 3.91 Â 10 -9 ). Transcriptome-wide association analysis with the 6 eQTLs identified to have association with hTERT expression showed significant associations with mRNA expressions of additional 29 genes (P < 4.75 Â 10 -6 ; Table 2). No gene other than hTERT was identified to have significant association with rs35070061 (P > 4.75 Â 10 -6 ).

Discussion
The current study identified 6 genetic association signals that might regulate the expression of the gene encoding hTERT and explain partial variability of its expression. In particular, a signal located in the gene of RSPO3 was found. The association of RSPO3 with the expression of hTERT could be supported by previous studies. It has been shown that the R-spondin family can disrupt the inhibition of LRPs caused by DKK1, and thus LRPs activate the Wnt signaling pathway. [13,14] The RSPO3 can enhance angiogenesis and proliferation of human endothelial cells as a Wnt signaling regulator. [15] It has been reported that b-catenin, a critical intracellular signal transducer in the Wnt signaling pathway, can regulate the mRNA expression of hTERT. [16][17][18] Thus, RSPO3 can induce hTERT through the Wnt signaling pathway. One of the 6 association signals covered 26 nucleotide variants in strong linkage within LOC105370092. Some of these variants were corresponding to a transcription factor binding sites uncovered by ChIP-Seq (RegulomeDB; Supplementary Table 1, http://links.lww.com/MD/B366), suggesting that these variants might regulate the expression of the uncharacterized ncRNA gene. Significant associations of these variants with the expression of CTNND2 were found (P < 4.75 Â 10 -6 ; Table 2), indicating that the ncRNA might influence the expression of hTERT and CTNND2. The RegulomeDB showed that rs113285167 was found in lymphoblastoid cells as binding sites of transcription factor 12 (TCF12) and CCCTC-binding factor (CTCF) (Supplementary Table 1, http://links.lww.com/ MD/B366), which have opposite functions as a transcription factor and a repressor, respectively. The functionally contrasting proteins could bind the locus even in the same cell line (GM12878; ENCODE http://genome.ucsc.edu/ENCODE), but relative binding affinity remained unknown in these studies, especially by alleles. Further research on specific underlying mechanism on their interaction with binding sites is required to understand how gene expression is affected by alleles.
Transcriptome-wide association analysis revealed that the 6 signals associated with hTERT expression were further associated with the expression of many other genes. The most significant association was observed between prickle planar cell polarity protein 2 (PRICKLE2) and rs2636908. Since the PRICKLE2 gene and this SNP were both located in chromosome 3, a regulatory element for expression of PRICKLE2 might include this SNP or other SNPs strongly linked to it. The complex of PRICKLE2 with Vangl2 can regulate the Wnt/planar cell polarity (PCP) pathway that changes cytoskeleton. [19] Genes identified with eQTL of rs2636908 were enriched with cytoskeleton (Supplementary Table 2, http://links.lww.com/MD/B366). Acute withdrawals of TERT in mouse have triggered a rapid change in the expression of genes with the functions in the cytoskeleton, [20] suggesting that hTERT might be involved in cytoskeletal changes associated with the PCP pathway.
The current study dealt with expression of telomerase genes in lymphoblastoid cell lines transformed by Epstein-Barr virus (EBV). [8] The lymphoblastoid cell lines may be suitable for identifying eQTLs of telomerase genes since they demonstrated strong telomerase activity to be immortalized. [21][22][23] In EBVinfected B cell, expression of ectopic hTERT simultaneously increased with expression of basic leucine zipper ATF-like transcription factor (BATF) which maintains EBV latent status, and hTERT silenced by shRNA can induce B-cell death through transition into lytic cycle of EBV. [24] Thus, lymphoblastoid cells transformed by EBV can be pertinent to hTERT expression studies. Also, we suspect that the Wnt signaling pathway can be instrumental in EBV-transformed lymphoblastoid cells. This is because b-catenins are expressed in lymphoblastoid cells, and moreover, EBV-transformed lymphoblastoid cells increased the expression of b-catenin more than EBV-negative cells. [25][26][27][28] Any eQTLs for TERC, another important component of human telomerase, were not found in this study (P > 5 Â 10 -8 ). However, many false negatives are suspected because we employed a conservative significance threshold. Some "suggestive" eQTLs for TERC are presented in Fig. 1 as potential false negatives (10 -5 > P > 5 Â 10 -8 ).
The current GWAS identified 6 novel eQTLs for hTERT. Some of these eQTLs were involved in the Wnt signaling pathway critical to the production of tumors. Functional studies are needed in order to understand the underlying mechanisms of the Wnt signaling pathway influenced by hTERT. This will provide some evidence to unveil the reason why hTERT expression is associated with tumor.