RNA Sequencing Analysis of Changes in Pancreatic RNA in Tree Shrews with Traumatic Systemic Inflammatory Response Syndrome after Treatment with Umbilical Cord Mesenchymal Stem Cells

Systemic inflammatory response syndrome (SIRS) is a systemic nonspecific inflammatory response to severe injury caused by infectious or noninfectious factors, such as infection, trauma, burns, surgery and ischemia-reperfusion. SIRS presents as a group of clinical symptoms that ultimately lead to the body’s uncontrolled inflammatory response. Systemic reactions to severe infection include changes in body temperature, respiration, heart rate and white blood cell count (Anderson and Singh, 2017; Fang et al., 2016; Garcia-Lamberechts et al., 2018; Li et al., 2018; Mahassadi et al., 2018; Wagner et al., 2016; Ward et al., 2017; Xu, 2019).
       
A traumatic systemic inflammatory response syndrome model was established in 30 tree shrews by using the impact method to make unilateral femoral comminuted fractures and then intravenously injecting lipopolysaccharide. Another 10 normal tree shrews were taken as the control group and 10 were taken as the treatment after the model was established. The treatment group was infused with tree shrew umbilical cord mesenchymal stem cells via the tail vein. After 10 days of treatment, three tree shrews in each of the control group, model group and treatment group were killed and pancreatic tissue was collected for RNA sequencing analysis.

Construction of a tree shrew traumatic systemic inflammatory response syndrome model
 
Thirty tree shrews were given a unilateral femoral comminuted fracture by the weight impact method, Brieflyÿa hammer hit the femur head of a tree shrew under the force of gravity, causing a comminuted fracture and then 0.5 mg lipopolysaccharide was injected intravenously to create a traumatic systemic inflammatory response syndrome model. The body temperature of the tree shrews began to rise one day after the model was made and reached a maximum of 41.3°C 10 days after the model was made. The tree shrew white blood cell count began to increase 1 day after modeling and reached a maximum of 8.91 x 10⁹/L 4 days after modeling. The two markers used to confirm SIRS model establishment.
       
Another 10 tree shrews were used as the normal control group. For model identification, see published literature (Ruan et al., 2021).
       
The experiment was conducted at the Basic Medical Laboratory of the 920^th Hospital of the Joint Logistics Support Force. The experiment will take place from 2020 to 2021.
 
Treatment of the tree shrew traumatic systemic inflammatory response syndrome model
 
Of 30 tree shrew, the SIRS model was successfully constructed in 25 and 10 were selected as the treatment group. One hundred microliters of umbilical cord mesenchymal stem cell suspension (1 x 10⁶ cells) was reinfused into the tail vein. Umbilical cord mesenchymal stem cell suspension was provided by the National and Local Joint Engineering Laboratory of Biomedical Technology. All 10 tree shrews in the treatment group received cell transplantation.
       
The preparation and identification of umbilical cord mesenchymal stem cells were described in published literature (Ruan et al., 2016).
 
Collection of pancreatic tissue from three groups of tree shrews
 
After 10 days of treatment, 3 tree shrews were killed in each of the control group, the model group and the treatment group. The pancreas was removed washed with PBS, quickly placed into liquid nitrogen and sent to the company for RNA sequencing analysis. Three pancreases were collected from each of the three groups for RNA sequencing analysis.
 
RNA-sequencing analysis method to analyze the differences in transcripts in pancreatic tissue of the three groups of tree shrews
 
RNA extraction, library construction and sequencing
 
Total RNA was extracted using a TRIzol kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions and checked using RNase-free agarose gel electrophoresis. After extraction of total RNA, mRNA was enriched by Oligo(dT) magnetic beads. The library was constructed using the KAPA Stranded RNA-Seq Library Prep Kit (Illumina). The process was as follows: First, the enriched mRNA was fragmented into short fragments with RNA fragmentation buffer and the first-strand cDNA was reverse transcribed with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTPs and buffer. End repair was then performed. A bases were added and Illumina sequencing adapters were ligated. The library quality was assessed with an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) after PCR amplification and finally, the cDNA was sequenced using an Illumina Novaseq6000 (Gene Denovo Biotechnology Co.).
 
Bioinformatics analysis
 
Expression statistics
 
Based on the alignment results of HISAT2, StringTie was used to reconstruct transcripts and the expression levels of all genes in each sample were calculated. The expression level is displayed by the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) value. The formula for calculating the FPKM value is: In the formula, for the expression level FPKM value of gene A. C = Number of fragments aligned to gene A.
N = Number of all fragments aligned to the reference genome in the sample.
L = Number of bases in gene A.
 
Analysis of differences between groups
 
The input data of gene differential expression analysis were FPKM data obtained in gene expression level analysis, which was analyzed by ballgown software. Based on the difference analysis results, genes with P<0.05 and |log2FC|>0.585 were identified as significantly differentially expressed genes.
 
