Changes in Nuclear Gene Expression Related to Mitochondrial Function Affect Extracellular Matrix, Collagens, and Focal Adhesion in Keratoconus

Abstract

Purpose: Mitochondrial DNA (mtDNA) abnormalities were previously found to be causative in the pathogenesis of various diseases. Here, comprehensive mitochondrial and nuclear sequence and transcript analyses, along with analyses of the methylation aspects of nuclear genes related to mitochondrial function, were performed in patients with keratoconus (KTCN) to evaluate their contribution to the KTCN pathogenesis.

Methods: Blood mtDNA of 42 KTCN and 51 non-KTCN individuals was Sanger sequenced and analyzed along with the previously obtained corneal RNA-sequencing data of 20 KTCN and 21 non-KTCN individuals. In addition, the expression and methylation of mtDNA genes and 1223 mitochondria-related nuclear genes were evaluated.

Results: The mtDNA sequence alterations detected in blood coincided with variants identified in transcripts of the matched corneal tissues. In KTCN corneas, 97 mitochondria-related genes were deregulated, including TGFB1, P4HB, and BCL2, which are involved in the extracellular matrix (ECM) organization, collagen formation, and focal adhesion pathways. No changes in the expression of mtDNA transcripts and no differentially methylated genes among the assessed mitochondrial–nuclear gene sets were found.

Conclusions: The absence of corneal-specific mtDNA variants indicates that there is no direct relationship between mitochondrial sequence variability and KTCN phenotype in the studied individuals. However, the identified KTCN-specific transcriptomic alterations of the nuclear genes directly related to the mitochondria functioning point to their possible involvement in the ECM organization, collagen formation, and focal adhesion.

Translational Relevance: The identification of abnormalities within nuclear genes regulating ECM formation, collagen synthesis, and/or focal adhesion may form the basis of future treatment strategies or predict the progression of corneal changes in KTCN.

Introduction

Due to its unique physiology and location, the cornea receives a significant amount of high-tension atmospheric oxygen and sunlight, including visible light and ultraviolet radiation.¹ The exposure to environmental stress induces generation of reactive oxygen species (ROS), leading to oxidative stress.¹ However, in response to oxidative stress, cells can produce antioxidant enzymes, including superoxide dismutase, catalase, or glutathione peroxidase, protecting the tissue from oxidative damage.² 

Increased oxidative stress has been reported in keratoconus (KTCN),^3–7 a corneal disease characterized by progressive thinning and anterior protrusion of the cornea leading to reduced vision and changed refractive power.⁸ The underlying mechanisms of oxidative changes in KTCN corneas are not fully understood.⁹ Nevertheless, because mitochondria play an important role in maintaining oxidative stress, it has been shown that the oxidative damages observed in KTCN could be related to mitochondrial dysfunction. 

In support of this hypothesis, Pathak et al.¹⁰ identified 84 variants in mitochondrial complex I genes (ND1–ND6) in 20 Indian patients with KTCN, and Abu-Amero et al.¹¹ revealed potential KTCN-related mitochondrial DNA (mtDNA) variants, using sequencing of the whole mitochondrial genome in a group of 26 Saudi patients with KTCN.¹¹ 

It has also been observed that relative mtDNA content and integrity were lower in KTCN corneas due to decreased levels of nuclear TFAM transcript.⁷ In contrast, higher expression of other nuclear genes, POLRMT and TFB2M, resulted in significantly increased transcript levels of three mtDNA genes: ND1, ND6, and COX1.⁷ Together, these findings indicate the involvement of increased mtDNA damage, increased expression in progressive respiratory chain dysfunction, and elevated ROS levels in the KTCN corneas,⁷ suggesting cooperation between mtDNA and nuclear-encoded genes. 

