Identification And Bioinformatics Analysis Of The M~6A-modification Differentially-expressed Genes In Keratoconus

Objective:Keratoconus（KC）is a progressive central thinning and conical protrusion of the cornea.It is one of the most common blinding corneal diseases in adolescents and the fourth indication of corneal transplantation in China and the United States.Currently,there is no medical treatment for KC in clinic,but only corneal transplantation can be available in the advanced stage.Profoundly understanding the pathogenesis of KC is therefore crucial for clinical prevention and treatment.At present,although the pathogenesis of KC mainly involves in inflammation,oxidative stress and biomechanics,the exact mechanism still remains largely unknown.N⁶-methyladenosine（m⁶A）modification is one of the most abundant epigenetic modifications in eukaryotic RNAs,and play critical roles in regulating RNA splicing,expression,decay,translation and stability.m⁶A modification reportedly is closely associated with various physiological and pathological conditions,suggesting the possible role of m⁶A modification in the KC pathogenesis.Therefore,we aimed to investigate whether m⁶A modification plays potential role in the development of KC.Methods:1.Clinical samples collectionEleven KC patients（11 eyes）undergoing corneal transplantation at the Eye Hospital of Shandong First Medical University were selected as the experimental group（KC）,and11 healthy corneas from eye bank donors as the control group（CON）.Among them,6pairs were used to determine the expression of m⁶A-modified associated proteins,and 5pairs were applied to perform m⁶A microarray analysis.2.Determination of m⁶A-modification associated proteinsWestern blot was used to examine the expression of m⁶A-modified proteins in 6 KC and 6 normal corneas,respectively.The m⁶A-modified proteins mainly included Wilms tumor-associated proteins（WTAP）,methyltransferase-like 3,Methyltransferase-like 14（METTL14）,Fat mass and obesity-related proteins（FTO）,Alk B Homolog 5（Alk B Homolog 5,ALKBH5）,YTH Domain family（YTHDF2）.3.RNA extraction and m⁶A microarray analysisTotal RNA was extracted from KC and CON samples using TRIzol method,and m⁶A-mrna epigenetic transcriptomic microarray analysis was then performed to construct m⁶A modification profiles in CON and KC corneas.m⁶A modification data on m RNA were thus obtained.4.The functional prediction of m⁶A-modified differentially-expressed genes using bioinformaticsThe differentially-expressed genes（DEGs）were selected based on the the threshold of m⁶A modification standard（Fold change≥2,p-value≤0.05）.The identified m6A-modification associated DEGs were then used to perform GO（Gene Ontology）and KEGG（Kyoto Encyclopedia of Genes and Genomes）analysis,to preliminarily clarify the role of DEGs in the pathogenesis of KC.5.The validation of m6A-modified DEGs by Me RIP-q PCRTo determine the reliability of m6A microaray data,the candidate DEGs were randomly screened out,and then verified via the combination of methylated RNA-protein immunoprecipitation and real-time fluorescence quantitative PCR（Me RIP-q PCR）.Results:1.The diordered expression of m⁶A-modification-related proteins in KC tissuesWhen compared CON group,the protein levels of m⁶A reader YTHDF2,m⁶A demethylase（Eraser）FTO and methyltransferase（Writer）WTAP in KC samples were significantly downregulated.However,the levels of methyltranferase METTL3 and METTL14 in KC samples remained no pronounced alteration.These findings revealed the abnormality of m⁶A modification related protiens in KC,indicative of the possible involvement of m⁶A modification in the pathogenesis of KC.2.RNA quality controlAfter extraction by Trizol mothed,the OD_260/280ratio of 1.8-2.1 and OD_260/230ratio of1.8-2.1 were used to evaluate RNA quality.All RNA samples achieved the quality control standard,paving the way for the subsequent m⁶A-m RNA epitranscriptomic microarray analysis.3.Profiles of m⁶A-modified DEGs in KCHierarchical cluster analysis showed that the transcriptional expression pattern of KC samples was quite different from the CON samples.A total of 17,297 transcripts were identified through m⁶A-m RNA epitranscriptomic chip,and 366 genes with differentially m⁶A modification（Fold change≥2,P-value≤0.05）were obtained,with 287hypermethylated genes and 79 hypomethylated.Through the four-quadrant diagram,we found 527 transcripts with hypermethylation and upregulated transcriptional levels,241transcripts with hypomethylation and downreulated transcriptional levels,and 49transcripts with hypermethylation but lowered RNA levels.4.Functional analysis of m RNA regulated by m⁶A modification in KCWe performed bioinformatic function analysis of 366 genes with differential m6A modifications correspoinding to 817 transcripts（|log2FC|≥2）.Firstly,Go analysis revealed that the biological processes（BPs）of hypermethylated genes involved mainly included the regulation of cardiac contractility,positive regulation of action potential,positive regulation of ATPASE activity,positive regulation of protein localization to the center of microtubule tissue and positive regulation of striated muscle contraction,and so on.While the BPs enriched in hypomethylated genes mainly covered pri-mi RNA transcribed by RNA polymerase II,and aortic development,etc.Moreover,KEGG analysis showed that the hypermethylated genes were mainly involved in porphyrin and chlorophyll metabolism,fluid shear stress and atherosclerosis,VEGF signaling pathway and c GMP-PKG signaling pathway,but the hyppmethylated genes were mainly associated with neuroactive ligands-receptor interaction,glutamate synapse and other signaling pathways.5.The validation of candidate m⁶A modified genes by Me RIP-q PCRA total of 11 candidate genes were verified by MERIP-QPCR.Among them,Cathepsin B（CTSB）and Cereblon,ATP-dependent Lon,CRBN,Zinc finger protein ZBTB38,Leucine receptor protein（SAR1B）,Abhydrolase Domain Containing 16A,ABHD16A)and acyl-Co A Binding Domain Containing 4（ACBD4）in KC samples were presented with increased methylation level and m RNA transcription level,when compared with the CON samples.By contrast,HMGCR,ADAMTS6 and DLGAP1showed hypomethylation and lower transcriptional levels in KC samples than in CON samples.The trends of the transcriptional and m⁶A modification levels of these candidates were consistent with the results obtained from m⁶A microarray,which consolidated the rliability of m⁶A chip data.Conclusion:These findings revealed the disordered expression of m⁶A-modificaton related proteins in KC samples,and a total of 336 genes with siginificant m⁶A modification were identified using m⁶A-m RNA epitranscriptomic microarray.Functional analysis revealed that those DEGs were closely associated with biomechanics,neuromodulation,and pri-mi RNA production.Overall,these findings provided the alternative clues for further investigating the pathogenesis of KC.