Expression and predictive functional pro�les on MicroRNAs data in Vascularized Composite Allotransplantation Acute Rejection

Aim: To investigate the mechanisms of acute rejection for vascularized composite allotransplantation (VCA) using microRNAs (miRNAs) differential expression in a VCA animal model. Methods: Three brown Norway rats were used as transplant donors and three Lewis rats as VCA receptors. The changes were divided into different stages before and after transplantation in Lewis rats. The appearance changes of all subjects were recorded. Also, histological evaluations were performed on all recipients, and the microRNAs expression was analyzed when acute immune rejection occurred. Then, GO and KEGG Pathway analyses were employed to predict miRNA targets. Finally, differentially expressed miRNAs were detected by RT-qPCR. Results: Compared to pre-operation, 22 miRNAs were screened for signi�cantly different expressions. Among them, nine were upregulated and 13 were downregulated in skin tissues. The RT-qPCR results revealed that rno-miR-340-5p and rno-miR-21-5p were signi�cantly upregulated and enriched in the PI3K-Akt signaling pathway. Moreover, rno-miR-145-5p and rno-miR-195-5p were signi�cantly downregulated, and most of their target genes were enriched in the Hippo signaling pathway. The histological evaluations showed that, after VCA, the skin tissue presented severe acute rejection. Conclusion: Overall, the miRNAs rno-miR-340-5p, rno-miR-21-5p, rno-miR-145-5p, and rno-miR-195-5p were signi�cantly regulated during VCA acute rejection. These miRNAs are potential biomarkers for objective, early, and minimally invasive rejection diagnosis.


Introduction
For some end-stage diseases patients, organ transplantation can be a life-saving choice. In recent decades, the use of immunosuppressive drugs signi cantly improved the overall survival rate of solid organ transplant recipients [1,2]. However, compared with other organ transplantation types, vascularized composite allotransplantation (VCA) has a higher incidence of acute rejection [3]. Moreover, few studies regarding why the VCA acute rejection process occurs and how it is regulated are available.
MicroRNAs (miRNAs) are single-stranded noncoding RNAs with 19-25 nucleotides [4]. They are stably present in various tissues and blood. Also, miRNAs can bind to the 3' or 5' untranslated region of a target messenger RNA (mRNA), thereby inhibiting its function or leading to its degradation, nally regulating the occurrence and development of many diseases [5]. Recently, many miRNAs have been related to the activation and regulation of the innate immune system and in ammatory response [6,7]. Additionally, some reports have shown that the expression of miRNAs in plasma changes during acute rejection after organ transplantation, including heart or kidney [8,9]. Other studies have reported that plasma miRNA-146a and miRNA-155 can be potential acute rejection biomarkers after VCA [10]. Also, this research team continued to select several differentially expressed miRNAs in plasma for histological diagnosis [11].
However, regarding VCA acute rejection, there is no study with its whole transcriptome sequence and screening differentially expressed miRNAs. Therefore, in the present study, we aimed to evaluate miRNA changes during acute rejection on VCA model rats.

Animals
The VCA model was constructed using a rodent hindlimb allograft model. The Shanghai Sippr-BK Laboratory Animal Co. Ltd. (Shanghai, China) provided male Lewis and male Brown Norwegian (BN) rats (200-250 g) for this study. Two different receptors received one hindlimb from the same donor. All animal procedures were approved by the ethics committee of the Ninth People's Hospital A liated with the Shanghai Jiao Tong University School of Medicine. This work was performed following the guidelines of the Laboratory Animal Manual of the NIH Guideline for the Care and Use of Laboratory Animals.

