The Telomere-Telomerase System Is Detrimental to Health at High-Altitude

The hypobaric-hypoxia environment at high-altitude (HA, >2500 m) may influence DNA damage due to the production of reactive molecular species and high UV radiation. The telomere system, vital to chromosomal integrity and cellular viability, is prone to oxidative damages contributing to the severity of high-altitude disorders such as high-altitude pulmonary edema (HAPE). However, at the same time, it is suggested to sustain physical performance. This case-control study, comprising 210 HAPE-free (HAPE-f) sojourners, 183 HAPE-patients (HAPE-p) and 200 healthy highland natives (HLs) residing at ~3500 m, investigated telomere length, telomerase activity, and oxidative stress biomarkers. Fluidigm SNP genotyping screened 65 single nucleotide polymorphisms (SNPs) in 11 telomere-maintaining genes. Significance was attained at p ≤ 0.05 after adjusting for confounders and correction for multiple comparisons. Shorter telomere length, decreased telomerase activity and increased oxidative stress were observed in HAPE patients; contrarily, longer telomere length and elevated telomerase activity were observed in healthy HA natives compared to HAPE-f. Four SNPs and three haplotypes are associated with HAPE, whereas eight SNPs and nine haplotypes are associated with HA adaptation. Various gene-gene interactions and correlations between/among clinical parameters and biomarkers suggested the presence of a complex interplay underlining HAPE and HA adaptation physiology. A distinctive contribution of the telomere-telomerase system contributing to HA physiology is evident in this study. A normal telomere system may be advantageous in endurance training.


Introduction
Hypobaric hypoxia at high-altitude (HA, >2500 m) triggers oxidative stress that activates or even exaggerates several physiological pathways, disturbing the normal physiology of the body [1][2][3][4][5][6][7]. Elevated reactive oxygen species (ROS) and reactive nitrogen species are the basis of oxidative stress and may influence hypoxia-inducible factor-1 (HIF-1) mediated pulmonary leakage and also damage the sodium transport system in alveoli epithelium promoting fluid retention [8][9][10][11] to exaggerate HA-related disorders such as high-altitude pulmonary edema (HAPE [MIM: 178400]) [1,12,13]. ROS damage DNA, which is further aggravated when exposed to the high UV radiation that prevails at HA. DNA damage may occur as early as the first day of HA exposure [14]. The telomeric DNA, due to its richness in guanine, is more prone to oxidative damage than other regions of the DNA [15,16]. Since the telomere-telomerase system protects the chromosomes, shorter telomeres may decrease the viability of cells, leading to cellular dysfunction and disturbed vascular homeostasis. Telomere length has evolved as a potential candidate as a disease risk marker, with various studies reporting the association between shorter telomere length in leukocytes and agerelated or stress-related diseases such as Alzheimer's disease, diabetes, cardiovascular and pulmonary diseases, myocardial infarction, hypertension, chronic obstructive pulmonary disease (COPD) and others [17][18][19][20][21]. It clearly suggests that a normal telomere system benefits healthy people in several ways [22][23][24][25][26]. Exercise training has been demonstrated to maintain active telomerase and keep telomere length under control; thus, a longer telomere will be beneficial at altitude [23][24][25][26]. As a result, altitude training in a hypoxic environment is known to increase athletes' exercise capacity and performance at sea level and HA [22,26,27]. However, the telomere system has yet to be investigated in detail at HA to substantiate the concept.
The telomere-telomerase system is comprised of telomeres, the shelterin complex and telomerase. It determines chromosome integrity and cell viability [28][29][30][31][32]. The shelterin complex masks the telomeric DNA from attrition-causing factors and regulates telomere extension and telomeric DNA replication [28][29][30][31]. Telomerase, a ribonucleoprotein, catalyzes telomere extension by adding six nucleotide repeats to telomeric ends [32]. Several components of the telomere-telomerase system, such as TERF1 interacting protein 2 (TINF2 [MIM: 604319]), telomerase reverse transcriptase (TERT [MIM: 187270]) and heat shock protein 90 (Hsp90/HSP90AA1 [MIM: 140571]) regulate ROS production and oxidative stress [33][34][35][36][37]. We propose that functional impairment of shelterins and telomerase because of HA environmental stress and genetic variations may cause telomere attrition and promote HAPE and other HA disorders ( Figure 1A). The present study hence comprehensively investigated the physiological and genetic significance of the telomere-telomerase system in HAPE and HA adaptation under the hypobaric hypoxic environment of HA. Physiological relevance was investigated by evaluating telomere length and telomerase activity along with the oxidative stress biomarkers, 8-iso prostaglandin F2α (8-isoPGF2α) and total antioxidant activity, and expression of the genes ( Figure 1B). The genetic study included eleven genes, including shelterin genes and telomere-associated genes (Supplementary Table S2). The telomere-telomerase system is influenced by oxidative stress under the hypobaric hypoxic environment of high-altitude. In a sequential process, they play a crucial role in HAPE pathophysiology and HA adaptation. (B) Study design. Study subjects were recruited, and their sampling was conducted at high-altitude (~3500 m above sea level). Socio-demographic and clinical data were collected at the time of recruitment. Plasma biomarkers namely, 8-isoPGF2α and total antioxidants, telomere length and telomerase activity were estimated, followed by SNP genotyping, gene expression and in silico analyses to evaluate the site-specific protein structure variations. All the data were analyzed using relevant statistical and bioinformatic tools and are described in the Methods section. ROS, reactive oxygen species; HAPE, high-altitude pulmonary edema; HAPE-p, HAPE patients; HAPE-f, HAPE-free healthy controls; HLs, high-land natives; SNPs, single nucleotide polymorphisms.