In this study, we performed a comprehensive analysis of mtDNA and mitochondrial corneal transcripts in individuals with KTCN. We aimed to determine whether KTCN patients have more significant mtDNA damage than individuals without KTCN and to evaluate the co-occurrence of mtDNA alterations and transcripts variants. We have also used our previously obtained RNA and DNA data derived from KTCN and non-KTCN corneas to assess the contribution of abnormal expression or methylation of nuclear genes that may be related to mitochondrial function in KTCN pathogenesis in Polish patients. 

Materials and Methods

Patients’ Involvement

Blood DNA samples obtained from 42 Polish patients with sporadic KTCN (29 males and 13 females; ages 18–69 years) and 51 individuals without KTCN phenotype (non-KTCN; 25 males and 26 females; ages 21–80 years) were assessed in this study. All individuals were ascertained and examined as previously described.¹² Briefly, the diagnosis of KTCN was based on best-corrected visual acuity testing, intraocular pressure assessment, a slit-lamp examination (Orbscan IIz; Bausch & Lomb Surgical, Rochester, NY, USA), and videokeratography data (Supplementary Table S1). Following an explanation of the nature and possible consequences of the study, written consent forms were obtained from all participants in accordance with the tenets of the Declaration of Helsinki. The research protocol was approved by the Institutional Review Board at Poznan University of Medical Sciences in Poland. 

Expression and Methylation Analyses of Nuclear Genes Related to Mitochondria Functioning

The expression level values of 1223 nuclear genes, including 1075 genes with strong evidence for mitochondrial location (listed in the MitoCarta 2.0),¹³ 103 genes from the oxidative phosphorylation pathway (OXPHOS; Kyoto Encyclopedia of Genes and Genomes: hsa00190),¹⁴ and 45 genes involved in mitochondrial function and/or disease based on previous literature data,^15–19 were obtained from our RNA-sequencing (RNA-seq) study.^20,21 

Methylation levels of the same 1223 nuclear genes related to mitochondria functioning were obtained from the reduced representation bisulfite sequencing (RRBS) experiment.²² 

The expression analysis also encompassed genes encoded by the mitochondrial genome. Differential expression analysis was performed using a previously published protocol²³ using the limma package.²⁴ Genes were determined to be differentially expressed based on the 0.01 false discovery rate threshold value, and 1.5-fold change cutoff. 

Whole Mitochondrial Genome Sequencing

DNA was extracted from the peripheral blood samples of all KTCN and non-KTCN individuals with the use of the Gentra Puregene Blood Kit (QIAGEN, Hilden, Germany) as previously described.¹² Mitochondrial DNA was amplified using Taq DNA Polymerase (Thermo Fisher Scientific, Waltham, MA) in 24 overlapping amplicons covering the whole mitochondrial genome. PCR products were cleaned using FastAP Thermosensitive Alkaline Phosphatase and Exonuclease I (Thermo Fisher Scientific) and Sanger sequenced using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA) according to the manufacturer's protocols. Following precipitation, samples were analyzed on an ABI Prism 3730xl genetic analyzer (Applied Biosystems), visualized using Sequencher 5.0 Software (GeneCodes Corporation, Ann Arbor, MI) and aligned with the reference sequence of mtDNA (GRCh37/hg19, NC_012920.1). 

