Whole mitochondrial genome analysis in Chinese patients with keratoconus.

Purpose: Keratoconus (KC) is a corneal disorder characterized by corneal ectasia, progressive corneal thinning, and conical protrusion. This study aimed to elucidate the mitochondrial gene profile in Chinese patients with KC, analyze the mitochondrial haplogroup and heteroplasmy, and further explore the association between mitochondrial genes and KC.
Methods: Mitochondrial sequencing was conducted on 100 patients with KC and 100 matched controls. Haplogroup analysis was conducted with logistic regression analysis. The heteroplasmy was analyzed with ANOVA (ANOVA) and Student t test. Sequence kernel association tests (SKATs) were performed to analyze the association between mitochondrial genes and KC. Mtoolbox, Mitoclass.1, and APOGEE were used to estimate the impact of the identified variants in protein-coding genes. PON-mt-tRNA was used to annotate the impact of the variants in tRNA. RNAstructure was used to predict the secondary structures of native and mutated tRNAs.
Results: We identified 689 variants in patients with KC and 725 variants in controls (with 308 variants shared by both). The mitochondrial haplogroups exhibited no statistically significant differences between the two groups. Based on the heteroplasmy analysis, the number of heteroplasmic variants in the complete mitochondrial genome, RNA coding regions, and noncoding regions were statistically significantly different in the KC cases and controls (p<0.05). The heteroplasmic levels of the m.16180_16182delAA, m.16182insC, and m.14569 G>C variants in the KC cases were statistically significantly higher than those in the controls (p<0.05). The SKAT analysis showed that the COX3 and TRNH genes were statistically significantly associated with KC (p<0.05). Among the nine variants of COX3 included in the SKAT analysis (m.9300G>A, m.9316T>C, m.9327A>G, m.9355A>G, m.9468A>G, m.9612G>A, m.9804G>A, m.9957G>A, and m.9966 G>A), m.9612G>A was predicted to be deleterious by Mtoolbox. The m.9316T>C, m.9327A>G, m.9355A>G, m.9612G>A, m.9804G>A, and m.9957G>A variants were predicted to be damaging by Mitoclass.1. The m.9355A>G and m.9804G>A variants were predicted to be pathogenic by APOGEE. All identified variants located in TRNH (m.12153C>T, m.12178C>T, and m.12192G>A) were predicted to be neutral by the PON-mt-tRNA website. Conclusions: This study presents the mitochondrial gene profile of Chinese patients with KC and demonstrated that the COX3 and TRNH genes were associated with KC.

Introduction

Keratoconus (KC) is a corneal disorder characterized by corneal ectasia, progressive corneal thinning, and conical protrusion [1], with a reported prevalence of 1.38 in 1,000 individuals [2]. The disease is progressive and requires keratoplasty treatment at an advanced stage [3]. Several studies have indicated that heredity plays an important role in the process of KC [4,5]. Genetic studies on the nuclear genome have identified numerous loci associated with KC [4,6]. However, these loci explain only a small proportion of the heritability, and additional heritable factors have yet to be discovered. As circular DNA independent from the nucleus, the mitochondrial genome also plays a vital role in the pathogenesis of KC [7,8].

Currently, the molecular pathogenesis of KC is poorly understood. Recent studies suggested that oxidative stress could be involved in the pathogenesis of the disease [9,10]. Variations in mtDNA have the potential to alter mitochondrial function, affect the generation of reactive oxygen species (ROS), and lead to increased oxidative stress [11]. Abu-Amero et al. [12] found that mitochondrial sequences were altered in Saudi patients with KC, which indicated that mtDNA mutation potentially contributed to KC progression through the oxidative stress mechanism. Moreover, mitochondrial complex I gene analysis in Indian patients with KC demonstrated that nonsynonymous mtDNA sequence variations might account for elevated ROS production, leading to oxidative stress and finally, resulting in the pathogenesis of KC [13]. Hao et al. [14] found that KC corneas had increased mtDNA damage and higher ROS levels compared to normal corneas in Chinese populations. In addition, Hao et al. [15] revealed that the mtDNA copy number might contribute to the pathogenesis of KC. Nevertheless, the mitochondrial genome profile has not been reported in Chinese patients with KC. In addition, no study has comprehensively elucidated the relationship between mtDNA variants and risks of KC in the Chinese population.