Difference comparison volcano chart
 
Volcano plot analysis was performed according to the significantly differentially expressed genes in each comparison group. The volcano plot visually displayed the differentially expressed genes between the comparison groups. In the graph, the closer the genes are to the two ends, the greater the difference.
 
Differential gene clustering heatmap
 
Hierarchical clustering of differential gene expression patterns was performed and a heatmap was used to present the clustering results. These genes with similar expression patterns may have common functions or participate in common metabolic pathways and signaling pathways. At the time of analysis, each gene was z score processed and then plotted.
 
KEGG (Kyoto encyclopedia of genes and genomes)       
Enrichment analysis
 
Pathway significant enrichment analysis uses KEGG Pathway as the unit and applies the hypergeometric test to find pathways that are significantly enriched in differentially expressed genes compared to background genes. The test assumes that N is the number of genes with pathway annotation in all background genes; n is the number of differentially expressed genes in N; M is the number of genes annotated as a specific pathway in all background genes; and m is the difference annotated as a specific pathway number of expressed genes. It is calculated as follows: After correction for multiple testing, pathways with a Q value < 0.05 were selected to be defined as pathways significantly enriched in differentially expressed genes. The Q-value here is the p value after FDR correction. The most important biochemical metabolic pathways and signal transduction pathways involving in differentially expressed genes were determined by pathway significant enrichment.
 
Statistical method
 
The statistical methods used in this experiment were calculation of hypothesis testing probability (p value) and multiple hypothesis testing correction to obtain FDR value (false discovery rate).

Among the genes upregulated in the model group compared to the control group and the genes downregulated in the treatment group compared to the model group, there were 11 genes in common through Venn diagram analysis (Fig 1A). These 11 common genes are shown in Table 1 and Table 2.






       
According to Table 1, a heatmap was drawn, as shown in Fig 1B. There are 3 genes that clustered well, namely, ENSTBEG00000013408, ENSTBEG00000006837 and ENSTBEG00000006937. These 3 genes were significantly elevated in the model group. Compared with the normal group, these three genes may be gene markers that are significantly elevated in the model group.
       
According to Table 2, a heatmap was drawn, as shown in Fig 2A. There are 5 genes that clustered well, namely, ENSTBEG00000009022, ENSTBEG00000013408, ENSTB EG00000017424, ENSTBEG00000006837 and ENSTBEG 00000003426. The expression levels of these 5 genes were significantly reduced after treatment. Compared with the model group, they may be the gene markers of reduction in the treatment group.


       
Among the genes downregulated in the model group and the control group and the genes upregulated compared with the treatment group and the model group, there were 11 shared genes identified through Venn diagram analysis (Fig 2B). These 11 shared genes are shown in Table 3 and Table 4.




       
According to Table 3, a heatmap was drawn, as shown in Fig 2C. There is one gene with better clustering/: ENSTBEG00000021330. This gene was significantly decreased in the model group. Compared with the normal group, this gene may be a gene marker that was significantly decreased in the model group.
       
According to Table 4, a heatmap was drawn as shown in Fig 2D. There is one gene with better clustering: ENSTBEG00000021330. This gene was significantly increased in the treatment group. Compared with the model group, it may be a meaningful gene marker for recovery after treatment.
 
KEGG results compared with the model group and the control group
 
3Compared with the control group, the KEGG pathway analysis results of downregulated differentially expressed genes (Fig 3A) were as follows:


       
Compared with the model group and the control group, the upregulated differentially expressed genes were subjected to KEGG pathway analysis (Fig 3B).
 
KEGG results of the comparison of the treatment and control groups
 
Compared with the control group, the KEGG pathway analysis results of downregulated differentially expressed genes (Fig 4A) were as follows.


       
KEGG pathway analysis of the upregulated differentially expressed genes compared with the control group (Fig 4B)
 
KEGG results in comparison of treatment and model groups
 
Compared with the treatment group and the model group, the KEGG pathway analysis results of down-regulated differentially expressed genes (Fig 5A) were as follows.


       
KEGG pathway analysis of the upregulated differentially expressed genes compared with the model group (Fig 5B).
 
Trend analysis results
 