Study Subjects
The present cross-sectional case-control study was carried out in three well-characterized study groups, i.e., HAPE-free healthy controls (HAPE-f; n = 210), HAPE-patients (HAPE-p; n = 183), and high-land natives (HLs; n = 200). Samples were collected through the general outpatient clinic of SNM Hospital, Leh at~3500 m after obtaining signed informed consent. The inclusion/exclusion criteria of subjects and data collection methods are detailed in the online Supplementary Methods, and the clinical characteristics of the study groups are summarized in Supplementary Table S1.

Ethics Approval
Ethical committees of the Council of Scientific and Industrial Research-Institute of Genomics and Integrative Biology (CSIR-IGIB), Delhi and Sonam Norboo Memorial (SNM) Hospital, Leh, Ladakh, Jammu and Kashmir approved the investigation. All procedures were performed in accordance with the Declaration of Helsinki.

Estimation of Relative Telomere Length
qRT-PCR was run with MESA BLUE qPCR Master Mix Plus for SYBR ® Assay No ROX (Eurogentec, Seraing, Belgium) on LightCycler 480 Real-Time PCR System (Roche Molecular Diagnostics, Pleasanton, CA, USA). The relative telomere length, T/S ratio, is the ratio of the telomere repeat copy number (T) to the copy number of a single copy gene (S) [38].

Estimation of Telomerase Activity and Oxidative Stress Markers
Telomerase activity was estimated by Quantitative Telomerase Detection Kit (Allied Biotech Inc., Ijamsville, MD, USA). The 8-isoPGF2α level was estimated using an ELISA kit (Enzo Life Sciences Inc., Farmingdale, New York, NY, USA) and total antioxidant activity spectrophotometrically (Cayman Chemical, Ann Arbor, MI, USA).

Genotyping of Single Nucleotide Polymorphisms
Single nucleotide polymorphisms (SNPs) were selected based on location, clinical and functional relevance, and thorough literature survey. HapMap and Indian Consortium data were referred for further selection and inspection for the presence of any tagged SNPs. The sixty-five selected SNPs included 6 TERF1, 5 TERF2, 7 POT1, 2 TINF2, 4 ACD, 6 TERF2IP, 16 TERT, 2 TERC, 9 TEP1, 5 HSP90AA1 and 3 PTGES3 (Supplementary Table S2). Fluidigm Dynamic Array Integrated Fluidic Circuit technique on a Biomark HD System (Fluidigm Corporation, South San Francisco, CA, USA) performed the genotyping.