In the present study, whole mitochondrial genome sequencing was used to elucidate the mitochondrial gene profile in Chinese patients with KC. Association analysis was also conducted to illustrate the relationship between mtDNA genes and the risks of KC.

Methods

Study population

A total of 100 patients with KC were recruited from Henan Eye Hospital & Henan Eye Institute. The diagnosis of KC was based on the following inclusion criteria: Corneal topography revealed an asymmetric bowtie pattern with or without skewed axes or signs of KC detected via a slit-lamp examination, such as localized stromal thinning, conical protrusion, Vogt’s striae, Fleischer’s ring, or anterior stromal scar [5]. Cases secondary to trauma, surgery, and other diseases were excluded from the study. The mean corneal curvature (steep keratometry, Ks) of all patients with KC (serious eye for each patient) was 56.3±8.30 D. One hundred matched individuals without any ocular disorder or previous ophthalmic surgeries were enrolled as controls. All the control subjects demonstrated a clear cornea on a slit-lamp examination, and a topographic map was within normal limits. The study adhered to the tenets of the Declaration of Helsinki and the ARVO statement on human subjects.

The study protocol was approved by the ethical committees of Henan Eye Hospital [ethical approval number: HNEECKY-2019 (5)]. Parents were informed of the study purpose for cases that were less than 18 years old. Written informed consent was obtained from all subjects.

Mitochondrial genome sequencing

EDTA anti-coagulated venous blood samples were collected from each subject in the current study. And the blood samples were preserved at −80 °C cryogenic refrigerator before use. Total DNA was extracted from peripheral blood samples with QIAamp DNA Blood kits (Qiagen, Hilden, Germany) according to the manufacturer’s recommendations. The DNA samples were quantified with NanoDrop 2000 and ascertained with electrophoresis.

The mitochondria genome was sequenced by Shanghai Genesky Biotechnologies. The mtDNA was amplified with six pairs of primers with long-range PCR. The six partially overlapping PCR products (about 6 kb each) were pooled and then purified with Agencourt AMPure XP-PCR Purification kits (Beckman, Germany). Then the PCR products were fragmented into 100–500 bp pieces with ultrasound (ME220, Woburn, MA). The fragmented DNA was subjected to end-repairing, dA-tailing, and adaptor ligation with the NEBNext® DNA Library Prep Reagent Set from Illumina® (New England Biolabs, Ipswich, MA). The adaptors ligated to each sample contained different indexes, which we subsequently used for pooled samples testing. Then the PCR reaction was performed in a total volume of 50 μl consisting of 15 μl genomic DNA with indexes, 5 μl pair end primers mix (Illumina, San Diego, CA), 25μl NEBNext UltrallQ5 Master mix (New England Biolabs) and 5 μl ddH2O. The cycling condition was 98 °C for 30 s, followed by 10 cycles of amplification at 98 °C for 10 s, 65 °C for 75 s, and a final extension at 65 °C for 60 s. After PCR amplification, the amplicons were size-checked and quantitated with BioAnalyzer 2100, and then subjected to 2 × 150 bp paired-end massively parallel sequencing using Illumina NovaSeq System (Illumina).

Sequencing data analysis and functional annotation