Establishment of the hind limb transplantation model
Before operations, rats were anesthetized by continuous inhalation of iso urane. Operations were performed as previously described. Brie y, the two hind limbs of donor rats were dissected for the femoral vessels and related nerves, then amputated in the middle of the thigh. Before amputation, the limbs were perfused with cold heparinized Ringer's lactate solution, then stored in the same solution until transplantation. Each hind limb harvested was transplanted in situ to isolated recipient rats. The femur was xed with an 18G intramedullary needle. The dorsal muscle and skin were sutured with a 3 − 0 silk thread. Then the femoral artery, vein, and sciatic nerves were anastomosed with a 12 − 0 nylon thread under an operating microscope. After con rming that the anastomotic vessels were unobstructed, the ventral muscles and skin were sutured with a 3 − 0 silk thread (Fig. 1). Motor and sensory functions were not evaluated in the current study.
Three male BN rats were used as allograft donors and transplantations were performed into three male Lewis rats. The preoperative skin tissue of rats was used as controls. The postoperative skin tissue comprehended the experimental group and was used to analyze the preliminary mechanisms of acute immune rejection.

Tissue collection and preservation
Pre-and post-operation changes and special conditions observed in mice were daily recorded. After the 7th day, mice were sacri ced. Sterilized scissors were used to collect the skin of experimental mice. Skin tissue pieces were collected from each mouse. Removed skin tissue specimens were divided into two parts and stored. Half of the skin specimen was cut into pieces and frozen in liquid nitrogen for future use. The other half was placed in an EP tube and stored on dry ice for molecular detection. Each skin sample was marked with a corresponding label.
Then, sections were depara nized, rehydrated in an ethanol gradient, and stained with H&E trichrome according to standard protocols and using commercially available reagents.

RNA sequencing (RNA-seq)
The Hiseq Single-End sequencing was used in the present study. First splices were removed and cuts were performed according to the sequence quality. The original sequence was searched using a 5 base length as the window. When the average sequencing base quality in the window was lower than 20, the front end of the window was truncated and discarded. Filtered sequences were further used. Then, the number of clean reads (total reads) with a sequence length of 18-36 nt was counted. Identical sequences in each sample were duplicated and the sequence abundance was assessed. These sequences were called unique reads and were used for subsequent analyses.

Data collection and analyses
The miRNAs levels were evaluated using sequencing counts and normalized as counts per million of total aligned reads (CPM). Differentially expressed miRNAs were screened based on count values. Hierarchical clustering and Volcano plots were generated using R or Perl environments for statistical computing and graphics.

Target prediction and functional analyses
MiRNAs can bind to target sites mainly through complementary pairing. To analyze miRNA binding in animals, we predicted the target genes of differentially expressed miRNAs using miRanda and the 3 'UTR' sequence of mRNAs of the species as the target sequence. Top GO was used for GO enrichment analysis.
During this analysis, the list of miRNA target genes and the number of miRNA target genes of each term were calculated using differential miRNA target genes annotated by GO term. Then, the p-value was calculated by the hypergeometric distribution method (the standard for signi cant enrichment as p-value ≤ 0.05) to nd GO terms with signi cant enrichment of differential miRNA target genes compared to the whole genome background and to determine the main biological functions of these target genes. KEGG pathway enrichment analysis was also performed and the differentially expressed miRNA target genes top 20 pathways with the smallest p-values (most signi cant enriched) were selected.

Quantitative reverse transcription polymerase chain reaction (RT-qPCR)
First, RNA quality was assessed. Then, total RNA was reverse transcribed into cDNAs (PrimeScript TM 1st stand cDNA Synthesis Kit). Brie y, the reverse transcription reaction system 1 was added to the test tube in ice. Then, 10 µL of RNase-free dH 2 O was added. After mixing, samples were incubated at 65°C for 5 min, then quickly placed in an ice bath. Next, the reverse transcription reaction system 2 was added to the test tube in the ice. The reaction was performed at 42°C for 30-60 min. Samples were heated at 95°C for 5 min to end the reaction, then placed into ice for subsequent experiments or cryopreservation. For each target and housekeeping gene, the cDNA template of the sample was selected for PCR reactions (reaction A). The PCR reaction solution, con gured according to the reaction system A, was placed on the real-time PCR instrument to react. For uorescence signal measurements, the reaction was performed at 95°C for 5 min, then cycled 40 times at 95°C for 15 s and 60°C for 30 s.

