The Network of miRNA–mRNA Interactions in Circulating T Cells of Patients Following Major Trauma – A Pilot Study

Abstract

Purpose

Following major trauma, genes involved in adaptive immunity are downregulated, which accompanies the upregulation of genes involved in systemic inflammatory responses. This study investigated microRNA (miRNA)-mRNA interactome dysregulation in circulating T cells of patients with major trauma.

Patients and Methods

This study included adult trauma patients who had an injury severity score ≥16 and required ventilator support for more than 48 h in the intensive care unit. Next-generation sequencing was used to profile the miRNAs and mRNAs expressed in CD3+ T cells isolated from patient blood samples collected during the injury and recovery stages.

Results

In the 26 studied patients, 9 miRNAs (hsa-miR-16-2-3p, hsa-miR-16-5p, hsa-miR-185-5p, hsa-miR-192-5p, hsa-miR-197-3p, hsa-miR-23a-3p, hsa-miR-26b-5p, hsa-miR-223-3p, and hsa-miR-485-5p) were significantly upregulated, while 58 mRNAs were significantly downregulated in T cells following major trauma. A network consisting of 8 miRNAs and 22 mRNAs interactions was revealed by miRWalk, with three miRNAs (hsa-miR-185-5p, hsa-miR-197-3p, and hsa-miR-485-5p) acting as hub genes that regulate the network. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis suggested that “chemokine signaling pathway” was the predominant pathway.

Conclusion

The study revealed a miRNA-mRNA interactome consisting of 8 miRNAs and 22 mRNAs that are predominantly involved in chemokine signaling in circulating T cells of patients following major trauma.

Keywords:

• miRNA
• next-generation sequencing
• T cells
• adaptive immunity
• major trauma
• critical illness

Introduction

Major trauma can lead to a systemic inflammatory response and disrupted immune system homeostasis.^1,2 The cytokine/genomic “storm” that occurs after major trauma is well described^3,4 and involves the simultaneous upregulation of genes involved in systemic inflammatory responses and the downregulation of genes involved in adaptive immunity.³ These immunoinflammatory responses greatly impact the prognosis of patients with major trauma.⁵ Prolonged immune suppression may lead to increased susceptibility to secondary infections.⁷ A retrospective analysis of 917 patients with an injury severity score (ISS) ≥16 revealed that the absolute lymphocyte counts on day 3 and the ratio of absolute lymphocyte counts on day 3 and day 1 after trauma are independent risk factors for sepsis and patient mortality.⁸ In addition, the extent of immune dysfunction following major trauma is strongly correlated to patient morbidity and mortality.⁶

Immunosuppression after major trauma is characterized primarily using T cell populations of the adaptive immune system.^7,9 Weighted gene co-expression network analysis of trauma patients revealed that modules with high T cell activation and low neutrophil activation indicate better survival of patients with sepsis.¹⁰ Genome-wide expression analysis of whole-blood leukocytes in 167 adult patients with severe blunt trauma revealed widespread transcriptome changes, with alterations in up to 80% of the leukocyte transcriptome.³ The genomic response to traumatic injury involves the upregulation of several genes involved in the mediation of inflammation, pattern recognition, and antimicrobial activity, along with downregulation of genes engaged in antigen presentation, natural killer (NK) cell function, and T cell proliferation and apoptosis.³ Of the interleukin (IL)-10, IL-6, and p38 mitogen-activated protein kinase (MAPK) pathways are the most upregulated pathways associated with complicated recovery, whereas T cell regulation and antigen presentation are among the most downregulated pathways.³

MicroRNAs (miRNAs) are endogenous small non-coding RNAs that regulate gene expression through post-transcriptional modulation,¹¹ primarily by binding to the 3′ UTR of target mRNAs, which leads to mRNA degradation or translation inhibition.¹² miRNAs are important regulators of almost all cellular processes,¹³ including the immune system.^14–16 The loss of certain miRNAs or genetic ablation of miRNA machinery can severely compromise immune system development and immune responses, thus leading to immune dysfunction.¹⁴ A high-throughput study identified 69 dysregulated circulating miRNAs in patients with major trauma 12 h after intensive care unit (ICU) admission compared to healthy controls.¹⁷ Among these dysregulated miRNAs, 14 were correlated with innate immune responses involving the toll-like receptor (TLR) 3, TLR4, myeloid differentiation primary response 88 (MYD88), and Toll-receptor-associated molecule (TRAM) pathways. However, little is known about the role of miRNAs in major trauma, especially in T cells of the adaptive immune system. Therefore, we investigated dysregulated miRNAs and potential mRNA targets in circulating T cells of major trauma patients. To reduce patient-to-patient variations, we used samples from the same patients during the recovery stage as control samples.