Fastq data were acquired for a preliminary quality check. Burroughs-Wheeler Aligner (BWA) [16] was used to align the sequence data to the mitochondrial reference sequence (revised Cambridge Reference Sequence, rCRS) and generate sequence alignment files. Variant calling was performed using the Genome Analysis Toolkit (GATK) Mutect2 and Haplotype Caller [17]. Then variants were filtered according to the threshold value (mapping quality ≥20 and base quality ≥20). Haplogroup analysis was conducted on https://haplogrep.i-med.ac.at/app/index.html. Variants exhibiting frequency values within the 0.1–0.9 range were considered heteroplasmic. Moreover, the mean number of heteroplasmic mtDNA variants was obtained by the average number of individual heteroplasmic variants present in the specified location from all samples within the group. The symbols and full names of the total 37 genes in human mitochondria are listed in Appendix 1. The schematic representation of the sequencing data analytical workflow is shown in Figure 1. For associated variants located in the protein-coding gene, Mtoolbox [18], Mitoclass.1 [19], and APOGEE [20] were applied to estimate the impact of the identified variants. For associated variants located in the tRNA gene, PON-mt-tRNA [21] was used to annotate the impact of the variants. Moreover, RNAstructure [22] was used to predict the secondary structures of native and mutated tRNAs.

Figure 1

Schematic representation of the sequencing data analytical workflow.

Statistical analyses

Statistical analyses were performed using the SPSS 21.0 package (SPSS Inc., Chicago, IL). Age was expressed as mean ± standard deviation and analyzed with Student t test, whereas gender was analyzed with the chi-square test. Haplogroup analysis was conducted through logistic regression analysis. The number of heteroplasmic variants was compared through ANOVA (ANOVA) and Student t test. The Mann–Whitney U test was used to compare the heteroplasmic level between KC cases and controls in variants with heteroplasmic samples ≥3 in each group. Sequence kernel association test (SKAT) was conducted with the software Efficient and Parallelizable Association Container Toolbox v.3.2.6 (EPACTS v.3.2.6) [23]. SKAT, which applies a multivariate regression model to uncover genes associated with a specific disease, could evaluate the cumulative effect of rare and common variants at the gene level [23]. Variants with a genotype missing rate greater than 0.1 were included in the SKAT. P values less than 0.05 were considered statistically significant.

Results

Spectrum of mitochondrial genome variants

Whole mitochondrial sequencing was performed in 100 patients with KC (73 men and 27 women, mean age of 19.7±3.50 years) and 100 controls (72 men and 28 women, mean age of 19.1±2.40 years) to explore the profile of the mitochondria genome and the role of mtDNA genes contributing to the susceptibility to KC. No statistically significant difference was found in the mean age and gender distribution between patients with KC and controls (all p>0.05). The sequencing data were generated from all subjects with a mean read depth of 3,139.50±524.160 (range from 2,122.15 to 6,140.10). As shown in Figure 2A, 368 variants were shared between patients with KC and controls, 321 were found only in patients with KC, and 357 were found only in controls. The distribution of the variants in the protein-coding genes identified in the study population is shown in Figure 2B. The results showed that the ND5 (OMIM:516005; Gene ID:4540) gene contained 107 variants, followed by the CYTB (OMIM:516020; Gene ID:4519), COX1 (OMIM:516030; Gene ID:4512), and ND4 genes. Furthermore, we categorized the variants in the protein-coding genes in patients with KC and controls separately. A total of 68 variants were located in the ND5 gene, followed by 54 variants in the CYTB gene in patients with KC. The ND5 gene contained 68 variants in the control group, followed by 57 variants in the CYTB gene (Figure 2D, Appendix 2).

Figure 2

mtDNA variants in patients with KC and controls. A: Venn diagram showing the number of identified variants in patients with keratoconus (KC) and controls. B: Number of identified variants in protein-coding genes identified in the study population. C: Number of synonymous and nonsynonymous variants identified in the study population. D: Number of synonymous and nonsynonymous variants in protein-coding genes identified in the study population. E: Number of synonymous and nonsynonymous variants in patients with KC and controls. F: Number of nonsynonymous variants in protein-coding genes identified in patients with KC and controls.

Distribution of synonymous and nonsynonymous variants