Genetic Analyses
The Hardy-Weinberg equilibrium (HWE) for SNPs was checked using a χ 2 goodnessof-fit test using EpiInfo™ ver. 6 (Center for Disease Control, Atlanta, GA, USA). Haploview 4.2 checked the Linkage-disequilibrium (LD) and the tagging efficiency of SNPs. Two or more SNPs were considered in LD with D ≥ 0.9 and r 2 ≥ 0.9. Multivariate logistic regression analysis established the association of SNPs. The power of the study was calculated using the online tool OSSE-An Online Sample Size Estimator (http://osse.bii.astar.edu.sg/calculation2.php (accessed on 14 May 2017)).
To avoid false-positive results due to ethnicity difference between HLs and HAPEf/HAPE-p and to ascertain SNPs associating with HA adaptation, HLs were compared with two more populations, i.e., Han Chinese (CHB, Han Chinese in Beijing, China and CHS, Southern Han Chinese, China) and Japanese (JPT, Japanese in Tokyo, Japan). Genotype distribution data for CHB, CHS and JPT were procured from the 1000 Genomes Phase 3 project (www.ensembl.org (accessed on 2 June 2017)) [39]. Haplotypes were deduced from the genotypes of genes located at five chromosomes by the PHASE v. 2.1.1 program. The SNPs at the forward strand were recognized according to the contig position in the haplotypes (Supplementary Figure S1). An OR > 1 depicts risk, and OR < 1 depicts protection or adaptation.

Gene-Gene Interactions in HAPE and Adaptation
Gene-gene interactions for respective associated SNPs were evaluated by MDR 2.0 beta 8.4. The significant best models were selected based on higher scores of testing accuracy (TA ≤ 0.55) and cross-validation consistency (CVC = 9/10 or 10/10). A p-value ≤ 0.05 was considered statistically significant.

DNA-Protein Docking
Transcription factors (TFs) were identified using the dbSNP tool and PASTA program. The Hex docking program was used to dock TFs on the promoter region and to determine the binding affinities. PyMOL was used to visualize the docked structure. The 3D structure of TFs was retrieved from Protein Data Bank (PDB). The docking analysis was further cross-checked with the HADDOCK and CLusPRO.

Protein-Protein Interactions
STRING ver. 11.0 identified the protein-protein interactions of the respective 10 genes. TERC was excluded because it does not code a protein but is a constituent of the RNA component of telomerase.

Expression Analysis of Telomere Maintaining Genes
Differential gene expression was performed using qRT-PCR as detailed in the Supplementary Methods. The ∆∆Ct calculated the relative change in transcript with 18S rRNA as the reference gene. Primer sequences are presented in Supplementary Table S3.

Statistical Analyses
Statistical tests were applied using Statistical Package for Social Sciences version 16.0 (SPSS 16.0 [IBM Inc., Armonk, NY, USA]), EpiInfo™ ver. 6 (Center for Disease Control, Atlanta, GA, USA) or Simple Interactive Statistical Analysis (SISA) online tool (http:// www.quantitativeskills.com/sisa/ (accessed on 25 September 2017)). Unpaired Student's t-tests compared the differences in baseline characteristics, demographic features, and gene expression between the two groups. Biomarkers were expressed as mean ± standard error (SE) and SISA calculated the statistical power. A general linear model (GLM) was applied to calculate age, gender and BMI-adjusted p-values. Multivariate logistic regression analyses of single locus and haplotype/multilocus yielded p-value, chi-square (χ 2 ), odds ratio (OR) and 95% confidence interval (CI). The significance was attained after FDR correction (BenjaminiHochberg.xlsx calculator). Pearson's Correlation was used to ascertain the correlations between clinical parameters and biomarkers and within biomarkers. Alleles and genotypes of SNPs associated with HAPE susceptibility and HA adaptation were correlated with clinical parameters and biomarkers. Bonferroni's correction was applied for multiple comparisons. A p-value of ≤0.05 was considered statistically significant after adjusting with confounders.

Low SaO 2 Indicates HAPE
HAPE-p had lower SaO 2 levels than HAPE-f. The HLs of Tibeto-Burman origin had SaO 2 levels close to HAPE-f but higher than HAPE-p ( Figure 2A). In addition, HAPEp had increased mean arterial pressure (MAP) compared to HLs ( Figure 2B) and was non-significantly higher than HAPE-f. The increased stressor and decreased antioxidant complement the two clinical parameters and the shortened telomere length and the telomerase activity in HAPE patients compared to the reverse trend in the two control groups suggesting their contribution to HAPE pathophysiology. The data are represented as mean ± SD. SPSS 16.0 was used to obtain p-values after adjusting with age, gender, and BMI by performing UNIANOVA. A p-value of ≤ 0.05 was considered statistically significant. * represents p ≤ 0.05, ** represents p ≤ 0.01, *** represents p ≤ 0.001 and **** represents p ≤ 0.0001. HAPE-f, HAPE-free healthy controls; HAPE-p, HAPE-patients; HLs, High-land natives; SaO 2 , arterial oxygen saturation; MAP, mean arterial pressure; T/S, ratio between telomere repeat copy number (T) and the copy number of single copy gene (S); 8-iso PGF2α, 8-isoprostaglandin F2α; p, p-value.

Aberrant Telomere-Telomerase System Inclined with HAPE Susceptibility
Among the three study groups, HAPE-p had the shortest relative telomere length ( Figure 2C), which complies with the lower telomerase activity ( Figure 2D). In contrast, longer telomere length and higher telomerase activity were observed in HLs. The relative telomere length was the longest in the HLs compared to the two sojourn groups, i.e., HAPEp and HAPE-f ( Figure 2C). However, the difference between the two healthy groups was modest. The telomerase activity was in agreement with longer telomere length in HLs compared to HAPE-p and HAPE-f ( Figure 2D).

The Oxidative System Shows Stronger Influence
The levels of 8-isoPGF2α were elevated ( Figure 2E), whereas the total antioxidant activity was slightly lower ( Figure 2F) in HAPE-p compared to HAPE-f. The HLs present an interesting status of 8-isoPGF2α and total antioxidant activity. Similar to HAPE-p, the 8-isoPGF2α levels were higher in HLs than HAPE-f ( Figure 2E) and were supported by thẽ 20% increase in the total antioxidant activity than in the two sojourn groups ( Figure 2F). This suggests that HLs are under continuous stress due to an active oxidative process.

Genetic Association with HAPE and HA Adaptation
Sixty-five SNPs with six genes of telomere-associated shelterins, namely TERF1, TERF2, POT1, TINF2, ACD and TERF2IP (also known as RAP1) and five genes of telomerase system, namely TERT, TERC, TEP1, HSP90AA1 and PTGES3 were genotyped. Sixteen SNPs with <95% call rate or being non-polymorphic and six SNPs deviating from HWE in HAPE-f were excluded from the study. The remaining 43 SNPs were not consistently in LD or tagged among each other (D ≥ 0.90, r 2 ≥ 0.90) in the three study groups (Supplementary Figure S2A-G). These were screened and analyzed for genotype and allele distributions and are detailed in Supplementary Tables S4-S13; however, genotypic and allelic distributions of the significant SNPs are presented in Table 1. Subtables B in Supplementary Tables S4-S12 represent comparisons between the HLs against the Han Chinese (CHB+CHS) and Japanese (JPT) populations (1000Genomes).
Three genes with respective alleles such as TERF2IP rs59297469T, TERC rs2293607G and TEP1 rs2228036T, associated with HAPE risk (Table 1 and Supplementary Tables S8A, S10A and S11A), while ACD rs72556537A associated with protection against HAPE (Table 1 and Supplementary Table S7A). Eight alleles TERF1 rs2975843T and rs10099824G, TERF2 rs3785073G and rs153058C, POT1 rs7794637C, ACD rs6979A and TEP1 rs2184282G and rs4246977T were abundant in HLs and hence associated with HA adaptation (Table 1 and Supplementary Tables S4A, S5A, S6A, S7A and S11A). A comparison of the HLs with Han Chinese (CHB+CHS) and Japanese (JPT) populations provided mixed results, with few of the genes/SNPs having near similar distribution in the three populations (Subtables B in Supplementary Tables S4-S12).

Haplotypes Distribution Is Specific in HAPE and HA Adaptation
Seventy-two haplotypes with frequency > 2% were observed in eight telomere-maintaining genes across five chromosomes. Three haplotypes, namely, GAAA of TERF1, TCCC-GACTTA and TTCCGATCTA of TEP1+HSP90AA1, were observed in HAPE protection ( Table 2).  In highlanders, nine adaptive haplotypes were observed (Table 2). TERT was observed in the adaptive haplotype GATCGCCAG. The AAGG haplotype of TERF1 included adaptive alleles TERF1 rs2975843A and TERF1 rs10099824G located at the second and fourth contig position, respectively. Five haplotypes of ACD+TERF2+TERF2IP were adaptive, i.e., AAAAATTCGCA, AAGTATCCGCA, AAGTATTCGCA, AAGTGTCCGCA and TG-GTATCCGCA; these were comprised of adaptive alleles ACD rs6979A or TERF2 rs3785073G or both at second and third position, respectively. Further, TCCCGGCTCA and TTCCG-GCTCA in TEP1+HSP90AA1 were adaptive, having adaptive alleles TEP1 rs2184282C and TEP1 rs4246977T at the seventh and eighth positions, respectively. The haplotypes could not be generated in CHB+CHS and JPT populations as the genotype data of only 49% of SNPs were available.  and HLs (high-land natives). (A) HAPE-associating shelterin genes in a 2-locus best model. The TT-TT (TPP1/ACD rs72556537A/T-RAP1/TERF2IP rs59297469G/T) interaction was observed as risk in which the protective allele TPP1 rs72556537A was absent. (B) HAPE-associating telomerase genes in a 2-locus best model. A risk-predicting TT-AG (TEP1 rs2228036G/T-TERC rs2293607A/G) interaction was observed in which risk allele TEP1 rs2228036T was present as homozygous genotype TEP1 rs2228036TT and risk allele TERC rs2293607G was present as heterozygous rs2293607AG genotype. (C) HAPE-associating genes of the telomere-telomerase system in a 4-locus best model. The TT-GT-AA-TT (RAP1 rs59297469G/T-TEP1 rs2228036G/T-TERC rs2293607A/G-TPP1 rs72556537A/T) risk-predicting interaction consisted of risk allele RAP1 rs59297469T as homozygous rs59297469TT genotype and risk allele TEP1 rs2228036T as heterozygous rs2228036GT genotype. The GG-GT-AG-TT interaction consisted of risk allele TEP1 rs2228036T as heterozygous rs2228036GT genotype and risk allele TERC rs2293607G as heterozygous rs2293607AG genotype. (D) HA adaptation-associating shelterin genes and HA adaptation-associating genes of the telomere-telomerase system in a 3-locus best model. Two interactions, GG-AA-TC (TERF1 rs10099824G/A-TPP1 rs6979G/A-TERF1 rs2975843T/C) and GG-AA-TT, predicted HA adaptation consisting of adaptive allele TERF1 rs10099824G as homozygous rs10099824GG genotype and adaptive allele TPP1 rs6979A as homozygous rs6979AA genotype. The adaptive allele TERF1 rs2975843T was present as heterozygous rs2975843TC genotype in GG-AA-TC and as homozygous rs2975843TT genotype in GG-AA-TT. (E) HA adaptation-associating telomerase genes in a 2-locus best model. Three adaptive interactions, i.e., AA-TT, AG-TT. GG-TT (TEP1 rs2184282A/G-TEP1 rs4246977T/C) were obtained. The adaptive allele TEP1 rs4246977T as homozygous rs4246977TT genotype is present in all the three adaptive interactions revealing a strong association with HA adaptation. In (A-C) presentation, the dark grey cells associate with HAPE risk, and the light grey cells associate with HAPE protection; the first bar in a cell represents the number of HAPE-p and the second bar represents the number of HAPE-f. In (D,E), the dark grey cells represent interactions enriched in HAPE-f and the light grey cells associate with HA adaptation; the first bar in a cell represents the number of HAPE-f and the second bar represents the number of HLs. * marked cells represent the significant genotype interactions, which are summarized in respective tables. The genotypes of the SNPs appear in alphabetical order. Chi-square test or Fisher's exact test (when the number of samples was ≤5) was applied to calculate p-value, OR (odds ratio) and 95% CI (confidence interval) using SISA online tool.

Best Interaction Models for HA Adaptation
The HA adaptation predicting MDR model among shelterin genes was represented by a 3-locus model comprising TERF1 rs10099824G/A, TERF1 rs2975843T/C and ACD rs6979G/A (TA = 0.77, CVC = 10/10, p < 0.0001) ( Figure 3D). The telomerase genes were represented by a 2-locus best model comprising SNPs TEP1 rs2184282A/G and rs4246977T/C (TA = 0.65, CVC = 10/10, p < 0.0001) ( Figure 3E). Interestingly, when all the HA adaptationassociating SNPs were combined, the same 3-locus model of shelterin genes emerged as the best model with the same genotype combinations suggesting a strong interaction among the three SNPs (TA = 0.77, CVC = 10/10, p < 0.0001) ( Figure 3D); thus, a preference for shelterin genes was visible. MDR analysis among a total of 43 SNPs did not reveal any significant interacting model for HAPE susceptibility and HA adaptation.

Expression Profile Differed for Telomere Maintaining Genes
The expression profile of the 11 telomere-maintaining genes was examined in 18 HAPEf, 10 HAPE-p and 17 HLs subjects. A significant differential expression was observed for genes TERF2IP, TERT and HSP90AA1. TERF2IP was down-regulated by 3.1 fold in HAPE-p and up-regulated by 1.4 fold in HLs ( Figure 4A). The contrasting expression of HSP90AA1 in HAPE-p and HLs identifies its importance in HA physiology. TERT was down-regulated by 1.8 fold in HAPE-p and 1.3 fold in HLs ( Figure 4A). HSP90AA1 was down-regulated 2.1 fold in HAPE-p but up-regulated by 1.1 fold in HLs ( Figure 4A), which was in line with longer telomere length, higher telomerase activity and high levels of NO in HLs compared to HAPE-p. Further HSP90AA1 was down-regulated 2.1 folds in HAPE-p. * represents p ≤ 0.05, ** represents p ≤ 0.01 and **** represents p ≤ 0.0001. (B) Docking image of TFs, USF1 and N-myc on RAP1 or TERF2IP with wild type allele, rs59297469G (WtRAP1). (C) Docking image of TFs, USF1 and N-myc, on RAP1 or TERF2IP with variant/risk type allele, rs59297469T (VtRAP1). Figure (B,C) differentiates the binding affinity for the same locus in the presence of risk and the protective allele that would influence the expression and function of the gene, which is depicted in the subsequent presentation. RMS, root mean square. (D) RAP1 expression is significantly down-regulated with a higher ∆Ct in the presence of risk allele RAP1 rs59297469T compared to the protective allele rs59297469G. (E) Interactions among the telomere-maintaining proteins. (F) Telomere system shows stronger interactions among its proteins contributing to the regulation of other pathways such as the nitric oxide synthase (NOS) (marked red and blue) and cellular response to stress (marked green and yellow) which are depicted in subsequent presentation. (G) Enriched biological processes such as the nitric oxide synthase (NOS) and cellular response, which are important pathways in HAPE pathophysiology. (H) Interactions among the proteins of the telomere-maintaining and the hypoxia-maintaining systems in HA physiology influence several pathways. (I) Prominent biological processes as a result of the interaction of the proteins of the two major systems. In silico tools were applied to gain insight into the differential physiological system.

In Silico Analysis Underlined the Association of Decreased TERF2IP Expression with Risk Alleles in HAPE
Among the genes with differential gene expression, the TERF2IP T allele of the promoter SNP rs59297469G/T emerged as prominent with HAPE risk. Hence, in silico tools were applied to determine the influence of this allele on TERF2IP (also designated as RAP1) expression, especially through TFs. Identifying TFs, and their binding to a gene are key steps to understanding the transcriptional regulation in the diseased condition. Here, two TFs, Upstream stimulatory factor 1 (USF1) and N-myc proto-oncogene protein (N-myc), emerged as significant. The TFs were further checked for the change in binding affinity in the presence of the wild type (Wt, rs59297469G) and variant/risk type (Vt, rs59297469T) allele on the TERF2IP promoter. TFs were docked on the specific promoter sites, and the Hex program identified the binding affinities, whereas PyMOL was used to visualize the docked structure. The 3D structure of USF1 (PDB ID: 1AN4) and N-myc (PDB ID: 5G1X) had a resolution of 2.90 Å and 1.72 Å, respectively. The Hex docking tool identified 650 clusters from 1000 solutions in 5.25 s with a radius cut-off of 20.0. A 250 × 250 × 250 grid was generated around the TERF2IP-TF complex with a prefix root mean square of −1.00. HADDOCK and CLusPRO supported the docking by Hex.
The binding affinity of WtRAP1-USF1 was −693.9 kcal/mol ( Figure 4B), and of VtRAP1-USF1 was −681.1 kcal/mol ( Figure 4C), while the binding affinity of WtRAP1-N-myc was −675.5 kcal/mol ( Figure 4B) and of VtRAP1-N-myc was −605.5 kcal/mol ( Figure 4C). Elevations in binding affinity scores of both the TFs when interacting with VtRAP1 suggested weaker intermolecular forces between TFs and TERF2IP as compared to the WtRAP1-TF complex, leading to shorter residence time at the binding site. Thus, in silico analysis suggested apparent transcription of TERF2IP in the presence of Wt rs59297469G, while down-regulation of TERF2IP in the presence of risk allele rs59297469T.
To further ascertain the association of risk allele rs59297469T with the down-regulation of TERF2IP in HAPE-p, we correlated the ∆Ct values of each patient in the HAPE-p group with the rs59297469 genotypes. It was observed that the presence of risk allele rs59297469T in the GT genotype increased the ∆Ct ( Figure 4D) which meant TERF2IP expression decreased, which is in agreement with the findings of the in silico analysis. The TT genotype was not observed in the HAPE expression samples.

Protein-Protein Interactions
The interactions among the respective 10 shelterin proteins highlighted close bonding ( Figure 4E). Biological functions related to telomere extension and its maintenance were among the top priority highlighting the relevance in the telomere-telomerase system. Regulation of NOS3 activity and stress response were inclined towards TERT, HSP90AA1 and TERF2; additionally, the stress response was also inclined towards ACD, TERF2IP and PTGES3 ( Figure 4F,G).
To further identify interactions and protein networks associated with the telomeretelomerase system, we included a few of the most specific proteins of the vital HA physiology systems, namely, HIF-1-alpha (HIF1A), Egl-9 family hypoxia-inducible factor 1 (EGLN1), Endothelial PAS Domain Protein 1 (EPAS1), NOS3 and endothelin-1 (EDN1). These proteins bonded well with TERT, HSP90AA1 and PTGES3 of the telomere-telomerase system ( Figure 4H). The enriched biological processes included telomere maintenance, telomerase activity, homeostasis, oxidoreductase activity, and stress response, which are essential for HA acclimatization and adaptation ( Figure 4I).

Correlation between Biomarkers and Clinical Parameters
Correlation analyses among various clinical and biochemical parameters were mixed. The correlations were not encouraging between biomarkers and clinical parameters in HAPE-f and HAPE-p but they were in HLs. In HLs, total antioxidant activity correlated negatively with MAP (r = −0.25, p = 0.003) and positively with SaO 2 (r = 0.33, p = 5.35 × 10 −5 ), suggesting that with increased antioxidant activity, highlanders would have low MAP and heightened oxygenation, vital to the hypobaric hypoxic environment (Supplementary Table S14). Correlation analyses among the various biomarkers revealed an inverse correlation between total antioxidant activity and 8-isoPGF2α in HLs (r = −0.23, p = 0.005), indicating a healthy state under the given environment (Supplementary Table S15).

Discussion
The interplay between environmental, biochemical and genetic factors manipulates several physiological pathways; among these, the telomere-telomerase system is one of the crucial and sensitive systems. We expectedly observed its aberration under the hypobaric hypoxic environment of HA with shorter telomere length and lower telomerase activity influencing the health and disease cycle, HAPE susceptibility and HA adaptation. The telomere attrition may be because of increased ROS and UV radiation at HA. ROS are known to cause multiple types of DNA damage, including oxidized bases, and singlestrand and double-strand breaks throughout the genome [40,41]. The consequence of such a development is chromosomal instability and decreased cellular viability, thereby increasing cellular dysfunction and diverting the body's physiology toward disorder. The overall effect could be the dysfunction of proteins of the telomere-telomerase system contributing to disturbed body physiology. However, literature is scant on the permanent residents of HA or the sojourners. Nonetheless, it is known that shorter telomere length is associated with COPD risk, mortality in COPD and reduced pulmonary function [20,42]. Thus, shorter telomere length in HAPE patients may associate with poor pulmonary function. We could relate this phenomenon to the significantly differential SaO 2 levels and ROS in HAPE and health. SaO 2 drops on exposure to HA or a hypobaric environment but rises to near normal with acclimatization [12,14]. In individuals who fail to acclimatize, such as HAPE patients, the SaO 2 levels are decreased significantly and could even lead to mortality if not treated immediately. SaO 2 in both healthy sojourners and highlanders is generally near normal and not as low as in the HAPE patients [43][44][45]. Among the various factors influencing the body, oxygen saturation, and oxidative stress seems the most effective, other than the hypoxia-sensing pathway [13,43].
The relevance of oxidative stress biomarkers, i.e., 8-iso PGF2α and total antioxidant activity, was apparent in HAPE-p with elevated 8-isoPGF2α levels suggesting increased oxidants in these individuals, who had decreased enzyme activity and as an upshot, shorter telomeres. ROS are initially produced to stimulate a series of mechanisms allowing acclimatization and adaptation to environmental changes. However, the accumulation of ROS leads to oxidative stress and causes the oxidation of several biomolecules in a vascular system [46,47]. ROS scavenges NO and reduces its production, accompanied by a reduction in other vasodilators, thereby impairing normal vascular physiology [13,48]. We and others had previously reported decreased NO levels in HAPE patients and the association of NOS3 variants with the HAPE risk [1,7,13,[49][50][51][52]. The observation also suggested that similar antioxidant activity in sojourn groups is optimal for maintaining ROS in healthy individuals but unable to normalize exaggerated ROS in HAPE, thus sustaining oxidative stress.
Genetic analyses of the telomere-telomerase system made equally credible revelations, with genes clearly associating with adaptation and maladaptation (disease). The promoter SNP rs59297469G/T of TERF2IP, the exonic SNPs TERC rs2293607A/G and TEP1 rs2228036G/T strongly pointed to the shelterin and telomerase dysfunction and thereby the risk of HAPE. The subdued TERF2IP expression with the risk allele suggested its regulatory influence. We subsequently predicted allele-specific binding of TFs in the case of TERF2IP SNP rs59297469G/T which was corroborated by in silico analyses. The analysis identified two TFs, USF1 and N-myc, for the SNP, wherein each allele had a varied binding affinity for the TFs and thus varied transcription regulation, i.e., apparent transcription of TERF2IP in the presence of rs59297469G, while down-regulation in the presence of risk allele rs59297469T. This was further ascertained with the association analysis of the ∆Ct values of rs59297469 genotypes for each patient that revealed elevated ∆Ct in the presence of risk allele rs59297469T in the GT genotype of each patient, consequent to this decreased gene expression. Incidentally, the TT genotype was not observed in the evaluated samples to ascertain the double allelic effect. Needless to say, TERF2IP regulates telomere length and the expression of other genes [53,54]. Among the other relevant polymorphisms of genes, TERC SNP rs2293607A/G is present in the H/ACA domain of the TERC RNA, which is required for the accumulation and stability of the telomerase [55]. Whereas the ACD/TPP1 3'UTR SNP rs72556537T/A is associated with HAPE protection. TPP1 appears a crucial protein for acclimatization because it assists in the recruitment of telomerase onto telomeres and promotes telomerase processivity [29,31,56]. Likewise, the multiallelic involvement in protection was apparent with three protective haplotypes involving TERF1, TEP1 and HSP90AA1, and also by gene-gene interactions.
Although the genetic association of TERT and HSP90AA1 was not very prominent in HAPE susceptibility, both genes were significantly down-regulated in HAPE-p along with TERF2IP. Down-regulation of telomerase subunits TERT and HSP90AA1 supported the decreased telomerase activity observed in HAPE-p. TERT prevents mitochondrial DNA damage, increases mitochondrial efficiency, and regulates ROS production under hypoxia; thus, its down-regulation not only suggests telomere shortening and telomerase insufficiency but also increased oxidative stress [34,35]. Down-regulated TERT may weaken protection against oxidative stress by increasing mitochondrial ROS production, reducing the protection of mitochondrial DNA from oxidative damage, and undermining mitochondrial function [34,35]. TERT also has transcriptional efficiency in influencing downstream genes that together may be crucial for HA physiology [57]. Whereas HSP90 is essential for the maintenance of telomere length, stabilization of telomerase structure, telomerase activity and regulation of NOS3 to produce NO or ROS [36,37,58,59]. Thus, HSP90AA1 can influence HA physiology by regulating the telomere-telomerase system, vascular homeostasis and oxidative stress pathways.
The protein-protein networking also revealed a complex interplay via association among various pathways. We observed deviations in highlanders compared to HAPE patients with longer telomere length, higher telomerase activity, 8-isoPGF2α levels, total antioxidant activity and SaO 2 and normal MAP, suggesting an adaptive physiology [43,44]. The expression profile also suggested an altered physiology for HLs compared to HAPE. Up-regulation of TERF2IP and HSP90AA1 and down-regulation of TERT expression were observed in HLs. Another important shelterin, POT1, represented by the intronic SNP rs7794637C/T, is associated with HA adaptation. It is a genetic marker for longevity, regulates telomere elongation by inhibiting telomerase, maintains telomere structure and inhibits the DNA damage response [16,60,61]. We further observed two SNPs, the intronic rs2184282A/G and the promoter rs4246977T/C from the telomerase subunit TEP1 associating with adaptation. Interestingly, TEP1 can modulate telomere length, not via telomerase, but by alternative lengthening through associated telomere pathways [32]. Thus, the individual gene and each respective genotype emerged relevant to the telomere-telomerase system. In addition, the haplotypes and the multiallelic within genes and between genes interactions were amply visible in diseases and adaptation. Among the genes, TERF1, TERF2, POT1, ACD, TEP1, TERF2IP, TERT and HSP90AA1 were at the forefront, whereas the genegene interactions showed a stronger association between TERF1 and TPP1. Nevertheless, interactions of genes of the telomere-telomerase system with genes of other pathways such as vascular homeostasis may also exist in HA adaptation. A schematic representation of interactions of various molecules, including the telomere-telomerase system in HAPE pathophysiology, is proposed in Figure 5.
With the proposed role of shorter telomere length in HAPE pathophysiology ( Figure 5), shorter telomere length among endurance athletes is observed as a stressor in their skeletal muscles [62]. On the contrary, reports have associated higher physical activity levels and exercise training with longer telomere length [23][24][25][26][27]. Our present findings would correlate with athletes that utilize altitude training to increase their exercise capacity; where being in a hypoxic environment would improve their performance at sea level or acclimatize to sports events at HA. A recent study reported that repeated and continuous exposure to HA might significantly alter oxidative stress and antioxidant capacity among HA visitors [63]. Another review highlighting the effect on telomere length and telomerase activity due to physical exercise concluded that regular aerobic activity, of moderate to vigorous intensity, preserves telomere length among athletes [27]. Thus, longer telomere length and higher telomerase activity in response to more physical activity at HA have the potential to improve/provide therapy in better adaptation at HA for mountain sports and high-altitude mountaineering. Figure 5. The proposed pathophysiology of HAPE. Hypobaric hypoxia is a stimulant to stress, which may in turn stimulate a normal physiological system to sustain a normal status which may not be possible in a few subjects. The same is depicted in the given graphical sketch. The expression of various proteins and genes, which are essential to regulate ROS production is dysregulated on exposure to HA. The telomere system is one of the major systems, whose genes/proteins are significantly disturbed. Down-regulation of HSP90AA1, TERT and the telomerase is observed and contributes to telomere attrition and may also enhance ROS production in a reversible process. Overall, in this milieu, the exaggerated ROS and the down-regulated candidate markers lead to the NOS3-mediated decrease in NO production and enhanced vasoconstriction leading to the disturbance of vascular homeostasis causing pulmonary hypertension, a hallmark of HAPE.
The present study is not without limitations. First, despite the higher sample size of each cohort, replicating our findings in a much larger cohort will confirm our findings. Second, we could not investigate the HAPE patients post-recovery to obtain a veritable picture of the studied parameters and whether they were truly involved in the pathophysiology. It is, however, a difficult task to follow the patients, who are tourists, under mental trauma of becoming severely ill on vacations. They hence immediately leave the hospital upon recovery and refuse to participate in the follow-up study. A longitudinal study though would be beneficial. Third, we could not check the role of longer telomeres and telomerase activity in high-enduring sports persons; thus, replicating our findings in high-performance sports persons will confirm the role of the telomere system in athletes. Moreover, investigating multigenetic systems is desired, especially the haplotypes involving the short and long telomeres. Although we analyzed the variants (short-read telomeres) that provide an accuracy of greater than 99% (similar to long-read telomers) over a few megabases of the genome [64][65][66], investigating the long-range telomere haplotypes along with short-read telomeres in a disease state or high-endurance sports is desired. Lastly, additional experimental studies determining the cellular effect of significant genes with a human cell line would authenticate their functional consequence.

Conclusions
The present study reports the involvement of the telomere-telomerase system in HA physiology for the first time. The shorter telomere length and lower telomerase activity define HAPE; the longer telomere length and higher telomerase activity define health and adaptation at HA. Various regulatory, exonic or intronic SNPs and haplotypes of this system are crucial for HAPE susceptibility and HA adaptation. A complex interaction exists between these genes and genes of other pathways such as vascular homeostasis. A study involving more SNPs and genes from the telomere-telomerase system and different established pathways at HA might be helpful to unfold the proper framework of HAPE pathophysiology and HA adaptation. The telomere-telomerase system holds promise for future studies.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijerph20031935/s1, Figure S1: Position of each SNP and allele conversion followed in the haplotypes; Figure S2: Linkage disequilibrium plots of (A) TERF1, (B) TERF2, (C) POT1, (D) ACD, (E) TERF2IP, (F) TERT, and (G) TEP1; Table S1: Clinical characteristics of the three study groups; Table S2: Position of the selected SNPs on the genes; Table S3: qRT-PCR primers for gene expression profiling; Table S4 Table S13: Genotype and allele distribution and comparison of PTGES3 SNPs in HAPE-f, HAPE-p and HLs; Table S14: Correlation analysis between biomarkers and clinical parameters in the three study groups; Table S15: Correlation analysis between the biomarkers in the three study groups.
Author Contributions: Q.P.; conceived, designed, directed the study, analyzed, and interpreted the data, and wrote and revised the manuscript. M.R.; performed the estimation of telomere length and telomerase activity, gene expression, genotyping, and statistical analyses. M.R. and S.T.; estimated oxidative stress markers. H.K.; performed in silico analyses. G.M., T.T. and Q.P.; recruited patients and gathered detailed clinical information for the study. M.R. and Q.P.; wrote the manuscript with critical inputs from all authors. All authors have read and agreed to the published version of the manuscript.