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. The ethics committee of the Ninth People’s Hospital, affiliated with the Shanghai Jiao Tong University School of Medicine, approved all animal procedures. This work followed the guidelines of the Laboratory Animal Manual of the NIH Guideline for the Care and Use of Laboratory Animals.

Before operations, rats were anesthetized by continuous inhalation of isoflurane, and operations were performed as previously described [12]. Briefly, 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 fixed with an 18G intramedullary needle. and 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 confirming 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.

Fig. 1.

Surgical procedures. (A) Prepare the donor site for amputation. (B) Amputation. (C) Prepare the receptor. (D) Perform vascular anastomosis. (E, F) Close muscles and skin.

Male BN rats were used as allograft donors and transplantations were performed on male Lewis rats. The preoperative skin tissue of rats was used as the control. The postoperative skin tissue comprehended the experimental group and was used to analyze the primary mechanisms of acute immune rejection.

Pre- and post-operation changes and special conditions observed in rats were recorded daily. Rats were sacrificed after the 7th day. Sterilized scissors were used to collect the skin of each experimental rat. Skin tissue pieces were collected from each rats. 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 further 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.

First, fixed skin tissue samples were embedded in paraffin, sectioned (4–5 $\mu$m), and mounted on slides. Then, sections were deparaffinized, rehydrated in an ethanol gradient, and stained with H&E trichrome according to standard protocols and using commercially available reagents.

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.

The miRNA 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.

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 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 significant enrichment as p-value $\le$ 0.05) to find GO terms with significant 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 significantly enriched) were selected.

First, we assessed the RNA quality, and then total RNA was reverse transcribed into cDNAs (PrimeScript TM 1st stand cDNA Synthesis Kit). Briefly, the reverse transcription reaction system 1 was added to the test tube in ice. Then, 10 $\mu$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, configured according to the reaction system A, was placed on the real-time PCR instrument to react. For fluorescence 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.

Data are presented as means $\pm$ standard deviations (SDs). The miRNA expression followed a discrete distribution. Accordingly, significant differences between groups were compared using the negative binomial distribution. Differentially expressed miRNAs were screened based on count values. The significance level was set at p $\le$ 0.05.

The cumulative number of recipients with visual skin changes is shown in Fig. 2. The allogeneic hind limb transplantation 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 blister areas appeared, and skin keratolysis was observed. On the 14th day, all transplanted hind limbs gradually became blackened and necrotized upward from the fingertip.

Fig. 2.

Allograft images at post-operative days 3 (A), day 7 (B), day 10 (C) and day 14 (D).

On the 7th day, early epidermis changes, necrosis, keratinization, epidermal thickening, lymphocyte infiltration, 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 slight cell infiltration, increased basal cell vacuolation, and that the epidermis fell off. On the 14th day, the skin edema was severe, and mixed cell infiltration and cell necrosis reached Grade III rejection (Fig. 3 and Table 1). The grading was based on the classification system of Büttemeyer et al. [13].

Fig. 3.

Hematoxylin and eosin staining of skin samples. Allograft skin tissues at days 0 (A), day 7 (B), day 10 (C) and day 14 (D). Bar; 250 $\mu$m.

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 their expression pattern 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 presented in Fig. 4A.

Fig. 4.

Screening differentially expressed miRNA. (A) The horizontal line in the middle of the box is the median. The upper and lower edges of the box are 25% and 75% respectively, and the upper and lower limits are 10% and 90% respectively. The external shape is kernel density estimation. (B) Differential expression miRNA clustering. Horizontal represents miRNA, and each column is a sample. Red represents high expression miRNA and green represents low expression miRNA. (C, D) Trend analysis. The gray line in the figure shows the expression pattern of genes in each cluster, and the blue line represents the average expression of all genes in the cluster in the sample. (E) Volcano diagram of differentially expressed genes. The abscissa is –log10 (p-value) and the ordinate is –log10 (p-value). The two vertical dashed lines in the figure are twice the expression difference threshold; The horizontal dotted line is p-value = 0.05 threshold. The red dot indicates the up-regulated genes, the blue dot indicates the down-regulated genes, and the gray dot indicates the nonsignificant differentially expressed genes.