When we categorized all the variants according to their effects, 462 synonymous variants were identified, while there were 201 nonsynonymous variants (Figure 2C). Among all variants in protein-coding genes, there were more synonymous variants than nonsynonymous variants except the ATP6 (OMIM:516003; Gene ID: 4538) gene (Figure 2D). Among all the nonsynonymous variants, 135 variants were identified in patients with KC, whereas 130 were identified in controls (Figure 2E). No statistically significant difference was found by comparing the number of synonymous and nonsynonymous variants between the two groups (p=0.316). Afterward, we compared the number of nonsynonymous variants in protein-coding genes identified in the two groups. Of note, the ATP6, ATP8 (OMIM:516060; Gene ID: 4508), COX3 (OMIM:516050; Gene ID: 4514), ND2 (OMIM:516001; Gene ID: 4536), ND3 (OMIM: 516002; Gene ID: 4537), ND4L (OMIM:516004; Gene ID: 4539), ND5, ND6 (OMIM:516006; Gene ID: 4541) genes had higher number of nonsynonymous variants in KC patients, whereas the COX1, COX2 (OMIM:516040; Gene ID: 4513), CYTB, ND1 (OMIM:516000; Gene ID: 4535), and ND4 genes had lower number of nonsynonymous variants (Figure 2F, Appendix 3).

Mitochondrial haplogroup analysis

Mitochondrial haplogroups were analyzed in 100 patients with KC and 100 controls. The distribution of the mitochondrial haplogroups in the KC cases and controls is shown in Table 1. No statistically significant difference was found between KC patients and controls for all mitochondrial haplogroups tested (p>0.05).

Table 1

Distribution of mitochondrial haplogroups in KC cases and controls.

Haplogroup	Control (n,%)	KC (n,%)	P	OR (95%CI)
R	3(3%)	7(7%)	0.331	2.423(0.534–14.954)
C	5(5%)	2(2%)	0.445	0.389(0.036–2.449)
D	23(23%)	28(28%)	0.517	1.300(0.655–2.600)
G	11(11%)	8(8%)	0.631	0.705(0.234–2.027)
F	16(16%)	13(13%)	0.689	0.785(0.326–1.862)
M	15(15%)	17(17%)	0.847	1.160(0.508–2.674)
A	7(7%)	7(7%)	1	1.000(0.287–3.486)
B	9(9%)	9(9%)	1	1.000(0.335–2.987)
H	1(1%)	1(1%)	1	1.000(0.013–79.233)
N	4(4%)	3(3%)	1	0.743(0.106–4.521)
T	2(2%)	1(1%)	1	0.497(0.008–9.684)
Z	4(4%)	4(4%)	1	1.000(0.181–5.533)

Mitochondrial heteroplasmy analysis

For this analysis, heteroplasmic variants were compared according to different mitochondrial regions in the KC cases and controls (Figure 3). Among the protein-coding regions, RNA coding regions, and noncoding regions, the noncoding regions had the highest number of heteroplasmic variants in the two groups. By comparing the number of heteroplasmic variants in different mitochondrial regions, the average number of total heteroplasmic variants in the KC cases was higher than that in controls (p<0.05) in the complete mitochondrial genome. Regarding heteroplasmic variants in the protein-coding regions, no statistically significant differences was found between the KC cases and the controls (p>0.05). However, the number of heteroplasmic variants in the RNA coding regions was lower than that in the controls (p<0.05). The number of heteroplasmic variants in the noncoding regions was higher than that in the controls (p<0.05). After comparing the heteroplasmic levels between the KC cases and the controls, we found that the heteroplasmic levels of the m.16180_16182delAA, m.16182insC and m.14569 G>C variants in the KC cases were significantly higher than those in the controls (p<0.05, Appendix 4).

Figure 3

Comparison of heteroplasmic mtDNA variants in protein-coding regions, RNA coding regions, noncoding regions, and the complete mitochondrial genome in patients with keratoconus (KC) and controls. # denotes statistical significance among the protein-coding regions, RNA coding regions, and noncoding regions; * denotes statistical significance between the KC cases and the controls.