Statistical analyses
Data are presented as means ± standard deviations (SDs). The miRNA expression followed a discrete distribution. Accordingly, signi cant differences between groups were compared using the negative binomial distribution. Differentially expressed miRNAs were screened based on count values. The signi cance level was set at p ≤ 0.05.

Skin visual examinations
The cumulative number of recipients with visual skin changes is shown in Fig. 2. The allogeneic transplantation hind limb led to increasingly worsening edema and swelling within a few days. Erythema and blisters appeared on the 7th day after transplantation. On the 10th day, large erythema and blisters areas appeared, and skin keratolysis was observed. On the 14th day, all transplanted hind limbs gradually became blackened and necrotized upward from the ngertip.

Histological evaluations
On the 7th day, early epidermis changes, necrosis, keratinization, epidermal thickening, lymphocyte in ltration, and basal cell vacuolation were observed in the upper dermis of post-transplant samples. Therefore, Grade I rejection occurred on the 7th day after transplantation. On the 10th day, the HE staining showed mixed small cell in ltration, increased basal cell vacuolation, and that the epidermis fell off. On the 14th day, the skin edema was severe, and mixed cell in ltration and cell necrosis reached Grade III rejection ( Fig. 3 and Table 1). The grading was based on the classi cation system of Büttemeyer et al. (1996).

Differentially expressed miRNAs
Page 6/17 The miRNA-seq was used to identify miRNA expression levels in control and experimental groups. Precise sample extraction, detection, and quality control processes were performed to control the sample quality through all steps. Then, the miRNA density distribution was used to investigate the expression pattern of all miRNAs in the sample. Overall, most miRNAs were averaged expressed, and low and high miRNAs expressions accounted for a small part of the total. The miRNA expression characteristics in each sample are shown in Fig. 4A.
We used the R heatmap software package for bidirectional cluster analysis of all miRNAs and samples. Clustering was performed according to the expression level of the same miRNA in different samples and the expression pattern of different miRNAs in the same sample. The Euclidean method was used to calculate the distance, and the hierarchical clustering longest distance method (complete linkage) was used for clustering (Fig. 4B).
According to hierarchical clustering results, differential genes were divided into different expression patterns. They were divided into different clusters (the gene expression trend in the same cluster is similar) to obtain gene clustering results.

Bioinformatic analyses
Moreover, miRNAs can bind to target sites mainly through complementary pairing. To analyze miRNAs bindings in rats, we used miRanda to target the 3 'UTR' sequence of their mRNAs. Target genes were predicted for differentially expressed miRNA sequences.
The GO enrichment analysis was performed using Top GO. During analysis, the list of miRNA target genes and the number of miRNA target genes of each term were calculated using the differential miRNA target genes annotated by GO terms. Then, the p-value was calculated using the hypergeometric distribution method (the standard for signi cant enrichment was p-value < 0.05) to nd the GO terms of different miRNA target genes with signi cant enrichment compared to the whole genome background and to determine the main biological functions performed by these miRNA target genes. Results were classi ed according to molecular functions (MF), biological processes (BP), and cellular components (CC). The rst 10 GO term items with the smallest p-values (most signi cant enriched) were selected for each GO classi cation (Fig. 8).
According to the GO analysis results, the enrichment degree was measured by the rich factor, FDR value, and the number of miRNA target genes enriched on the GO term. The rich factor refers to the ratio. The larger the rich factor, the greater the enrichment degree. The general FDR value range was 0-1. The closer it is to zero, the more signi cant the enrichment is. The rst 20 GO term entries with the lowest FDR values (most signi cant enriched) were selected for display ( Fig. 5A-B).
The enrichment analysis results provided direct acyclic diagrams of the three GO ontologies (cellular components, molecular functions, and biological processes). In this diagram, the closer to the root node, the more general the GO term is. Besides, the GO term branching down is the term annotated at a more precise level. By default, the program sets the top 10 GO terms with the highest signi cance to be squares, and the other terms to be circles. The darker the color, the more signi cant the GO term is. The colors from light to dark are: colorless, light yellow, dark yellow, and red ( Fig. 5C-H). Additionally, the KEGG pathway enrichment analysis was carried out according to the target gene results. The top 20 pathways with the smallest p-value (most signi cant enriched) were selected for display ( Fig. 5I-J).

