Molecular Vision 2020; 26:76-90
Received 22 May 2019 | Accepted 21 February 2020 | Published 25 February 2020
Last two authors contributed equally to this work.
1Department of Ophthalmology & Visual Science, Eye & ENT Hospital, Shanghai Medical College, Fudan University, Shanghai, China; 2Bioinformatics, School of Life Sciences and Biotechnology, Tongji University; 3Key Laboratory of Myopia, NHFPC (Fudan University), Shanghai, China; 4Shanghai Key Laboratory of Visual Impairment and Restoration (Fudan University), Shanghai, China; 5State Key Laboratory of Medical Neurobiology, Institutes of Brain Science and Collaborative Innovation Center for Brain Science, Fudan University, Shanghai, China
Correspondence to: Junyi Chen, Department of Ophthalmology, 83 Fenyang Road, Shanghai, China; email: firstname.lastname@example.org
Purpose: This study investigates the impact of aging on the miRNA expression profile in porcine angular aqueous plexus (AAP) cells, which are the porcine equivalent of human Schlemm’s canal endothelial cells.
Methods: AAP endothelial cells were isolated and cultured in physiologic (5% O2) or hyperoxic condition (40% O2) for 14 days to induce cell senescence. miRNA and protein expression profiles of control and senescent cells were analyzed with miRNA microarray and isobaric tags for relative and absolute quantification (iTRAQ), respectively.
Results: The miRNA microarray identified 33 differentially expressed miRNAs in senescent cells compared with controls (p<0.05), and quantitative real-time PCR (qRT-PCR) confirmed 12 of them (p<0.05). iTRAQ analysis identified 148 upregulated and 222 downregulated proteins (p<0.05, fold change>1.2). Bioinformatics analysis of miRNA microarray and proteomics data predicted that six out of seven miRNAs are associated with aqueous humor outflow by targeting integrin and the downstream pathways (Src/Rho kinase, focal adhesion kinase (FAK)/NO-cGMP), and one miRNA might influence gap junction by targeting the Inositol trisphosphate receptor (IP3R) /Protein kinase C (PKC) pathway.
Conclusions: This study identified miRNAs in senescent AAP cells that might regulate aqueous humor outflow by targeting proteins involved in focal adhesion, cytoskeleton, NO-cGMP signaling, and gap junction.
Aging is a strong risk factor for developing primary open angle glaucoma (POAG) [1-7]. Schlemm’s canal (SC) inner wall and juxtacanalicular tissue are the main sites of aqueous outflow resistance  which are responsible for the intraocular pressure (IOP) elevation in patients with POAG .
Senescence of SC and trabecular meshwork (TM) cells have been implicated in the initiation and progression of glaucoma [9-13]. The oxidizing species produced and accumulating during cell metabolism can cause molecular damage to the TM cells, which could increase aqueous humor outflow resistance [10,14]. Further, signiﬁcantly higher levels of DNA oxidation products were observed in the TM of patients with glaucoma, and the level of DNA damage strongly correlated with IOP elevation and visual field defects [13,15].
Hyperoxia leading to senescence is an established model of aging which has been widely used in aging research [16-18]. This model has also been investigated in the cells of the outflow pathway in the context of glaucoma pathogenesis [10,14,19,20]. This model is based on one of the most accepted theories of aging, stress-induced premature senescence (SIPS) . This model provides a constant increase in reactive oxygen species without the need for additional chemical or cellular treatments and has relatively low toxicity in that sustained confluence, and thus, normal post-mitotic conditions were assured throughout the length of the experiment [22-24].
Our previous studies showed that cell senescence influenced cell monolayer permeability through cytoskeletal proteins and cell adhesion proteins in porcine angular aqueous plexus (AAP) cells . Moreover, we demonstrated that senescence reduced the mechanotransduction sensitivity of AAP cells; specifically, the senescent cells were less able to respond to shear stress which should downregulate junctional proteins to facilitate aqueous humor drainage and maintain IOP homeostasis . Cell senescence also altered eNOS and phosphorylated eNOS, and resulted in reduced NO production . However, the mechanism by which these responses were regulated remains unknown.
MicroRNAs (miRNAs) are highly conserved non-coding RNAs. They typically negatively regulate protein expression  by binding to their target sites in the 3′-untranslated region (UTR) of the mRNA. miRNAs are key modulators of cellular senescence. miRNAs are differentially expressed in senescent cells, which further implicates them in the implementation of the senescent phenotype [28-30]. The dysregulation of miRNA-governed senescence underlies age-associated diseases [31-33] and cancer . Researchers reported the miRNA profile in the aqueous humor of patients with glaucoma and identified several differentially expressed miRNAs [35-37]. Further, microRNAs were differentially expressed in the retinas of eyes with advanced glaucomatous damage compared with normal controls . This evidence suggests that miRNA may be important in the pathogenesis of glaucoma and merits further investigation. However, the miRNA expression profile is not known in senescent Schlemm’s canal endothelial cells.
This study aimed to characterize miRNA expression profile in senescent AAP cells, which are the porcine equivalent of human SC cells. These results were further correlated to the protein expression profile with isobaric tags for relative and absolute quantification (iTRAQ) analysis to identify their target proteins. Bioinformatics analysis revealed the differentially expressed miRNAs may regulate critical pathways involved in aqueous humor outflow.
This study adhered to the ARVO Statement for Use of Animals in Research. And our research was proved by Laboratory Animal Management and Ethic Committee of Eye & ENT Hospital. Porcine AAP cells were isolated and cultured. Cell senescence was induced by hyperoxia for 2 weeks. miRNA microarray and proteomics analyses of normal and senescent AAP cells were performed to identify differential miRNAs and proteins. miRNA and protein data were correlated with bioinformatics analysis. The workflow of this study is shown in Figure 1.
Cell culture AAP endothelial cells were isolated, cultured, and characterized according to an established method . Briefly, for each cell line, outflow tissue from ten porcine eyes (4- to 6-month-old pigs) was collected. The cells were allowed to multiply for 8 days and then treated with puromycin (4 µg/ml, InvivoGen, San Diego, CA) for 2 days. The puromycin-selected cells were similar to human Schlemm’s canal endothelial cells in that the cells exhibited contact inhibition and expressed the same surface markers as those seen in cultured human SC cells and whole porcine tissue (Appendix 1). Three independent cell lines (each cell line came from ten porcine eyes) were used as one biological replicate for microarray assay or iTRAQ analysis.
Cell senescence was induced by a model of chronic oxidative stress [20,25]. Cells were cultured under normobaric hyperoxia condition (40% O2, 5% CO2) for 14 days in a triple-gas incubator (Smart cell, Shanghai, China). Control cultures were grown under physiologic oxygen condition (20% O2, 5%CO2) for 14 days in parallel with the experimental group of the same cell line. AAP cells were stained positive for endothelial cell specific markers VE-cadherin and eNOS, and senescence of AAP cells was confirmed with β-galactosidase staining (Appendix 1). The expression levels of p53 and p21 and cell cycle analysis were also investigated to ensure cell senescence (Appendix 1).
The miRNA microarray analysis was performed three times using three independent cell samples (each individual cell line came from ten porcine eyes). Total RNA was isolated with an miRNeasy kit (QIAGEN, Shanghai, China). The RNA integrity number (RIN) was assessed using a Bioanalyzer (Agilent 2100, Shanghai, China) according to the manufacturer’s instructions. It allows visual inspection of RNA integrity to determine the RNA quality. Only samples that achieved an RIN score greater than 9.5 were used for further microarray analysis. The miRNA microarray analysis was performed three times. In each one of the biological replicates, miRNA expression was measured using eight technical replicates.
The miRNA microarray analyses were performed using different samples by LC Sciences (Houston, TX). (The sequences are listed in Appendix 2) Four to eight micrograms of total RNA samples were fluorescence labeled and hybridized overnight. Fluorescence signals were collected using a laser scanner (GenePix 4000B, Molecular Devices, San Jose, CA) and digitized (Array-Pro, Media Cybernetics, Rockville, MD). After the background signals were subtracted, data were normalized using a locally weighted scatterplot smoothing (LOWESS) filter, which is a preferred method for miRNA data sets . Signals between groups were compared with the Student t test. To include as many miRNA candidates as possible, a p value of less than 0.05 without fold-change limitation was regarded as statistically significant.
Quantitative real-time PCR (qRT-PCR) of differentially expressed miRNAs was performed using qTOWER 2.2 (Analytik Jena, Jena, Germany) to validate differentially expressed miRNAs (p<0.05) to ensure statistical significance. RNA was extracted from three experimental samples and three control samples which were different from those used in the microarray analysis. cDNAs were synthesized in 40 μl reaction volume. The primer sequences are listed in Appendix 3. Amplification and detection were performed on 96-well plates. Each 10 μl Bestar™Real time PCR Master Mix (DBI, München, Germany) contained 5 μl 2×SYBR® Green Supermix, 1 μl cDNA reaction mixture, 0.5 μl reverse and sense primers, as well as 3 μl distilled deionized water (ddH2O). The reaction conditions were as follows: initial denaturation at 95 °C for 3 min, followed by 40 cycles of 95 °C for 10 s, 58 °C for 30 s, and finally 60 °C~95 °C,+1 °C/cycle holding time 4 s for melting curve analysis. The expression levels of mRNA were normalized to reference genes U6. Relative miRNA levels were calculated using the Pfaffl method  with the installed software qPCRsoft3.0. (The Pfaffl method is preferred for large differences between the amplification efficiency of the target gene and the internal gene.)
iTRAQ technology is a widely used method in the field of quantitative proteomics with high accuracy and reliability. It is an isobaric labeling method by tandem mass spectrometry to determine the number of proteins from different sources in a single experiment [42,43]. Proteomics of control cells and senescent cells were performed with iTRAQ, which detected the whole proteome in porcine proteins. iTRAQ analysis was performed three times.
Cell lysates were processed with radioimmunoprecipitation assay (RIPA) solution. Then, the protein concentration was measured using the Bradford method. One hundred micrograms of protein from each sample was processed for iTRAQ labeling. Then, the proteins were denatured, reduced, alkylated, trypsin digested, and labeled with iTRAQ labeling (8-plex iTRAQ, AB SCIEX, Framingham, MA). The pooled iTRAQ-labeled peptides were fractionated with strong cation exchange (SCX) chromatography (Shimadzu LC-20AB HPLC Pump system) using an SCX column containing 5 mm particles (4.66250 mm Ultremex column, Phenomenex, Torrance, CA).
The eluted peptides were pooled as 20 fractions, desalted with the Strata X C18 column (Phenomenex), and vacuum-dried. The elute from the Nano Liquid Chromatography system (Shimadzu LC-20AD) was coupled to a tandem mass spectrometer in an LTQ Orbitrap Velos (ThermoFisher Scientific, San Jose, CA), through an electrospray ionization source equipped with a 15 µm ID emitter tip.
Similar to previous publications ([44,45]), protein expression ratios of >1.2 or <0.83 and a p value of less than 0.05 according to the Student t test were considered differential proteins. The abundance ratio was calculated using the unique peptide strength corresponding to each protein (the sum of the ion strength of the identification spectrum label) to calculate the average ratio. Protein coverage is the length of the sequence identified in the protein divided by the length of the total protein sequence. If the unique peptide segment identified in the protein met the false discovery rate (FDR) filtering criteria, the protein was identified. Complete experimental details for iTRAQ are in Appendix 4.
Differentially expressed miRNAs of the control and experimental groups were identified with the t test using TMev bioinformatics software. To get the gene targets of differentially expressed miRNAs, porcine 3′-UTR sequences of mRNA were downloaded from the Ensemble database and complementarily linked with miRNA sequences. The final predicted targets and their corresponding proteins were determined by the intersection of results from three different databases (TargetScan, PicTa, and miRanda), which makes the prediction more reliable and highly accurate. Clustering analysis was routinely performed using the hierarchical method when the p value was less than 0.01. Average linkage and Euclidean distance metric were performed in cluster analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations were performed against porcine proteins using the DAVID gene annotation tool.
Bioinformatics analysis on the raw tandem mass spectrometry (MS/MS) data was performed according to a standard protocol (Appendix 4). Taking the amino acid sequences as input, differentially expressed proteins (iTRAQ) were further annotated with GO and Cluster of Orthologous Groups (COG) analysis. GO functional classifications were analyzed with Blast2GO software, while COG information was retrieved by blasting the sequences on the COG database . To explore the biological meanings of these proteins, GO enrichment analysis was performed to identify GO terms that were statistically significantly enriched in differentially expressed proteins (p<0.05), and the enriched metabolic pathways of these identified proteins were screened via the KEGG database (p<0.05).
Target proteins of differentially expressed miRNAs were identified from the iTRAQ analysis. These differentially expressed proteins were further mapped to KEGG pathways, if more than three proteins are related to known pathways that are important to the pathogenesis of glaucoma.
Cell lysates were prepared according standard protocols. Briefly, cells were prepared using RIPA solution. After protein concentrations were measured with the Bradford method, equal amounts of protein (30 μg protein/lane) were separated with 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE). Then the resolved proteins were transferred to nitrocellulose ﬁlters which was then blocked with 5% nonfat dry milk in Tris-buffered saline with 0.05% Tween-20 for 2 h. Filters were probed using primary antibodies, ITGAV (1:1,000, Abcam), ITGB3 (1:1,000, Abcam) pMLC (1:1,000), sGC (1:3,000), and eNOS (1:1,000), followed by incubation with peroxidase-linked secondary antibodies. GAPDH was used as a loading control. Four independent replicates were used for statistical analysis with the Student t test.
We identified differentially expressed miRNAs in senescent AAP cells with miRNA microarray analysis. Among the 407 miRNAs identified, 102 miRNAs were statistically significantly expressed in senescent cells compared to controls (p<0.05). The microarray raw data are shown in Appendix 5. The miRNA microarray hybrid signal ranged from 0 to 51,949; those that had a value greater than 500 were further analyzed (in our experience, a low hybrid signal was more likely to result in false positives). The heat map of the miRNAs shown in Figure 2 demonstrated 17 statistically significantly downregulated miRNAs and 16 statistically significantly upregulated miRNAs (p<0.05, signal>500, n=3, Student t test; detailed information is listed in Appendix 6). Then differential miRNAs were validated with qRT-PCR, which confirmed 12 miRNAs (Table 1).
Target genes and proteins of differentially expressed miRNAs are listed in Appendix 7. GO analysis revealed 504 GO terms with p values smaller than 0.05 with the submission of predicted target genes. GO terms of top-ranked classification gene numbers are listed in Appendix 8. Some are relevant to aqueous outflow, for example, extracellular matrix organization (GO:0030198, p=0.0111, gene number: 33), cell–cell junction organization (GO:0045216, p=0.0095, gene number: 10), focal adhesion (GO:0005925, p<0.0001, gene number: 139), cytoskeleton (GO:0005856, p=0.0056, gene number: 88), cell–cell junction (GO:0005911, p=0.0043, gene number: 65), microtubule cytoskeleton (GO:0015630, p=0.0436, gene number: 52), adherents junction (GO:0005912, p=0.0268, gene number: 14), and actin binding (GO:0003779, p=0.0316, gene number: 92).
KEGG analysis identified 138 pathways with p values smaller than 0.05. The top 20 KEGG pathways ranked by predicted gene number are listed in Appendix 8. Some are possibly involved in aqueous humor drainage.
iTRAQ analysis identified a total of 4,165 proteins and 17,603 unique peptides. The top ten up- and down-expressed proteins are presented in Table 2, and raw data are shown in Appendix 9 and Appendix 10. The protein abundance distribution is shown in Appendix 11. A total of 370 proteins were differentially expressed (p<0.05, fold change>1.2 or <0.83), of which 148 proteins were upregulated, and 222 proteins were downregulated. Protein fold-difference was expressed in log form with the base equaling 2. A log ratio greater than 0 indicates upregulation; less than 0 is downregulation. Most of the values fell within −5 and 5 (Figure 3). GO classification and enrichment analysis revealed the functional groups of differentially expressed proteins (Figure 4, Appendix 8). Intermediate filament (2.0%), extracellular matrix (3.8%), and cell–cell contact zone (1.5%) were identified. The KEGG pathway enrichment analysis is shown in Appendix 8. The COG protein analysis showed the protein distribution in different cell functional groups (Figure 5).
Differentially expressed miRNA and protein pairs are presented in Table 3. A total of seven miRNAs and 13 proteins were further mapped to KEGG pathways after selection. These miRNAs were computationally predicted to negatively or positively regulate their target proteins (Figure 6). Six miRNAs were predicted to directly target ITGA/ITGB, the two subunits of integrin, which could regulate IOP through the downstream pathways related to focal adhesion, cytoskeleton, NO-cGMP, and gap junction. Western blot (WB) verification of ITGAV and ITGB3 show no statistically significant changes in Figure 7. The possible reasons are analyzed in the Discussion section. However, the proteins (pMLC, eNOS, and sGC) of the downstream pathways showed a statistically significant difference.
This study identified a range of differentially expressed miRNAs and proteins in senescent AAP cells, which are the porcine equivalent of SC cells. The senescence of SC cells was thought to contribute to the elevation of IOP [11,12]. By correlating the miRNA microarray data with proteomics data, we identified those key miRNAs whose target proteins are also statistically significantly altered in AAP cells. The integration and analysis of differential miRNA and proteomics data may help us to better understand the responses occurring in SC cells under senescence.
miRNA microarray analysis identified 33 differentially expressed miRNAs (Appendix 6). However, the results of the miRNA microarray assays presented an inconsistent pattern within groups. Our speculation is that the difference between primary cells might be the possible reason for the inconsistency within groups. Although the variation is large, which may lead to omission, the detected miRNAs were unaffected. In addition, these miRNAs were confirmed with PCR to ensure a statistically significant difference. Twelve were verified with PCR. All 12 miRNAs were previously found to be altered or participate in cell senescence. Three miRNAs (miR-146a , miR-146b, and miR-15b ) were also detected in senescent TM cells. The other nine miRNAs (miR-34a , miR-99a , miR-24–3p , miR-181a , miR-23a , let-7g , let-7d-5p , miR-129a-3p , and miR-184 ) were also found to be involved in senescence in other cell types.
The total differential expressed miRNAs produced a total of 504 GO terms and 138 KEGG pathways. Some were identified to regulate IOP by actin binding, cytoskeleton, actin binding, cell–cell junction, extracellular matrix (ECM) organization, focal adhesion, and so on. In SC cells, stiffness  and junctions are related to the formation of pores. Reorganization of extracellular matrix  alters basement membrane permeability. The present findings are consistent with those of previous studies. For example, the let-7 family is believed to inhibit fibrosis by repressing expression of collagen genes . miR-24-3p accelerated the migration and invasion of bladder cancer cells . miR-15b directly targets tissue inhibitor of metallopeptidases 2 (TIMP2) to increase the migration and invasion of human lung cancer cells . However, overexpression of miR-15b aggravates IL-1beta-induced ECM degradation  in nucleus pulpous cells. In other words, these miRNAs were found to be involved in ECM regulation and cell junctions. The differential expression of 12 miRNAs was confirmed with PCR (Table 1). miRNA-protein pairs were identified for these miRNAs. A total of 4,165 proteins and 17,603 unique peptides were identified. GO and COG protein analysis further classified these differential proteins according to their functions. Binding, catalytic activity, and transporter activity were included in top molecular functions. The top biological processes were the cellular process, metabolic process, and single-organism process. In senescent cells, these top changes in molecular function and biological processes were reasonable and consistent with the cellular senescence hallmark [64,65]. In the senescent endothelial cells, signaling, biological adhesion, and localization might relate to alteration of the endothelial cell barrier and interaction with ECM.
Bioinformatic analysis further revealed that the seven miRNAs might regulate aqueous humor outflow in AAP cells through three signaling pathways which are important to glaucoma pathogenesis (Figure 6). The three pathways form a network through integrin and PKG. In the following, we discuss each signaling pathway in detail.
First, miR-23a, miR-146a-5p, miR-146b, and let-7g could regulate aqueous humor outflow by integrin/Src/Rho kinase/MLC/F-actin. ITGA and ITGB are direct targets of miR-23a, miR-146a-5p, miR-146b, and let-7g. Upon engagement of integrin receptors with extracellular ligands, the focal adhesion kinase (FAK)-Src complex is activated. Src transiently inhibits RhoA activity through RhoGAP and GRLF1 (p190RhoGAP) [66,67]. RhoA and Rac stimulation can increase the transendothelial resistance of human SC cells, by phosphorylating MLC , and MLC phosphorylation can trigger F-actin organization which impacts focal adhesion and cell permeability . It has also been reported that Rho-associated protein kinase inhibitor significantly decreases outflow facility by targeting junction or calcium ion transport of SC cells [70,71]. Rho kinase and MLC are important pharmaceutical targets for ocular hypertension. Rho kinase inhibitor (netarsudil, NDA:208254) was recently approved by the U.S. Food and Drug Administration (FDA) as a novel glaucoma drug to increase conventional outflow and treat ocular hypertension. MLC kinase inhibitor is also a potential new therapy for lowering IOP . It might be possible to use their upstream miRNA to regulate outflow through Rho kinase and MLC.
Second, let-7d-5p and miR-184 could regulate outflow by integrin, FAK, PI3K, Akt, eNOS, and NO-cGMP. The activation of ITGA and ITGB phosphorylates FAK, which is upstream of the phosphatidylinositol 3-kinase (PI3K)/Akt signal pathway . In a previous study, we showed that PI3K/Akt facilitates aqueous humor drainage through activation of eNOS , an enzyme that catalyzes the release of NO . NO-cGMP is an important pathway for IOP regulation. The NO-donating prostaglandin analog latanoprostene bunod (LBN) was developed as a novel drug for ocular hypertension .
Third, miR-24–3p could regulate aqueous humor outflow by IP3R, Ca2+, and PKC which indirectly controls gap junction. Gap junction is related to endothelial cell stiffness , which may affect IOP homeostasis potentially by vacuole and pore formation in SC cells [57,77]. Gap junction blocker (e.g., carbenoxolone) significantly increases outflow facility . Activation of IP3R by miR-24–3p results in the release of Ca2+ from the endoplasmic reticulum (ER) . The increase in the Ca2+ concentration activates connexins, a family of transmembrane proteins that constitute the gap junction channels . Connexin phosphorylation can change the unitary conductance and open probability of gap junctions .
It is interesting that six out of seven miRNAs directly target integrin and the expression of integrin. Integrin is a transmembrane heterodimer with two subunits (α and β). It serves as a ligand binding site with various downstream signaling in cytoplasmic side . Integrin is the binding site for structural adaptors (e.g., talin and tensin) which link directly to cytoskeleton . Integrin is also the binding site for scaffolding adaptors (e.g., kindlin and paxillin) which contribute to focal adhesion . As the binding site for catalytic adaptors, integrin interacts with Src and FAK which act as a signal transducer from adhesion molecules . ECM proteins are the primary ligands for integrin; integrin activation can enhance ECM deposition by integrin-mediated matrix assembly and therefore increase aqueous humor outflow resistance. However, an agonist of the integrin subtype α4β1 could increase cell–ECM detachment in SC cells and increase outflow facility .
In the present study and our previous study, although the downstream expression levels of proteins (pMLC, eNOS, and sGC) and cellular functions [11,12,26] were validated to be altered, the levels of ITGAV and ITGB3 were not statistically significantly changed, accompanied by the increased expression of some miRNAs (miR-23a, 146a-5p, miR-146b, and let-7g) and the decreased expression of some miRNAs (let-7d-5p and miR-184). As the pathway showed, the ITGAV and ITGB3 expression levels depended on the integrated results of multiple regulation. Therefore, the combined effect was hard to define when not only upregulated miRNAs were considered. As the accumulated effect of miRNAs on integrin in the present study remains elusive, further studies are needed to improve our understanding of the associations between these miRNAs and integrin. A limitation of this study is that it did not validate the effect of the miRNAs on aqueous humor outflow function, and this might be done in future studies. In conclusion, this study identified seven differentially expressed miRNAs in senescent AAP cells, which could contribute to IOP elevation by regulating three pathways responsible for cytoskeleton, focal adhesion, NO-cGMP signaling, and gap junction.
Appendix 1. Cell characterization of normal and senescent AAP cells.
Appendix 2. Microarray layout.
Appendix 3. qRT-PCR primer sequences.
Appendix 4. ITRAQ method in detail.
Appendix 5. Microarray raw data.
Appendix 6. Differentially expressed miRNAs.
Appendix 7. Target prediction.
Appendix 8. Bioinformatics annotations.
Appendix 9. Upregulated proteins by iTRAQ analysis.
Appendix 10. Downregulated proteins by iTRAQ analysis.
Appendix 11. Protein distribution.
This research was founded by State Key Program of National Natural Science Foundation of China (81430007), General Program of National Science Foundation China (81100662, 81371015, 81870661, 81470623), 211 Project of Fudan University (EHF158351), Young scientists program of EENT Hospital of Fudan University and BrightFocus Foundation (G2018112). Prof. Yuan Lei (email@example.com) and Prof. Junyi Chen (firstname.lastname@example.org) contributed equally to this research, and were regarded as co-corresponding authors.