Association of mitochondrial genes with KC

Given the observation that synonymous mutant alleles had a neutral effect on protein function and disease susceptibility, we selected nonsynonymous variants in the protein-coding genes for further analysis. The SKAT was applied to assess the contribution of nonsynonymous variants and variants over rRNA and tRNA genes in regard to KC susceptibility. The results revealed that the COX3 and TRNH (OMIM:590040; Gene ID: 4564) genes were statistically significantly associated with KC (p<0.05; Table 2). Nine variants in the COX3 gene and three variants in the TRNH gene were included in the SKAT analysis. Functional annotations of these variants are listed in Table 3 and Table 4. Among the nine variants located in the COX3 gene (m.9300G>A, m.9316T>C, m.9327A>G, m.9355A>G, m.9468A>G, m.9612G>A, m.9804G>A, m.9957G>A, and m.9966 G>A), m.9612G>A was predicted to be deleterious by Mtoolbox. The m.9316T>C, m.9327A>G, m.9355A>G, m.9612G>A, m.9804G>A, and m.9957G>A variants were predicted to be damaging by Mitoclass.1. The m.9355A>G and m.9804G>A variants were predicted to be pathogenic by APOGEE. All three variants located in the TRNH gene (m.12153C>T, m.12178C>T, and m.12192G>A) were predicted to be neutral via the PON-mt-tRNA website. The secondary structures of native tRNAHis and mutated tRNAHis predicted with RNAstructure are presented in Figure 4.

Table 2

Association between mitochondrial genes and KC risks based on SKAT analysis.

Region type	Gene	SKAT
		N	P
protein-coding genes	ATP6	39	0.190
	ATP8	5	0.880
	COX1	19	0.270
	COX2	8	0.846
	COX3	9	0.023*
	CYTB	32	0.069
	ND1	19	0.636
	ND2	18	0.433
	ND3	5	0.711
	ND4	10	0.325
	ND4L	2	0.225
	ND5	29	0.357
	ND6	10	0.553
rRNA genes	RNR1	127	0.406
	RNR2	34	0.225
tRNA genes	TRNA	4	0.685
	TRNC	5	0.204
	TRND	4	0.834
	TRNE	4	0.71
	TRNF	2	0.479
	TRNG	4	0.157
	TRNH	3	0.040*
	TRNI	2	0.479
	TRNK	1	0.485
	TRNL1	2	0.689
	TRNL2	4	0.834
	TRNM	-	-
	TRNN	-	-
	TRNP	136	0.405
	TRNQ	1	0.485
	TRNR	3	0.499
	TRNS1	2	0.479
	TRNS2	1	0.485
	TRNT	9	0.815
	TRNV	1	0.485
	TRNW	1	0.839
	TRNY	-	-

SKAT: sequence kernel association test; N: number of variants included in SKAT analysis for each gene. *p<0.05 denotes statistical significance.

Table 3

Site information and functional annotation of variants in COX3 identified from SKAT analysis.

Gene	Position	Refer allele	Alter allele	Amino acid change	Mtoolbox	Mitoclass 1	APOGEE	Published disease association
COX3	9300	G	A	A32T	neutral	neutral	neutral	sudden unexpected infant death [36]
COX3	9316	T	C	F37S	neutral	damaging	neutral	-
COX3	9327	A	G	T41A	neutral	damaging	neutral	-
COX3	9355	A	G	N50S	neutral	damaging	pathogenic	-
COX3	9468	A	G	T88A	neutral	neutral	neutral	-
COX3	9612	G	A	V136M	deleterious	damaging	neutral	-
COX3	9804	G	A	A200T	neutral	damaging	pathogenic	LHON/ TOF/ MELAS/OAG [37–40]
COX3	9957	T	C	F251L	neutral	damaging	neutral	MELAS / NAION /CPEO/ HCM / gout [41–45]
COX3	9966	G	A	V254I	neutral	neutral	neutral	-