Patients and Methods

Patient Enrollment

Trauma patients were recruited as study participants after admission to the trauma ICU of a level I trauma center in Southern Taiwan^18–20 between December 2017 and December 2018. Only patients who fit the following criteria were included: (1) age of 20 years or above; (2) admission to the ICU due to trauma injury; (3) ISS ≥ 16;^21–23 and (4) the use of ventilator support for more than 48 h. Patients who fit the following criteria were excluded: (1) immunocompromised condition; (2) malignancy; or (3) unwillingness to participate in the study. Forty trauma patients with critical illness who stayed in the ICU were enrolled in the study (Figure 1). Written informed consent was collected from each participant. Patient medical information, including sex, age, Glasgow Coma Scale (GCS), Abbreviated Injury Scale (AIS) in each body region, and ISS, was collected. The AIS is the basis of many severity scoring systems that assess the severity of an anatomical injury on a six-point ordinal scale ranging from 1 to 6 as minor, moderate, serious, severe, critical, and un-survivable, respectively.

Figure 1 Enrollment of the patients and flowchart illustrating the collection and processing of samples for the experiments.

Abbreviations: ICU, intensive care unit; ISS, injury severity score; QC, quality control; NGS, next-generation sequencing; miRNAs, microRNAs; RT-qPCR, real-time quantitative reverse transcription polymerase chain reaction; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Specimen Collection

Peripheral blood samples were collected from each patient within 72 h after admission to the ICU. Fasting blood samples (10 mL) were collected from the forearm at this time point and were defined injury samples. After the patients successfully left the ICU, another peripheral blood sample was collected immediately before discharge from patients who stayed in the hospital. The blood samples drawn at this time point were defined recovery samples. If the patients did not survive, collection of the second blood sample could not be performed. This study only compared miRNA and mRNA expression in injury and recovery samples from the same patients. Vacutainer Serum Separator Tubes (BD Diagnostics, Franklin Lakes, NJ, USA) were used to obtain serum from the blood samples without using additives. The collected blood samples were preserved in ice and were immediately sent to the laboratory. Ficoll–Paque Premium (17–5442-02, Merck, Kenilworth, NJ, USA) was used to separate the peripheral blood mononuclear cells (PBMCs) from whole blood by the density gradient centrifugation. Four-milliliter Ficoll solution was added to the bottom of a 15 mL tube before adding 8 mL whole blood. PBMCs were collected from the Ficoll/plasma interface after centrifugation at 400 g for 40 min at 25 ℃. The samples were diluted by mixing equal amounts of phosphate-buffered saline. After centrifuging the diluted samples at 240 g for 5 min, the supernatant was removed and the cell pellets were resuspended in 1X BD IMag buffer (BD Diagnostics). Cell numbers were counted. The cell suspensions were centrifuged at 200 g for 10 min. For every 10⁷ cells in the supernatant, 50 μL of BD IMag anti-human CD3 magnetic particles (552,593, BD Bioscience) was added. The mixtures were then incubated at 37 ℃ for 30 min. After washing with an equal amount of BD IMag buffer, CD3⁺ cell pellets were collected and resuspended in 500 μL QIAzol lysis reagent from an RNeasy Mini Kit (74,104, Qiagen, Venlo, Netherlands). RNeasy Mini Kits (Qiagen) was used to extract total RNA, which was subsequently quantified using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA) and a Qubit RNA Assay Kit (Q10210, Thermo Scientific). A Caliper LabChip Analyzer (PerkinElmer, Waltham, MA, USA) was used to determine the RNA quality based on the RNA integrity number (RIN).

