MIR31 May Mediate Multiple Biology Processes To Promote Vocal Fold Wound Healing

The aim of this study was to explore the role of MIR31 in vocal fold wound healing (VFWH) and its possible mechanism. MIR31 expression was studied in rat vocal fold tissue at multiple time points (1, 4, and 8 weeks) after vocal fold injury. Co-expression analysis, pathway analysis, and literature-based network analysis were conducted to explore the possible mechanism underlying the relationship between MIR31 and VFWH. MIR31 expression was signicantly elevated after vocal fold injury (p<5.65e-5), which also presented decreased expression in the late stage of VFWH process compared to the early and middle stages (p<5.40e-3). MIR31 also presented strong co-expression with 17 VFWH-signicant genes (absolute value of ROH ∈ (0.63, 0.83)), which were mainly involved in collagen production. Literature-based network analysis showed that MIR31 could suppress two inhibitors (gene SMAD1 and HDAC2) of wound healing and activate one promoter (adenosine triphosphate). MIR31 could also mediate multiple biology processes that were associated with wound healing, including keratinocyte proliferation, collagen production, and inammation. This study supports a strong association between MIR31 and the process of vocal fold wound healing, which may be related to collagen synthesis and other biological processes that need further study.


Introduction
The structure of the extracellular matrix (ECM) is essential for the functionality of the vocal fold to produce sounds (Gray et al., 2000). Vocal folds injury could lead to abnormal massive hyperplasia of the brous tissue and ECM disorder in the lamina propria. Vocal folds injury caused brous scars could bring persistent pararthria and irreversible changes (Choi et al., 2012;Hu et al., 2014;Tateya et al., 2006). So far, it is still a challenging clinical topic to treat vocal fold scars effectively.
Vocal fold healing could take from two weeks to several months (Tateya et al., 2006;Ling et al., 2010;Gurtner et al., 2008). Multiple genes and proteins have been shown to play roles in the wound healing process (Luo et al., 2020;Spallotta et al., 2013). For instance, SMAD1 has been suggested as a therapeutic molecular target in skin wound healing that plays an active role in wound repair and regenerative medicine (Luo et al., 2020). Spallotta et al.'s study showed that HDAC2 inhibition could signi cantly enhance skin wound repair (Spallotta et al., 2013).
MIR31 is a short non-coding RNA that is involved in post-transcriptional regulation of gene expression in multicellular organisms by affecting both the stability and translation of mRNAs (Wang et al., 2014;Li et al., 2015). Multiple previous studies suggested that MIR31 plays a role in the pathology of skin wound healing, including the stimulation of wound contraction to enhance the wound closure (Chen et al., 2019) and the promotion of skin wound healing by enhancing keratinocyte proliferation and migration (Gerlach and Vaidya, 2017). However, so far, it is not clear whether MIR31 contributes to the vocal fold wound healing (VFWH).
In this study, a microarray was used to detect the expression changes of MIR31 in different stages of vocal fold wound healing. In addition, a co-expression analysis was conducted between MIR31 and the genes that presented signi cant expression changes in the vocal fold after injury. Our results indicated that the expression of MIR31 was stimulated by the wound and was signi cantly related to the vocal fold's wound healing process. This study guarantees further study to explore the role of MIR31 and VFWH.

mRNA and miRNA data extraction
The animal experiments in this study were performed in accordance with the regulations of Peking University Animal Care and Use Committee (LSC-ZhangY-1) and license permit: 1116012800123. Male Sprague-Dawley (SD) rats aged 14 to 16 weeks and weighing 400 to 450 g underwent unilateral or bilateral vocal fold injury using procedures described in an earlier study (Tateya et al., 2005). The vocal fold was injured by separating and removing the lamina propria from the thyroarytenoid muscle. First, 20 animals were randomly divided into four experimental groups ( ve animals in each group) based on time of sacri ce: uninjured control, and 1, 4, and 8 weeks after injury. Following sacri ce, the larynx was harvested, and the bilateral vocal folds were dissected under magni cation. Each specimen was immediately snap-frozen in liquid nitrogen and stored at − 80°C for subsequent study.
Total RNA, including both mRNA and miRNA, was extracted from vocal folds using TRIzol Reagent (Life Technologies, Carlsbad, CA, USA) and puri ed with an RNeasy mini kit (Qiagen, Valencia, CA, USA).
Biotinylated cDNA was prepared according to the standard Affymetrix protocol from 250 ng of total RNA using an Ambion® WT Expression Kit. Following labeling, fragmented cDNA was hybridized for 16 h at 45°C with the Clariom™ S Assay (rat, Affymetrix). GeneChips were washed and stained in the Affymetrix Fluidics Station 450. All arrays were scanned using an Affymetrix® GeneChip Command Console (AGCC), which was installed in a GeneChip® Scanner 3000 7G.
The row data (.cel) were normalized using TAC software (Transcriptome Analysis Console, Version 4.0.1) with the Robust Multichip Analysis (RMA) algorithm using Affymetrix default analysis settings and global scaling as normalization method. Values presented are log2 RMA signal intensity. The microarray data discussed in this article are submitted to NCBI Gene Expression Omnibus (GEO) with accession number GSE139383.