Discussion
Compared with other organ transplantations, VCA has a higher acute rejection rate, but most patients can achieve long-term survival through timely treatment [12]. However, VCA mild rejection can damage many normal tissues, leading to dysfunction and even graft loss. Therefore, early diagnosis and timely treatment are very important for acute rejection. At present, rejection is mainly evaluated by subjective examination and blood indexes. These methods are not enough to directly re ect the current graft state and can not accurately predict the future development trend.
Moreover, miRNAs are non-coding RNAs that can be used as biomarkers of immune responses, in vivo index regulation, and tumor development [13,14]. Also, it has been reported that there is a difference in the expression of plasma miRNAs in acute rejection induced by VCA in rats [15], and the TaqMan miRNA reverse transcription kit has been used to detect the expression of known miRNAs in the skin and other tissues [11]. However, the miRNAs detected by these methods are very limited, and whether other miRNAs are important has not been fully clari ed. Therefore, we used miRNA-seq to detect miRNAs that are important during VCA acute immune rejection.
Also, rno-mir-340-5p and rno-mir-21-5p have many target genes enriched in the PI3K-Akt signaling pathway, which is closely related to corneal transplantation rejection. Additionally, rno-mir-340-5p plays a role in the progression of osteosarcoma through the PI3K Akt pathway [19]. However, the mechanism of rno-mir-340-5p/PI3K-Akt is rarely reported. Whether rno-mir-340-5p and its most abundant target gene regulate VCA acute immune rejection via the PI3K-Akt signaling pathway needs to be further explored.
Signi cantly downregulated rno-mir-145-5p and rno-mir-195-5p were enriched in the Hippo signaling pathway. However, little is known about the role of the Hippo signaling pathway in immune rejection. Loss of Hippo tumor suppressor activity and hyperactivation of Yap are commonly observed in cancers [20]. Tumor and liver-related Hippo's research has been very in-depth. In the future, we can focus on immunity mechanisms.
This study also has some limitations. First, we detected the expression of only four miRNAs. More miRNAs should be detected and the combination of different miRNAs should be analyzed to improve the speci city and sensitivity of rejection detection. Second, only rat miRNAs were used in this study and, although rats share many genes with humans, they are not the same. Hence, it would be better to have clinical specimens to detect human miRNAs to further assist clinical diagnosis and treatment.

Declarations
Funding This work was supported by a grant from National Major Science and Technology Projects of China (2018ZX10303502) and National Natural Science Foundation of China (82002065).

Competing Interests
The authors state that the study was conducted without any con ict of interest.

Author Contributions
YF and SW performed the research. YF and SW designed the research study. BS and JC contributed essential reagents or tools. XL, YX and YZ helped to analyze the data. YF wrote the article. SL and JZ revised the article. All authors read and approved the nal manuscript.

Data Availability
The datasets generated during or analysed during the current study are available from the corresponding author on reasonable request.

Ethical statement
No human studies were carried out by the authors for this article.    relations. The function range de ned from top to bottom is becoming smaller and smaller. The box represents the GO term with the enrichment degree of TOP10, and the darker the color, the higher the enrichment degree. (I-J) Histogram of KEGG pathway enrichment results in up-regulated (I) and downregulated (J) miRNA. The abscissa is the name, and the ordinate is -log10 (p-value) of each pathway.