Using the intersection of all differential genes in the normal group compared to the model group and the intersection of all the differentially expressed genes in the treatment group compared with the model group, a total of 35 genes were obtained. Trend analysis of these 35 genes was carried out and the results are as follows. The red box is a significant difference between the model group and the normal group and the treatment group. There were 11 genes in total. The blue box is the significant difference between the model group and the normal group and the treatment group, with a total of 7 genes. The lower figure (Fig 6) is the corresponding trend analysis result.


       
In modern local wars, all kinds of explosive weapons have high explosive power, many shrapels that can be projected in a fan shape or three-dimensional, large killing area and accurate targeting ability (Li et al., 2018). The killing effect of modern weapons has the characteristics of high speed, high efficiency, high intensity and soft kill (three high and one soft). These characteristics will cause serious injuries, mainly in the form of more serious injuries, multiple injuries and multiple burns. There are many psychological disorders and physiological imbalances, resulting in a high attrition rate, high shock rate and high operation rate (four more and three high), which complicates modern war injuries and brings greater difficulty to rescue (Hadziahmetovic, 1995; Klausner and Rozin, 1995). Among all kinds of war wounds, firearm injuries caused by high-speed and small-mass weapons are the highest proportion of war wounds in modern local wars. The wound is complicated and the infection is serious (Patzkowski et al., 2012; Yee et al., 2017). In research on the injury of important organs, it was found that after the maxillofacial injury of high-speed steel balls, the heart, lungs and other important organs of animals exhibited small-scale flaky bleeding spots. This is the pathological basis for the occurrence and development of serious complications, such as acute respiratory distress syndrome, disseminated intravascular coagulation and multiple organ failure (Del Sorbo and Slutsky, 2011; Huang et al., 2013; Lee et al., 2011; Lyons, 2010). This series of complications leads to posttraumatic syndrome and is an important reason for the high mortality rate of troops. Rapid and effective control and treatment of posttraumatic syndrome is a hot research topic at home and abroad. Therefore, it is of great military and scientific significance to establish a reproducible animal model for war-traumatic infection, systemic inflammatory response syndrome, shock and multiple organ failure and a new treatment method based on the animal model. The use of the tree shrew to establish a systemic inflammatory response syndrome model has the following advantages: the tree shrew is a nonhuman primate surrogate animal with abundant resources, low cost and a relatively close relationship with humans. In recent years, it has received increasing attention and the systemic inflammatory response syndrome model of tree shrews is rarely reported at home or abroad.
       
We established a tree shrew systemic inflammatory response syndrome model and treated it with umbilical cord mesenchymal stem cells. Previous studies have shown that the model was successfully established and that the treatment effect was obvious. In this paper, had control, treatment and model groups. RNA sequencing analysis of the specimens showed that some genes had statistically significantly altered expression. The two most significantly altered genes were ENSTBEG00000013408 (gene name TMEM211) and ENSTBEG00000006837 (gene name C1QL3), which had elevated expression in the model group. The decreased expression of these genes after treatment represents the change in inflammatory factors, indicating that the model group has a significant inflammatory response. Furthermore, the treatment of umbilical cord mesenchymal stem cells has the effect of reducing inflammatory factors, which further proves the anti-inflammatory response of umbilical cord mesenchymal stem cells. Another gene with obvious changes in expression was ENSTBEG00000021330 (gene name ATP8), which was downregulated in the model group (Fold change=0.360939) and upregulated after treatment (Fold change=4.476301), indicating that ATP8 is a gene that contributes to the recovery of systemic inflammatory response syndrome. In the future, we may use this gene to treat systemic inflammatory response syndrome.

These results all show that tree shrew umbilical cord mesenchymal stem cells are effective in the treatment of systemic inflammatory response syndrome. The experimental results of tree shrews have laid a solid foundation for the treatment of human beings. In the future, in-depth research is needed before further extension to modern military medicine.

We thank American Journal Experts for assisting with the English preparation of this manuscript.

All experimental protocols were approved by the Experimental Animal Ethics Committee of the 920^th Hospital of Joint Logistics Support Force of PLA (2020-022-01). All methods were performed in accordance with the relevant guidelines and regulations.

This work was supported by grants from the Yunnan Science and Technology Plan Project Major Science and Technology Project (2018ZF007), the Yunnan Province Applied Basic Research Program Key Project (202101AS070039), the 920^th Hospital of the PLA Joint Logistics Support Force In-hospital Technology Plan (2019YGB17,2020YGD12).

Guang-Ping Ruan and Xiang Yao made substantial contributions to the study conception and design, data acquisition and data analysis and interpretation.
Guang-Ping Ruan, Mei-yu Lin and Xiang Yao conducted the experiments.
       
Xiang Yao and Shu-Qian Lin agree to be accountable for all aspects of the work and ensure that questions related to the accuracy or integrity of any part of the study will be appropriately investigated and resolved.
       
Xing-Hua Pan and Guang-Ping Ruan provided final approval of this version of the manuscript for publication. Guang-Ping Ruan, Xing-Hua Pan and Shu-Qian Lin were involved in drafting the manuscript or revising it critically for important intellectual content.
       
All authors read and approved the final manuscript.

All data generated or analysed during this study are included in this article. Further enquiries can be directed to the corresponding author.

The authors declare that they have no competing interests.