SKAT: sequence kernel association test; LHON: Leber hereditary optic neuropathy; TOF: tetralogy of Fallot; MELAS: mitochondrial encephalomyopathy, lactic acidosis and stroke-like episodes; OAG: open-angle glaucoma; NAION: non-arteritic ischemic optic neuropathy; CPEO: chronic progressive external ophthalmoplegia; HCM: hypertrophic cardiomyopathy.

Table 4

Site information and functional annotation of variants in TRNH identified from SKAT analysis.

Gene	Position	Refer variant	Alter variant	PON-mt-tRNA		Published disease association
				Score	Prediction
TRNH	12,153	C	T	0.09	neutral	thiamine deficiency accompanied mitochondrial myopathy [46]
TRNH	12,178	C	T	0.17	neutral	KC [12]
TRNH	12,192	G	A	0.12	neutral	CM/deafness [47,48]

SKAT: sequence kernel association test; KC: keratoconus; CM: cardiomyopathy.

Figure 4

Secondary structures of native tRNAHis and mutated tRNAHis predicted by RNAStructure. A: native tRNAHis; B: mutated tRNAHis of m.12153C>T; C: mutated tRNAHis of m.12178C>T; D: mutated tRNAHis of m.12192G>A.

Discussion

In the present study, we provided a comprehensive understanding of the mitochondrial genome profile in Chinese patients with KC. Haplogroup analysis of the KC cases and controls revealed no statistically significant difference. In addition, heteroplasmy analysis demonstrated that the number of heteroplasmic variants in RNA coding regions, noncoding regions, and complete mitochondrial genome were statistically significantly different between the two groups. Heteroplasmic level analyses indicated the heteroplasmic levels of the m.16180_16182delAA, m.16182insC, and m.14569 G>C variants in the KC cases were statistically significantly higher than those in the controls. Of note, the COX3 and TRNH genes were statistically significantly associated with KC according to the gene-based SKAT analysis.

Previous studies identified several mitochondrial genome variants in KC [12,13]. However, the mitochondrial genome profile in Chinese patients with KC was unclear, and no studies have comprehensively elucidated the relationship between mtDNA genes or variants and risks of KC in the Chinese population. Therefore, we conducted whole mitochondrial sequencing to elucidate the genetic profile of mitochondria in Chinese patients with KC and analyzed the association between mitochondrial genes or variants and KC. Among the variants identified in the KC cases, the ND5 gene had the highest number of variants in mitochondrial protein-coding genes, which was consistent with the findings of Pathak et al. [13]. Generally, there were more synonymous variants than nonsynonymous variants [24]. In this study, more synonymous variants were identified in both groups. Comparison of synonymous and nonsynonymous variants between the KC cases and the controls revealed no statistically significant difference. Nevertheless, Pathak et al. [13] found a higher number of nonsynonymous variants in patients with KC than in controls in their study. The controversial results could be attributed to the differences in populations and sample sizes.

Previous studies have identified an association between mitochondrial haplogroups and various diseases, including ophthalmic diseases such as glaucoma [25,26] and Leber hereditary optic neuropathy (LHON) [27]. In the present study, we evaluated the association between mitochondrial haplogroup and KC. However, no statistically significant differences were found for all mitochondrial haplogroups tested. Abu-Amero et al. [28] found that mitochondrial haplogroups H and R were statistically significantly associated with KC in Saudi patients with KC. Thus, the differences in populations and sample sizes potentially account for the controversial results.

At present, there is a paucity of information regarding the extent of mtDNA heteroplasmy in patients with KC and controls. We found a relatively high number of heteroplasmic mtDNA variants in the noncoding regions of both groups, which might be attributed to various factors, including the variable D-loop within noncoding regions and negative selective pressures associated with mtDNA variation in protein and RNA coding regions [29,30]. Moreover, the number of heteroplasmic variants in RNA coding regions and noncoding regions were statistically significantly different between the KC cases and controls, implicating that the RNA coding and noncoding regions may influence KC. In addition, we found the heteroplasmic levels of the m.16180_16182delAA, m.16182insC, and m.14569G>C variants were higher in the KC cases, indicating that those variants might have a potential role in the disease. However, the accurate molecular mechanisms remain unclear. Further exploration of the potential role of mtDNA heteroplasmy on KC is warranted.