Next-Generation Sequencing (NGS) Analysis of miRNAs

Total RNA from the injury samples (n = 10) and the corresponding recovery samples (n = 10) were pooled for every five patients to create two pooled injury samples (n = 2) and two pooled recovery samples (n = 2). GeneTech Biotech Co., Ltd (GeneTech, Taipei, Taiwan) performed the small RNA cloning and NGS analysis. Briefly, miRNA 15–30 nucleotides long were passively eluted from polyacrylamide gel, precipitated with ethanol and melted in double-distilled water. An Illumina TruSeq Small RNA Sample Prep Kit was used to ligate small RNA linkers onto prepared bar-coded cDNAs. 3′ and 5′ adapter-ligated RNA in 1 μg total RNA was reverse-transcribed with Invitrogen SuperScript II Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA) and was then amplified by 15 cycles of polymerase chain reaction (PCR). To reduce sample bias, 15 barcoded reverse primers were ligated directly to the miRNAs. A BioAnalyzer 2100 (Agilent Technologies) was used to analyze individual libraries to evaluate the presence of linked cDNA and 15 bar-coded libraries (135–165 bp).

RT-qPCR for miRNA Expression

The expression of significantly upregulated miRNAs identified by NGS analysis of the injury and recovery samples (n = 28 each) was determined by real-time quantitative reverse transcription polymerase chain reaction (RT-qPCR). TaqMan MicroRNA Reverse Transcription Kits (Applied Biosystems, Foster City, CA, USA) were used to reverse transcribe the RNA into cDNA. Then, TaqMan Universal PCR Master Mix (No UNG, PN 4324018, Applied Biosystems) and specific miRNA primers from TaqMan MicroRNA Assays (Applied Biosystems) were used to amplify the cDNA. About 25 fmol single-stranded cel-miR-39 (Invitrogen) was added to each sample for use as an internal control to determine miRNA expression. RT-qPCR was performed using a 7500 RT-qPCR system (Applied Biosystems) and SYBR Green (Applied Biosystems). miRNA expression levels were quantified using the 2^−ΔΔCt method using normalized cycle threshold (Ct) values. Beta-actin was used as the internal control. When the mean value of all the samples was different from that of its control by more than two-fold and p < 0.05, the change in miRNA expression was considered statistically significant.

Next-Generation Sequencing of mRNAs

mRNA expression in the injury and recovery samples (n = 28, each) was determined by GeneTech Biotech Co., Ltd (GeneTech). NGS libraries were constructed using a NEBNext Ultra Directional RNA Library Prep Kit (Illumina). The rRNA was removed from each sample using a Ribo-Zero rRNA removal Kit (Illumina). Then, the rRNA-depleted RNA was fragmented and reverse-transcribed. ProtoScript II Reverse Transcriptase with random primers and actinomycin-D (Illumina) and Second Strand Synthesis Enzyme Mix (Illumina) were used to synthesize the first- and second-strand cDNA, respectively. An AxyPrep Mag PCR Clean-up kit (Axygen, New York, NY, USA) was used to purify the double-stranded cDNA. End Prep Enzyme Mix (Axygen) was used to repair both ends, add a dA-tail, and ligate adaptors to both ends of the cDNA fragments. An AxyPrep Mag PCR Clean-up kit (Axygen) was used to select the size of adaptor-ligated DNA (~360 bp). Uracil-Specific Excision Reagent enzyme (New England Biolabs, Ipswich, MA, USA) was used to digest the dUTP-marked second-strand. Each sample was amplified for 11 PCR cycles using P5 and P7 primers. Libraries were sequenced using a HiSeq 2000 sequencing system (Illumina) using paired-end reads. HiSeq Control Software and OLB with GAPipeline-1.6.0 (Illumina) were used for image analysis and base calling. GENEWIZ (South Plainfield, NJ, USA) processed and analyzed the sequencing data. Trimmomatic v0.30 was used to filter the data in FASTQ format to get high-quality clean data. Fragment counts were determined using Hisat2 v2.0.1 and human reference genome sequences hg19 obtained from the UCSC website. Gene and isoform expression levels were determined from cleaned data using HTSeq v0.6.1.