We used the “heatmap” R 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).

The hierarchical clustering divided differential genes into various expression patterns. These genes were divided into different clusters (the gene expression trend in the same cluster is similar) to obtain gene clustering results. The clustering of upregulated and downregulated miRNAs is shown in Fig. 4C,D, respectively. Moreover, the Volcano plots provided a quick visual identification of the miRNAs displaying large and statistically significant magnitude changes (Fig. 4E). The tRF & tiRNA selection criteria included a higher FC, lower q-value, and higher CMP. After summarizing the original data, 22 differentially expressed miRNAs were selected (9 significantly upregulated and 13 significantly downregulated) (Table 2).

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 and the number of miRNA target genes were calculated for each term using the differential miRNA target genes annotated by GO. Then, the p-value was calculated using the hypergeometric distribution method (the standard for significant enrichment was p-value 0.05) to find the GO terms of different miRNA target genes with significant enrichment compared to the whole genome background and to determine the main biological functions performed by these miRNA target genes. Results were classified according to molecular functions (MF), biological processes (BP), and cellular components (CC). The first 10 GO terms with the smallest p-values (most significantly enriched) were selected for each GO classification (Fig. 5).

Fig. 5.

Bioinformatic prediction. Bubble Diagram of GO enrichment analysis in up-regulated (A) and down-regulated (B) miRNA. (C–H) Up-regulated gene ontology DAG enrichment map of biological process (C), molecular function (D) and cell component (E). Down-regulated gene ontology DAG enrichment map of biological process (F), molecular function (G) and cell component (H). Each node represents a GO term, and branches represent inclusion relations. The function range defined 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 down-regulated (J) miRNA. The abscissa is the name, and the ordinate is –log10 (p-value) of each pathway.

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 in 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 significant the enrichment is. The first 20 GO term entries with the lowest FDR values (most significantly enriched) were selected for display (Fig. 5A,B).

The enrichment analysis results provided direct acyclic diagrams of the three GO ontologies (CC, MF, and BP). In this diagram, the closer the root node, the more general the GO term. 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 significance as squares, and the other terms as circles; The darker the color, the more significant the GO term. 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 significantly enriched) were selected for display (Fig. 5I,J).

GO and KEGG analysis were performed on 9 up-regulated and 13 down-regulated miRNAs. The ontology covers three domains: Biological Process, Cellular Component, and Molecular Function. Up-regulated miRNAs were mainly concentrated in Cell Periphery and Plasma Membrane. Fig. 5B indicates that down-regulated miRNAs were mainly concentrated in Cell Periphery and play a role in Cellular Response to Stimulus, Multicellular Organismal Process and Signaling (Fig. 5B).

Pathway analysis is a functional analysis whereby genes are mapped to KEGG pathways. Up-regulated miRNAs were mainly enriched in PI3-AKT signaling pathway (Fig. 5I). The Hippo pathway had the highest enrichment score in down-regulated miRNA (Fig. 4D). Therefore, we focused on this path. Among them, Hippo pathway not only had the most enriched genes, but also had a high degree of enrichment; accordingly, it became the focus of our follow-up experiments.

Finally, we verified the expressions of rno-mir-340-5p, rno-mir-21-5p, rno-mir-145-5p and rno-mir-195-5p by RT-qPCR (Fig. 6). The designed primer sequences are shown in Table 3 and their relative expression levels were analyzed by t-tests with a p $\le$ 0.05 as statistically significant. The rno-mir-340-5p and rno-mir-21-5p were significantly upregulated compared to pre-operative samples, while rno-mir-145-5p and rno-mir-195-5p were significantly downregulated. Then, we predicted their target genes and signaling pathways as described above.

Fig. 6.

Validation of the four selected miRNAs using RT-qPCR verification. Compared with the day 0th, (A,B) rno-miR-340-5p and rno-miR-21-5p were up-regulated; (C,D) rno-miR-145-5p and rno-miR-195-5p were down-regulated. The data were normalized using the mean $\pm$ standard error of the mean (SEM). ***p $\le$ 0.001, ****p $\le$ 0.0001.