Recent studies have shown that rare and common variants can contribute to the same disease [31,32]. SKAT, which applies a multivariate regression model to uncover genes associated with a specific disease, could evaluate the cumulative effect of rare and common variants at the gene level [23]. As synonymous variants have a neutral impact on protein function and disease susceptibility, we focused on nonsynonymous variants in mitochondrial protein-coding genes for further analysis. After analyzing the contribution of nonsynonymous variants in mitochondrial protein-coding genes and variants over rRNA and tRNA genes to KC susceptibility with SKAT at the gene level, we found the COX3 and TRNH genes may contribute to the occurrence of KC. The mitochondrial genome harbors 13 protein-coding genes, all of which encode the subunits of the electron transport chain, where oxidation phosphorylation (OXPHOS) occurs, and are responsible for most cell ATP production [33]. MtDNA variants often disrupt OXPHOS, compromising ATP production and releasing ROS. Cytochrome c oxidase III, encoded by the COX3 gene, is a subunit of respiratory complex IV. Notably, respiratory complex IV plays an essential role in the electron transport chain [33]. Therefore, we speculated that the abnormality of the COX3 gene might disrupt OXPHOS by altering the electron transport chain, consequently resulting in KC. The TRNH gene encodes tRNAHis which participates in protein translation. Gong et al. [34] constructed cybrids of a mutation in TRNH by transferring mitochondria from lymphoblastoid cell lines. After detecting the levels of mitochondrial ATP and membrane potential in mutant cybrids, Gong et al. found the mutant cybrids had increased ROS levels, inferring that the mutation in TRNH could lead to the increasing production of ROS. In addition, KC fibroblasts had increased basal generation levels of ROS [35]. Thus, we inferred that variants in the TRNH gene might elevate ROS levels, leading to the occurrence of KC.

In the present study, nine variants in COX3 and three variants in TRNH were included in the SKAT analysis. Among the nine variants located in the COX3 gene, the m.9300G>A variant was reported to have a potential role in some patients with sudden unexpected infant death [36]. The m.9804G>A variant has been reported in LHON, tetralogy of Fallot (TOF), mitochondrial encephalomyopathy, lactic acidosis and stroke-like episodes (MELAS), and open-angle glaucoma (OAG) [37-40]. The m.9957G>A variant has been reported in cases of MELAS, non-arteritic ischemic optic neuropathy (NAION), chronic progressive external ophthalmoplegia (CPEO), hypertrophic cardiomyopathy (HCM), and gout [41-45]. Among the variants located in the TRNH gene, m.12153C>T was reported in patients with thiamine deficiency–accompanied mitochondrial myopathy [46], m.12178C>T was previously reported in Saudi patients with KC [12], and m.12192G>A was reported in cardiomyopathy (CM) and deafness [47,48]. The predicted secondary structures of m.12153C>T and m.12192G>A were similar to native tRNAHis, whereas the predicted structure of m.12178C>T was statistically significantly different from native tRNAHis. This discrepancy of predicted structures indicates that the m.12178C>T variant might critically influence the occurrence of KC. Previously, Abu-Amero et al. [12] screened the mitochondrial genome in Saudi patients with KC and found nine potentially pathogenic variants in the study population. However, only m.12178C>T was identified in the present study, which is in accordance with Abu-Amero et al. findings. Of note, this might be attributed to different populations and different sample sizes, which warrants future studies with multicenter and larger population sizes for in-depth exploration.

In conclusion, this study presented the mitochondrial gene profile of Chinese patients with KC and identified that the COX3 and TRNH genes were associated with KC.