Analysis of miRNA–mRNA Interactions and Their Functions

miRWalk version 2.0 software (http://mirwalk.uni-hd.de),²⁴ a publicly available comprehensive resource for miRNA–mRNA interaction pairs, was used to explore the miRNA-mRNA interactome. The miRWalk hosts possible binding site interaction information between miRNAs and target mRNAs based on TarPmiR (http://hulab.ucf.edu/research/projects/miRNA/TarPmiR),²⁵ miRNA-target prediction data from TargetScan (https://www.targetscan.org)²⁶ and miRDB (http://www.mirdb.org),²⁷ and validated interaction data from miRTarBase (https://mirtarbase.cuhk.edu.cn).²⁸ All genes were input to the miRWalk website and analyzed using the default settings to produce miRNA–gene interaction output tables. Cytoscape version 3.40 (https://cytoscape.org) was used to visualize the miRNA-mRNA interactome. To explore the function of the mRNA targets identified in the miRNA-mRNA interactome, gene ontology (GO) enrichment analysis²⁹ was performed using the DAVID 2021 database (https://david.ncifcrf.gov/). The KEGG Orthology-Based Annotation System 3.0 (http://kobas.cbi.pku.edu.cn/)³⁰ was used to analyze and illustrate pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Statistical Analysis

R 3.3.3 software (The R foundation) was used to process and analyze the data. Differential expression analysis of miRNAs and mRNAs between the injury and recovery samples was performed using edgeR package (ver 3.15). Benjamini and Hochberg’s approach was performed to control the false discovery rate. miRNA expression from RT-qPCR was analyzed using non-parametric Mann–Whitney U-tests using the Wilcox. test function in R. Statistical significance was indicated by two-sided p < 0.05.

Results

Patient Characteristics

Of the 40 enrolled patients, six died and 34 survived (Figure 1). After excluding patients with failed RNA-seq quality control in either the injury or recovery samples (n = 6), 28 pairs of injured and recovery samples were used to quantify miRNA expression (by RT-qPCR) and mRNA (NGS). After excluding two patients with outlier RNA expression (defined as a data point having less than Q1 - 1.5 × interquartile range (IQR) or higher than Q3 – 1.5 × IQR) in either the injury or recovery samples, 26 patients with both injury and recovery samples were selected for further analysis. Among these 26 patients (Table 1), 14 were male and 12 were female with an average age of 55.8 ± 17.6 years (minimum 22 years and maximum 64 years). Based on injury severity (AIS ≥ 2), the head was the most frequently injured body region (80.8%), followed by the extremities (46.1%), thorax (42.3%), and abdomen (23.0%). The median consciousness level according to the GCS was 12 (Q₁-Q₃, 6–15) and the median ISS was 20 (Q₁-Q₃, 16–25). The time between the collection of injury samples and recovery samples was 24.2 ± 7.2 days.

Table 1 The Patient and Injury Characteristics of the Enrolled and Studied Cohorts

Variables	Enrolled Total Patients n = 40	Studied Patients n = 26
Age, years	60.0±17.8	55.8±17.6
Male, n (%)	22 (55.0)	14 (53.8)
Injured regions, AIS ≥ 2
Head, n (%)	33 (82.5)	21 (80.8)
Face, n (%)	2 (5.0)	1 (3.8)
Thorax, n (%)	19 (47.5)	11 (42.3)
Abdomen, n (%)	9 (22.5)	6 (23.0)
Extremity, n (%)	17 (42.5)	12 (46.1)
External, n (%)	1 (2.5)	0 (0)
GCS, median (IQR)	10 (5–14)	12 (6–15)
ISS, median (IQR)	23 (16–28)	20 (16–25)
16–24	16 (40.0)	14 (53.8)
≥25	24 (60.0)	12 (46.2)

Abbreviations: AIS, abbreviated injury scale; CI, confidence interval; GCS, Glasgow Coma Scale; IQR, interquartile range; ISS, injury severity score; OR, odds ratio.

miRNA and mRNA Expression

NGS analysis of miRNAs revealed that twelve miRNAs (hsa-let-7e-5p, hsa-miR-16-2-3p, hsa-miR-16-5p, hsa-miR-185-5p, hsa-miR-192-5p, hsa-miR-197-3p, hsa-miR-23a-3p, hsa-miR-26b-5p, hsa-miR-223-3p, hsa-miR-361-3p, hsa-miR-425-5p, and hsa-miR-485-5p) were more abundant in the circulating T cells of patients in the injury stage, compared to those in the recovery stage. As shown in Table 2 and Figure 2, RT-qPCR revealed 9 significantly upregulated miRNAs (hsa-miR-16-2-3p, hsa-miR-16-5p, hsa-miR-185-5p, hsa-miR-192-5p, hsa-miR-197-3p, hsa-miR-23a-3p, hsa-miR-26b-5p, hsa-miR-223-3p, and hsa-miR-485-5p) in the injury samples than the recovery samples. No miRNAs were downregulated in the injury samples compared to the recovery samples. In addition, 58 mRNAs were significantly downregulated in circulating T cells of patients in the injury stage compared to patients in the recovery stage (Table 3).

Table 2 The Up-Regulated miRNAs Identified from the Next-Generation Sequencing in the Injury and Recovery Samples. The Fold Expression of miRNAs Was Expressed as Mean (Standard Deviation)

miRNAs	Injured	Recovery	Adjusted p-value
hsa-let-7e-5p	1.6 (1.94)	1.13 (0.30)	0.193
hsa-miR-16-2-3p	3.29 (3.54)	1.10 (0.34)	0.002
hsa-miR-16-5p	4.07 (3.99)	1.21 (0.46)	0.005
hsa-miR-185-5p	3.26 (3.64)	1.26 (0.54)	<0.001
hsa-miR-192-5p	3.26 (3.21)	1.14 (0.30)	0.001
hsa-miR-197-3p	3.38 (3.25)	1.30 (0.39)	0.001
hsa-miR-23a-3p	3.80 (3.87)	1.33 (0.68)	0.002
hsa-miR-26b-5p	2.42 (2.69)	1.02 (0.05)	<0.001
hsa-miR-223-3p	3.42 (2.69)	1.02 (0.05)	<0.001
hsa-miR-361-3p	2.03 (3.35)	1.03 (0.05)	0.096
hsa-miR-425-5P	1.19 (0.24)	1.04 (0.07)	0.481
hsa-miR-485-5p	2.10 (1.33)	1.11 (0.13)	0.009

Table 3 Down-Regulated mRNA Targets Identified from the Next-Generation Sequencing Analysis in the Injury Samples Compared to Those in the Recovery Samples

No.	Gene Name	Full Name	NCBI Entrez Gene	Log2 Fold
1	BTNL9	Butyrophilin like 9	153579	−6.553
2	H3C15	H3 Clustered Histone 15	333932	−5.177
3	H3C14	H3 Clustered Histone 14	126961	−5.177
4	PRSS50	Serine Protease 50	29122	−4.238
5	SLED1	Proteoglycan 3, Pro Eosinophil Major Basic Protein 2 Pseudogene	643036	−3.293
6	LRRN3	Leucine Rich Repeat Neuronal 3	54674	−3.269
7	MZB1	Marginal Zone B And B1 Cell Specific Protein	51237	−3.125
8	PDE9A	Phosphodiesterase 9A	5152	−2.964
8	C10orf10	DEPP1 Autophagy Regulator	11067	−2.910
10	NOG	Noggin	9241	−2.846
11	CCL4L2	C-C Motif Chemokine Ligand 4 Like 2	9560	−2.838
12	PRRT2	Proline Rich Transmembrane Protein 2	112476	−2.727
13	PDGFRA	Platelet Derived Growth Factor Receptor Alpha	5156	−2.642
14	ELK2AP	ETS Transcription Factor ELK2A, Pseudogene	2003	−2.536
15	BIVM-ERCC5	BIVM-ERCC5 Readthrough	100533467	−2.490
16	TUBB2A	Tubulin Beta 2A Class IIa	7280	−2.433
17	TNFRSF17	TNF Receptor Superfamily Member 17	608	−2.417
18	SEC14L2	SEC14 Like Lipid Binding 2	23541	−2.274
19	CACNG6	Calcium Voltage-Gated Channel Auxiliary Subunit Gamma 6	59285	−2.214
20	ACHE	Acetylcholinesterase -Cartwright Blood Group	43	−2.171
21	CRACD	Capping Protein Inhibiting Regulator Of Actin Dynamics	57482	−2.154
22	MYBL2	MYB Proto-Oncogene Like 2	4605	−2.069
23	PCSK1N	Proprotein Convertase Subtilisin/Kexin Type 1	27344	−2.062
24	E2F1	E2F Transcription Factor 1	1869	−2.057
25	NPAS2	Neuronal PAS Domain Protein 2	4862	−2.015
26	PARM1	Prostate Androgen-Regulated Mucin-Like Protein 1	25,849	−1.980
27	SLC22A17	Solute Carrier Family 22 Member 17	51310	−1.913
28	PCYT1B	Phosphate Cytidylyltransferase 1B, Choline	9468	−1.890
29	CHCHD6	Coiled-Coil-Helix-Coiled-Coil-Helix Domain Containing 6	94303	−1.875
30	ALOX15	Arachidonate 15-Lipoxygenase	246	−1.810
31	LANCL3	GeneCards Symbol: LANCL3	347404	−1.805
32	SEC14L5	SEC14 Like Lipid Binding 5	9717	−1.794
33	LAMP5	Lysosomal Associated Membrane Protein Family Member 5	24141	−1.753
34	PARD3	Par-3 Family Cell Polarity Regulator	52688	−1.669
35	ADGRB1	Adhesion G Protein-Coupled Receptor B1	575	−1.643
36	ACSL6	Acyl-CoA Synthetase Long Chain Family Member 6	23305	−1.630
37	SIRPG	Signal Regulatory Protein Gamma	55423	−1.587
38	BMP6	Bone Morphogenetic Protein 6	654	−1.510
39	MMD	Monocyte To Macrophage Differentiation Associated	23531	−1.480
40	EGF	Epidermal Growth Factor	1950	−1.470
41	CTLA4	Cytotoxic T-Lymphocyte Associated Protein 4	1493	−1.457
42	MYLK	Myosin Light Chain Kinase	4638	−1.428
43	AFAP1L2	Actin Filament Associated Protein 1 Like 2	84632	−1.379
44	OLFM2	Olfactomedin 2	93145	−1.374
45	GP5	Glycoprotein V Platelet	2814	−1.286
46	GNG11	G Protein Subunit Gamma 11	2791	−1.224
47	TUBB1	Tubulin Beta 1 Class VI	81027	−1.222
48	AMIGO2	Adhesion Molecule With Ig Like Domain 2	347902	−1.208
49	CAVIN2	Caveolae Associated Protein 2	8436	−1.182
50	BEX3	Brain Expressed X-Linked 3	27018	−1.117
51	ENO2	Enolase 2	2026	−1.107
52	TAL1	TAL BHLH Transcription Factor 1, Erythroid Differentiation Factor	6886	−1.079
53	LTBP1	Latent Transforming Growth Factor Beta Binding Protein 1	4052	−1.067
54	GRAP2	GRB2 Related Adaptor Protein 2	9402	−1.061
55	ANKH	ANKH Inorganic Pyrophosphate Transport Regulator	56172	−1.056
56	GNB5	G Protein Subunit Beta 5	10681	−1.008
57	MSI2	Musashi RNA Binding Protein 2	124540	−0.980
58	AEN	Apoptosis Enhancing Nuclease	64782	−0.858

Figure 2 The real-time quantitative reverse transcription polymerase chain reaction (RT-qPCR) to determine the expression of the 12 up-regulated miRNAs, which were identified from the next-generation sequencing analysis in the injury samples compared to those in the recovery samples.

miRNA–mRNA Interactions and Functions

GO enrichment analysis of the 58 mRNAs revealed 6 biological processes, 2 cellular components, and 1 molecular function (Table 4). In the biological process category, “positive regulation of synapse assembly” and “cell adhesion” were the two most enriched terms in the respective categories. For cellular components, “plasma membrane” and “extracellular space” were the most enriched terms, while “transmembrane receptor protein tyrosine kinase activator activity” was the most enriched molecular function term. To explore miRNA–mRNA interactions, we analyzed the 9 upregulated miRNAs and the 58 downregulated mRNAs using miRWalk. We constructed a node graph of the interactions between the 8 miRNAs and 22 mRNAs (Figure 3). Three miRNAs (hsa-miR-185-5p, hsa-miR-197-3p, and hsa-miR-485-5p) were identified as hub miRNAs (defined by connections with more than three mRNAs) that regulate the miRNA-mRNA interactome. There were no connections between hsa-miR-26b-5p and any of the 58 mRNA targets. GO analysis of the 22 mRNAs involved in the miRNA–mRNA interactions indicated that “plasma membrane” (cellular components) was the only enriched term (Table 5). KEGG pathway enrichment analysis suggested that the “chemokine signaling pathway” (pathway ID: hsa04062) was the predominant pathway. Three genes (PARD3, CCL4L2, GNB5) were involved in the pathway (p = 3.6E-04).

Table 4 Gene Ontology Enrichment Analysis of 58 Differentially Expressed mRNAs Between the Samples of Patients in the Injured Stage Vs Recovery Stage

Category	Term	Count	Gene Name	P-value
CC	GO:0005886~plasma membrane	28	ACHE, PARM1, ALOX15, AMIGO2, SIRPG, ENO2, MYLK, LAMP5, CACNG6, ANKH, SLC22A17, GRAP2, CTLA4, TNFRSF17, BTNL9, PDGFRA, AFAP1L2, PRRT2, EGF, CAVIN2, ACSL6, GNG11, GP5, PARD3, ADGRB1, LANCL3, GNB5, PDE9A	5.16E-05
BP	GO:0051965~positive regulation of synapse assembly	3	LRRN3, AMIGO2, ADGRB1	1.28E-02
BP	GO:0007155~cell adhesion	6	ACHE, PARD3, AMIGO2, ADGRB1, SIRPG, GP5	1.93E-02
MF	GO:0030297~transmembrane receptor protein tyrosine kinase activator activity	2	TAL1, EGF	2.15E-02
CC	GO:0005615~extracellular space	11	ACHE, PCSK1N, OLFM2, LRRN3, EGF, CCL4L2, NOG, ADGRB1, ENO2, GP5, BMP6	2.83E-02
BP	GO:0042060~wound healing	3	PDGFRA, ALOX15, NOG	3.24E-02
BP	GO:0045893~positive regulation of transcription, DNA-templated	6	SEC14L2, AFAP1L2, TAL1, EGF, E2F1, NPAS2	3.80E-02
BP	GO:0007171~activation of transmembrane receptor protein tyrosine kinase activity	2	TAL1, EGF	3.88E-02
BP	GO:0000226~microtubule cytoskeleton organization	3	TUBB2A, PARD3, TUBB1	4.99E-02

Table 5 The Enrichment Analysis of the Gene Ontology (GO) Terms of the Down-Regulated mRNAs in the Constructed miRNA-mRNA Interactome. mRNA Interaction

Category	Term	Count	Gene Name	P-value
CC	GO:0005886~plasma membrane	13	PDGFRA, AFAP1L2, PRRT2, SIRPG, ACSL6, ENO2, GP5, MYLK, ANKH, PARD3, LANCL3, CTLA4, GNB5	1.5E-03

Figure 3 The node graph of the interactions between 8 miRNAs and 22 mRNAs identified from miRWalk.

Discussion

The present study indicates that following major trauma, nine miRNAs were significantly upregulated in the circulating T cells of patients in the injury stage compared to those of the same patients in the recovery stage. The miRNA-mRNA interactome consisted of 8 miRNAs and 22 mRNAs and was used to analyze the functions of the involved genes. The GO term, “plasma membrane”, and the KEGG pathway term, “chemokine signaling pathway”, were predominant.

Chemokines are small chemoattractant peptides that provide directional cues for cell trafficking and thus play a pivotal role in host response protection,^31,32 inflammatory responses, immune homeostasis, and cancer progression.³³ Some chemokines are induced during an immune response to recruit immune cells to infection sites and are thus considered pro-inflammatory, while some chemokines are considered to be homeostatic and are involved in regulating cell migration during development or tissue maintenance.^31,32 Impaired chemokine signaling pathway in T cells controlled by miRNA–mRNA interactions may partly explain the immunosuppression in patients following major trauma. Xiao et al found that even after different injuries, there is an apparently fundamental response to severe inflammatory stress characterized by common genomic signatures.³ Regarding the most significantly regulated pathways, injury led to early activation of innate immune responses and simultaneous suppression of adaptive immune responses. Injury severity, clinical outcomes, magnitude of physiological decline, and transfused blood volume only minimally affected these patterns.³ Interestingly, some studies demonstrated activation of chemokine signaling pathway genes in sciatic nerve injury,³⁴ traumatic brain injury,³⁵ and acute and chronic spinal cord injuries.^36,37 Notably, most patients (16 of 26 patients) in this study had traumatic brain injuries. The impact of injuries to various body parts on specific miRNA–mRNA interactions requires further investigation.

Three miRNAs (hsa-miR-185-5p, hsa-miR-197-3p, and hsa-miR-485-5p) were identified as hub miRNAs that regulate the miRNA–mRNA interactions. Hsa-miR-185-5p is involved in the development of Alzheimer’s disease^38,39 and acquired chemotherapy resistance in patients with gastric cancer.⁴⁰ Overexpression of hsa-miR-185-5p occurs in the vitreous of proliferative diabetic retinopathy patients,⁴¹ in the pericardiac adipose tissue of patients with coronary artery disease,⁴² and significantly enhances endothelial cell angiogenesis.⁴³ hsa-miR-197-3p is highly expressed in proliferative diabetic retinopathy patients,^44,45 is involved in the development of colorectal cancer⁴⁶ and non-small cell lung cancer,⁴⁷ and may act as a biomarker for follicular thyroid cancer.⁴⁸ Further, hsa-miR-485-5p is associated with lupus nephritis,⁴⁹ pulmonary tuberculosis,⁵⁰ and cancers such as lung cancer,⁵¹ hepatocellular carcinoma,⁵² oral tongue squamous cell carcinoma,⁵³ and glioblastoma.⁵⁴ None of these three miRNAs have so far been reported to be associated with T cell function or adaptive immunity. Therefore, investigating the impact of these miRNA–mRNA interactions on immune functions may be necessary to understand immune dysfunction in patients with major trauma.

This study has some limitations. First, a small, single-institution patient population was used to study a complex network of miRNA–mRNA interactions. Second, some factors such as damage control, resuscitation, surgery, and drug use may lead to bias. Third, the expression of these circulating molecules may be dynamic in nature and conclusions drawn from two single measurements may not reflect changes over time. Further, although the study population consisted of severe trauma patients with an ISS ≥ 16 and requiring more than 48 h ventilator support, these patients sustained injuries to different body parts, potentially leading to selection bias in the study. However, we believe that identifying the miRNA–mRNA interactions may provide valuable information for understanding immune dysfunction in patients with major trauma. Finally, the results derived from this study depend on which trauma is the predominant injury type, considering that patients with traumatic brain injuries comprised most of the study population.

Conclusion

Our results reveal that, following major trauma, nine miRNAs are significantly upregulated in the circulating T cells of trauma patients in the injury stage compared to those of the same patients in the recovery stage. A miRNA-mRNA interactome consisting of 8 miRNAs and 22 mRNAs is involved in regulating the chemokine signaling pathway after major trauma.

Ethics Approval and Informed Consent

The study was approved by the Institutional Review Board of Chang Gung Memorial Hospital (ref: 201600221B0 and 201701425B0) under the ethical guidelines of the 1975 Declaration of Helsinki. Written informed consent was obtained from all subjects involved in the study.

Disclosure

The authors report no conflicts of interest in this work.

Acknowledgments

We are grateful for support provided by the Biostatistics Center of the Kaohsiung Chang Gung Memorial Hospital.

Additional information

Funding

This work was supported by the Chang Gung Memorial Hospital (CMRPG8L0841 to Cheng-Shyuan Rau and CMRPG8K0362 to Ching-Hua Hsieh).