Screening of differentially expressed genes
For the mRNA and miRNA expression data, we rst used the Limma R package (version 3.36.5) to lter the differentially expressed genes (DEGs) for the injured groups of different weeks by comparing them with the uninjured control group. The Benjamini-Hochberg procedure was used for the multiple test false discovery rate (FDR) control (q = 0.05). The threshold set for up-and downregulated genes was fold change > 1.5 or <-1.5, FDR corrected P-value < 0.01. Then, we used one-way ANOVA to identify among these DEGs that also showed signi cant expression variance in the three injured groups (p < 0.01). After that, we calculated the co-expression between the signi cant miRNA(s) and the mRNA(s) with the purpose to explore the potential association between these miRNA(s) and mRNA(s). mRNA and miRNA that satisfy the following criteria were selected for further analysis: 1) when compared to the control group, it demonstrates signi cance in each injured groups (fold change > 1.5 or <-1.5, and p < 0.01); 2) when compared among the three injured groups, it also show signi cance (p < 0.01); 3) mRNAs showed signi cant correlation with the miRNA expression (ROH of Pearson Correlation coe cients > 0.6).

Functional network and pathway analysis
To explore the functional connection between the identi ed miRNA(s) and the mRNA, a literature-based functional network analysis has been conducted, which uncovers miRNA-gene relationships supported by previous studies. In addition, a gene set enrichment analysis (GSEA) has been performed against Gene Ontology (GO) terms and Human Protein Atlas Expression Ontology (HPAEO), which reveals the functional pro le of the mRNA corresponding genes. Pathways/GO items show signi cance were reported (FDR corrected p < 0.05 & overlap ≥ 2).

MIR31 expression
When compared to uninjured vocal fold tissues (control group, n = 5), MIR31 demonstrated signi cantly increased expression in injured vocal fold tissues at different time points (1, 4, and 8 weeks; n = 5 in each experimental group) with p < 5.65e-5, as shown in Fig. 1a. It was found that expression of MIR31 also varied signi cantly among injured vocal fold tissues at different time points (p < 0.0054; see Fig. 1b). The details of the expression pro le of MIR31 in the four groups were presented in supplementary data MIR31_WFWH◊Expression pro le.

Expression of 16 signi cant genes
Our analysis also identi ed 18 genes that present signi cance when compared among both four groups (p < 0.0099) and three injured groups (p < 0.0075). The expression heat map of these genes and MIR31 was presented in Fig. 2a.
Interestingly, 16 out of the 18 genes presented a strong positive correlation with MIR31 in terms of expression variation among the four groups (ROH (0.63, 0.83)), and two showed a robust negative correlation with MIR31 (ROH=-0.78 and − 0.70, respectively). These results suggested that MIR31 may ∈ regulate or co-function with these genes to in uence the wound healing process of the vocal fold. The details of the analysis results were presented in supplementary data MIR31_WFWH◊Expression pro le.

GSEA results
GSEA was performed to compare the 17 genes and MIR31 with the GO terms and HPAEO terms. Results showed that these 17 genes were mainly involved in collagen related GO terms, as shown in Fig. 3. No HPAEO gene group was identi ed. The details of these GO terms were presented in supplementary data MIR31_WFWH◊GSEA. Figure 3 presents the biological network connecting MIR31 and wound healing. The network were built based on previous study results, which was identi ed through literature data mining assisted by Pathway Studio (www.pathwaystudio.com). Each relationship (edge) within the pathways were supported by one or more references. In total, the two networks in Fig. 3 were built based on the results of 665 studies. The details of these references were presented in supplementary material MIR31_VFWH◊ Ref4_MIR31_Network, including titles and the sentences where a relationship has been identi ed.