In addition to Sanger sequencing, one mtDNA sample (KC42) was subjected to next-generation sequencing (NGS). Briefly, the mtDNA sample was enriched by long-range PCR in two overlapping amplicons, and DNA libraries were prepared using the Nextera XT Sample Preparation kit (Illumina, San Diego, CA) as described in the Illumina Human mtDNA Genome Guide. Paired-end indexed sequencing (2 × 150 bp) was performed on the MiSeq system (Illumina). Sequence readouts were initially analyzed by CASAVA software (Illumina) to generate reads in fastq format. After the quality control step, including adapter trimming and low-quality reads removal, reads were aligned to the mtDNA reference sequence (revised Cambridge Reference Sequence [rCRS], GenBank accession number NC_012920.1) with the Burrows–Wheeler Alignment Tool (http://bio-bwa.sourceforge.net/). They were further processed by the Picard and Genome Analysis Toolkit (http://broadinstitute.github.io/picard/) and then annotated using frequency and phenotype data from the MITOMAP database. 

Prediction of the Effect of Amino Acid Substitution on Protein Structure and Function

The consequences of amino acid substitutions in the mtDNA were analyzed by five different prediction algorithms: PolyPhen-2,²⁵ PROVEAN,^26,27 SNPs&GO,^28,29 MitoClass 1,³⁰ and APOGEE.³¹ 

Inter-Tissue mtDNA Sequence Variants Assessment

To determine possible inter-tissue differences of mtDNA sequence variants, raw mitochondrial transcriptome (RNA-seq) data available for the corneas of 20 Polish patients with KTCN and 21 non-KTCN Polish individuals^20,21 were re-analyzed and compared with the mtDNA sequence data of the blood samples obtained in this study. Mapped reads spanning introns from the RNA-seq were split separately using the module “SplitNCigarReads” in the GATK toolkit.³² Variants from Ensembl exons were called by FreeBayes.³³ The obtained variants were filtered according to the Utah Genome Project Variant Calling Protocol and then compared to mtDNA sequencing results. 

Haplogroup Analysis

Mitochondrial DNA sequencing data were used to determine the haplogroups of the 93 analyzed samples. The mtDNA haplotypes were affiliated with haplogroups by HaploGrep2³⁴ following Phylotree Build 17.^35,36 This tool first allocates mtDNA profiles to major haplogroups and then determines whether all haplogroup-specific mutations are observed in the profile of the analyzed sample. Finally, it connects the evolutionary pathways with the targeted sequence at the rCRS. 

Statistical Analyses of mtDNA Variant, Allele, and Haplogroup Distributions

Allele distributions of mtDNA variants and haplogroup distributions among patients with KTCN and non-KTCN individuals were evaluated using Pearson's χ² test. P < 0.05 was considered statistically significant. 

Results

Nuclear Genes Differential Expression

Of the 1223 nuclear-encoded genes related to mitochondrial function,^13–19 56 and 41 were significantly down- and upregulated, respectively, in KTCN corneas compared with non-KTCN corneas based on our previous RNA-seq study^20,21,37 (Supplementary Table S2). Seven of these 97 differentially expressed genes overlapped those from the top overrepresented molecular pathways that were previously described as deregulated in KTCN, including TGFB1, P4HB, and BCL2, which are involved in the extracellular matrix (ECM) organization, collagen formation, focal adhesion, non-integrin membrane–ECM interactions, and elastic fiber formation (Table 1). 

Table 1.

 

Differentially Expressed Nuclear-Encoded Genes Related to Mitochondrial Function Involved in Deregulated Molecular Pathways That Were Overrepresented in KTCN Samples, Based on the RNA-Seq Study^20,21

 

No changes in the expression of mitochondrial transcripts were observed in KTCN individuals. 

Methylation Profile of Mitochondria-Related Genes

The differentially methylated genes detected previously in KTCN corneas²² have been compared with the list of nuclear-encoded genes related to mitochondria function. No differentially methylated nuclear genes were detected among 1233 genes^13–19 from the mitochondria-related gene set. 

mtDNA Sequence Variants Identification

In total, 358 sequence variants in 93 mtDNA samples were identified. Among them, 143 were present in KTCN patients only, and 138 were exclusive for non-KTCN individuals. Seventy-seven single nucleotide polymorphisms (SNPs) were common for both analyzed groups: KTCN (out of 220 sequence variants) and non-KTCN (out of 215 sequence variants). Supplementary Table S3 presents all identified variants, localized within both the mtDNA genes and intergenic regions. 

The highest number of sequence variants (n = 54) was observed in the MT-ND5 gene, in which a similar distribution in KTCN and non-KTCN individuals was found. The lowest number of variants was observed in tRNAs genes, and no variants were found in MT-TI, MT-TQ, and MT-TM genes. The majority of sequence variants (n = 258) were identified in a single KTCN/non-KTCN sample only. The frequency of variants identified in this study did not differ significantly between KTCN and non-KTCN individuals (P > 0.05) (Supplementary Table S3). The average number of variants per individual was 15 and 13 for KTCN and non-KTCN individuals, respectively. The number of alterations per sample ranged from one in the sample KR5 (non-KTCN) to 57 in the sample KC8 (KTCN; see Supplementary Table S4). The number of sequence variants observed in individual mtDNA samples did not correlate with the patient's age (Fig.). We have not identified any mtDNA sequence variants that were significantly more represented in KTCN patients compared with individuals without KTCN. 

Figure.

 

Relationship between age at the clinical examination and numbers of sequence variants observed in mtDNA within the studied patients with KTCN (A) and non-KTCN individuals (B).

Analysis of the KC42 sample, assessed by NGS, revealed 39 sequence variants in the mitochondrial genome, including two indels and 37 nucleotide substitutions. All but one variant (m.8774A>G) were confirmed using Sanger sequencing analysis. In addition, NGS sequencing data indicated five heteroplasmic variants: m.8774A>G (variant allelic fraction [VAF] = 0.5), m.16179delCAA (VAF = 0.5), m.16182A>C (VAF = 0.5), m.16183A>C (VAF = 0.5), and m.16189T>C (VAF = 0.5). 

Amino Acid Substitutions and Their Impact on Protein Function

Out of 358 identified sequence variants, 97 resulted in amino acid missense substitutions (Supplementary Table S5). Analysis of the potential impact of these changes on protein function revealed that five, eight, and 15 SNPs were determined to be deleterious by four, three, and two algorithms, respectively (Table 2). In total, 52 variants were predicted to be damaging by at least one tool. The highest impact on protein structure and function was predicted for rs2298008 (m.8921G>A), localized within the MT-ATP6 gene and detected in one KTCN sample (KC24). Most of the variants, both possibly deleterious and benign, were similarly distributed in KTCN and non-KTCN samples (Table 2). 

Table 2.

 

Sequence Variants Causing Possibly Damaging Amino Acid Substitutions

 

Comparison of mtDNA Sequence Variants in Blood and Corneal Tissue Samples

A comparison of sequence variants found in blood mtDNA and corneal mitochondrial transcripts in 20 patients with KTCN and 21 non-KTCN individuals has revealed 174 variants (Supplementary Table S6). One hundred and forty-nine sequence variants were detected simultaneously in both mtDNA sequencing and mitochondrial transcriptome (RNA-seq) data. Eight variants were present in the RNA-seq data but not in the mtDNA sequencing data. On the other hand, 15 variants were observed in mtDNA sequencing results but not in the transcriptome data. 

Haplogroup Analysis

Fifteen haplogroups were revealed based on mtDNA sequencing results (Supplementary Table S7). Haplogroup H was the most common, as it was identified in 17 KTCN patients and 17 non-KTCN individuals. The second most frequent was haplogroup U, which was present more frequently in the non-KTCN individuals (n = 12) compared with the KTCN individuals (n = 5). The remaining 14 haplogroups were similarly represented in KTCN and non-KTCN samples. No significant differences in haplogroup distributions were observed. 

Discussion

Previous studies indicated that KTCN corneal fibroblast cells had increased ROS levels and lower mitochondrial membrane potentials compared with corneas derived from unaffected individuals,^5,38,39 suggesting that oxidative stress and mitochondrial abnormalities might be involved in KTCN. Although mitochondrial function depends on mitochondrial genes, it also requires a crosstalk between mtDNA and nuclear genes to maintain tissue homeostasis and proper oxidative metabolism.⁴⁰ Thus, to check if there were any abnormalities in the expression of mitochondria-related genes in KTCN patients compared with non-KTCN individuals, we analyzed a set of 1233 nuclear genes^13–19 that could be involved in mitochondrial metabolism. Based on our previous RNA-seq study,^20,21 we found that 97 nuclear encoded genes involved in mitochondrial function were down- or upregulated in KTCN corneas. Interestingly, three genes with altered expression levels, P4HB, TGFB1, and BCL2, belong to pathways regulating collagen synthesis, focal adhesion, or ECM function. These pathways influence corneal organization, and their disruption was indicated as a potential factor involved in KTCN pathogenesis.²⁰ 

BCL2 protein is localized in the outer mitochondrial membrane and plays an essential role in inhibiting the pro-apoptotic proteins.^41,42 It is also involved in the regulation of corneal homeostasis.⁴³ Transforming growth factor beta 1 (TGF-β1), through binding to its receptors, regulates the assembly of ECM and acts as a stimulator of a profibrotic response in the injured cornea.⁴⁴ Analysis of murine fibroblasts NIH/3T3 revealed that TGF-β1 induces differentiation of fibroblasts, which requires increased mitochondrial respiration and mitochondrial content.⁴⁵ P4HB encodes an enzyme involved in the hydroxylation of prolyl residues in preprocollagen, which serves to stabilize the collagen triple helices.⁴⁶ Our findings suggest that some nuclear genes, due to their involvement in ECM organization and collagen synthesis or focal adhesion pathways, which were previously indicated as being deregulated in Polish KTCN corneas based on DNA and RNA studies,^20,37,47 play a role in disease pathogenesis. However, the direct function of these genes and their impact on mitochondrial function or oxidative stress in KTCN remain to be elucidated. Of note, we did not detect differences in the expression of three mitochondria-related genes: TFAM, TFB2M, and POLRMT, previously indicated as being deregulated in KTCN by Hao et al.⁷ 

In addition to expression analyses, we also re-analyzed our previous RRBS data²² obtained for KTCN and non-KTCN corneas to determine whether the methylation of 1233 nuclear-encoded mitochondria-related genes plays a role in KTCN. As no differentially methylated genes were found among the assessed mitochondrial–nuclear gene sets, we do not postulate a role for methylation of mitochondria-related genes in the pathogenesis of KTCN. 

Oxidative stress may lead to an increased number of changes in the mtDNA, and previously alterations in the mitochondrial genome in KTCN patients were investigated.^11,48,49 Here, 52 of the observed mtDNA substitutions within protein-coding regions were predicted as potentially pathogenic based on at least one applied prediction tool. The highest impact on protein structure and function was predicted for the rs2298008 (m.8921G>A). The m.8921G>A results in G132D substitution in the mitochondrially encoded adenosine triphosphate (ATP) synthase membrane subunit 6. However, this SNP was observed in the KC24 sample only, and there are no data available on the impact of rs2298008 on phenotype. 

Previously, 11 of the observed non-synonymous substitutions were reported as being involved in various diseases, including seven variants (m.4216T>C, m.13708G>A, m.4917A>G, m.386-6T>C, m.6261G>A, m.8836A>G, and m.9016A>G) that were described in patients with Leber hereditary optic neuropathy (LHON) disease.^50–54 One of these LHON-associated variants, m.4216T>C (rs1599988),^52,55 was detected in eight KTCN and 11 non-KTCN individuals, but none of them was diagnosed with LHON during ophthalmic evaluation. The relatively frequent occurrence of m.4216T>C in our dataset excludes the involvement of this variant in KTCN. The negligible pathogenicity of this variant in Polish patients with KTCN was further supported by in silico analyses. Whereas m.4216T>C leads to Y304H substitution in the transmembrane domain of the NADH dehydrogenase subunit 1, four of the applied prediction tools determined that this variant is benign. Previous studies have also suggested a slight effect of this variant on the LHON development.^49,53 As for the other LHON-associated substitutions detected in our study, m.13708G>A was observed in seven non-KTCN and four KTCN samples, and variant m.4917A>G, was equally distributed in both study groups (five in KTCN and five in non-KTCN samples). Other four LHON-associated variants only occurred in single samples. Based on those observations, we conclude that these variants also were not involved in the pathogenesis of KTCN in the evaluated patients. 

To further assess the role of identified mtDNA sequence variants, we compared their distribution in Sanger sequencing results and our previous RNA-seq investigation²⁰ in the blood–cornea matched pairs. This allowed us to examine whether the variants generally seen in blood samples were consistent with the changes observed in individuals’ corneal tissues. Comparison of two datasets revealed that mtDNA sequencing data (blood) and RNA-seq data (corneas) corresponded, indicating that there were no inter-tissue differences. Therefore, alterations in KTCN corneas due to oxidative stress, previously postulated as a causative factor in the KTCN research, might be caused by reasons other than mtDNA sequence variation. 

It has been reported that 10 KTCN patients (38.5%) had potentially pathogenic non-synonymous mtDNA mutations, including variants in genes coding complex I and various tRNAs.¹¹ Among them, nine non-synonymous variants were homoplasmic, and only one was heteroplasmic.¹¹ Potentially pathogenic variant m.4218T>A in the ND1 gene that results in premature stop-codon was absent in healthy controls and has a high level of heteroplasmy (74%).¹¹ Our NGS analysis of the patient with KTCN did not confirm a high level of heteroplasmy for sequence variants. However, because only one mtDNA sample was sequenced with NGS, it was not possible to evaluate potential heteroplasmic variants in the entire group of patients with KTCN. Additionally, we did not observe any changes in the expression of the ND1 gene in patients with KTCN compared with non-KTCN individuals. 

Abu-Amero et al.⁵⁶ detected 12 different haplogroups, but none of them correlated with the KTCN phenotype. However, further studies performed in 114 Saudi patients with KTCN and 552 healthy controls indicated that mitochondrial haplogroups H and R were significantly overrepresented in patients with KTCN when compared with the control group,⁵⁶ suggesting that these haplogroups may increase the risk of developing KTCN.⁵⁶ Furthermore, increased mtDNA content (copy number) in 119 patients with KTCN compared with 208 controls without KTCN has been detected.⁴⁸ In contrast, the analysis of KTCN patients from the Han Chinese population showed that there was no association between KTCN and mtDNA haplogroups H and R.⁴⁹ We detected 15 different mtDNA haplogroups in our study. No particular mtDNA haplogroup was present significantly more frequently in our KTCN samples. The H haplogroup, the most frequently observed in our study groups, is the most common in Europe and accounts for over 40% of mtDNA variation in anatomically modern humans in Western Eurasia.⁵⁷ 

Conclusions

In this study, no direct relationship was found among mitochondrial sequence variation, mitochondrial gene expression, mtDNA methylation, and KTCN phenotype. However, we identified abnormalities within nuclear genes regulating ECM formation, collagen synthesis, and/or focal adhesion that were previously reported as involved in KTCN in the same population and could be also related to mitochondrial function. These findings suggest that, although KTCN is not directly linked to mtDNA abnormalities, mitochondrial functioning may be partially affected in KTCN by changes in mitochondria-related nuclear gene expression. 

Acknowledgments

Supported by grants from the National Science Centre in Poland (2011/03/N/NZ5/01470 to DMN; 2012/05/E/NZ5/02127 to MG). 

Disclosure: D.M. Nowak-Malczewska, None; J.A. Karolak, None; J. Swierkowska, None; M.M. Jaworska, None; K.I. Kulinska, None; P. Polakowski, None; M. Rydzanicz, None; R. Ploski, None; J.P. Szaflik, None; M. Gajecka, None