Functional network connecting MIR31 and wound healing
The networks presented in Fig. 3 showed that, MIR31 could drive multiple molecules and cell processes that play roles in pathological wound healing. According to the polarity of the network, most of the relationships supporting a positive role of MIR31 in wound healing. That is, elevated expression of MIR31 could promote the wound healing. However, we also noted that MIR31 could inhibit the production of two wound healing promoters, including nitric oxide (NO) and progesterone. This indicates the biological complexity of the wound healing process and the mixed in uence of MIR31 on wound healing.

Discussion
Multiple studies suggested that MIR31 plays a role in skin wound healing through the promotion of keratinocyte proliferation and migration (Chen et al., 2019; Gerlach and Vaidya, 2017). Results from this study showed that MIR31 was stimulated by vocal fold wound to demonstrate signi cantly high expression at the initial stage of vocal fold healing, which decreased in the late stage of the wound healing process. Moreover, MIR31 demonstrated a strong correlation with 17 genes showing signi cant expression variation during VFWH, which were mainly involved in collagen synthesis. Our results suggest that MIR31 is closely related to VFWH.
Expression analysis showed that MIR31 expression was signi cantly increased at the initial stage of VFWH (week 1), as shown in Fig. 1a. Upregulated MIR31 has been reported to alleviate in ammation in colon injury (Orlova et al., 2019) and cardiac injury (Miguel et al., 2018). At the early stage of wound healing, in ammation is the initial response to cellular injuries and is the key process in wound healing (Meng X.M., 2019). MIR31 has also been shown to stimulate wound contraction and thus enhance the wound closure (Chen et al., 2019). Therefore, our results indicate that increased MIR31 expression could help at the early stage of the VFWH.
The up-regulated expression of MIR31 lasted to week 3 without obvious decrease, as shown in Fig. 1b. MIR31 has been shown to promote skin wound healing by enhancing keratinocyte proliferation and migration, which may happen through suppressing its direct target gene, epithelial membrane protein 1 (EMP-1), during wound healing (Gerlach and Vaidya, 2017). MIR31 could also suppress the inhibitors of wound healing, including Mothers Against Decapentaplegic Homolog 1 (SMAD1) (Costa et al., 2019) and histone deacetylase 2 (HDAC2) (Wei et al., 2017). Moreover, up-regulation of MIR31 could lead to elevated cellular adenosine triphosphate (ATP) that are required for wound closure (Zhu et al., 2019;Handly and Wollman, 2017). Therefore, lasted overexpression of MIR31 could exert continued aid to the VFWH.
We also noticed that MIR31 expression was signi cantly dropped 8-weeks after vocal fold wound (p < 0.0054), which may be corresponding to the late stage of the wound healing process (Fig. 1b). However, this may also re ect balanced factors that in uenced VFWH. As shown in Fig. 4 Co-expression analysis showed that MIR31 was strongly related to 17 genes that presented signi cant expression changes during VFWH (Fig. 2). These genes were mainly involved in collagen biosynthetic and trimming, as shown in the GSEA results presented in Fig. 3.
Collagen has been shown to increase keratinocyte proliferation, positively acting on cell entry in the mitotic phase (Luparello et al., 2020). Thus, collagen synthesis and deposition into the wound is essential during wound healing (Han et al., 2012). Our results suggested that MIR31 may co-function with the 17 collagen production regulators to play a role in VFWH process. In addition, our study also showed that MIR31 was related to multiple biology processes that were linked to wound healing, including brosis and keloid, as shown in Fig. 4b.

Conclusion
This study supports a promotion role of MIR31 in the process of vocal fold wound healing, which may be through the regulation of multiple biological processes, including keratinocyte proliferation, collagen production, ATP generation, and in ammation.

Data Availability Statement
All data generated during this study are included in this article and its additional information les, and all the materials generated during this study are available upon request.

Author Contributions
HW developed the method, performed the experiments, and wrote the manuscript. All authors read and approved the manuscript.

Funding
This work was supported by programs of the National Natural Science Foundation of China (81870709).