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.[44]^1^,[45]^2 The cytokine/genomic “storm” that occurs after major trauma is well described[46]^3^,[47]^4 and involves the simultaneous upregulation of genes involved in systemic inflammatory responses and the downregulation of genes involved in adaptive immunity.[48]^3 These immunoinflammatory responses greatly impact the prognosis of patients with major trauma.[49]^5 Prolonged immune suppression may lead to increased susceptibility to secondary infections.[50]^7 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.[51]^8 In addition, the extent of immune dysfunction following major trauma is strongly correlated to patient morbidity and mortality.[52]^6 Immunosuppression after major trauma is characterized primarily using T cell populations of the adaptive immune system.[53]^7^,[54]^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.[55]^10 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.[56]^3 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.[57]^3 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.[58]^3 MicroRNAs (miRNAs) are endogenous small non-coding RNAs that regulate gene expression through post-transcriptional modulation,[59]^11 primarily by binding to the 3′ UTR of target mRNAs, which leads to mRNA degradation or translation inhibition.[60]^12 miRNAs are important regulators of almost all cellular processes,[61]^13 including the immune system.[62]^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.[63]^14 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.[64]^17 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[65]^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;[66]^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 ([67]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. [68]Figure 1 [69]Open in a new tab 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^7 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 ([70]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 ([71]http://mirwalk.uni-hd.de),[72]^24 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 ([73]http://hulab.ucf.edu/research/projects/miRNA/TarPmiR),[74]^25 miRNA-target prediction data from TargetScan ([75]https://www.targetscan.org)[76]^26 and miRDB ([77]http://www.mirdb.org),[78]^27 and validated interaction data from miRTarBase ([79]https://mirtarbase.cuhk.edu.cn).[80]^28 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 ([81]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[82]^29 was performed using the DAVID 2021 database ([83]https://david.ncifcrf.gov/). The KEGG Orthology-Based Annotation System 3.0 ([84]http://kobas.cbi.pku.edu.cn/)[85]^30 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 ([86]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 ([87]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[1]-Q[3], 6–15) and the median ISS was 20 (Q[1]-Q[3], 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) [88]Open in a new tab 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 [89]Table 2 and [90]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 ([91]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 [92]Open in a new tab Figure 2. [93]Figure 2 [94]Open in a new tab 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. 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 [95]Open in a new tab miRNA–mRNA Interactions and Functions GO enrichment analysis of the 58 mRNAs revealed 6 biological processes, 2 cellular components, and 1 molecular function ([96]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 ([97]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 ([98]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 [99]Open in a new tab Figure 3. [100]Figure 3 [101]Open in a new tab The node graph of the interactions between 8 miRNAs and 22 mRNAs identified from miRWalk. 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 [102]Open in a new tab 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,[103]^31^,[104]^32 inflammatory responses, immune homeostasis, and cancer progression.[105]^33 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.[106]^31^,[107]^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.[108]^3 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.[109]^3 Interestingly, some studies demonstrated activation of chemokine signaling pathway genes in sciatic nerve injury,[110]^34 traumatic brain injury,[111]^35 and acute and chronic spinal cord injuries.[112]^36^,[113]^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[114]^38^,[115]^39 and acquired chemotherapy resistance in patients with gastric cancer.[116]^40 Overexpression of hsa-miR-185-5p occurs in the vitreous of proliferative diabetic retinopathy patients,[117]^41 in the pericardiac adipose tissue of patients with coronary artery disease,[118]^42 and significantly enhances endothelial cell angiogenesis.[119]^43 hsa-miR-197-3p is highly expressed in proliferative diabetic retinopathy patients,[120]^44^,[121]^45 is involved in the development of colorectal cancer[122]^46 and non-small cell lung cancer,[123]^47 and may act as a biomarker for follicular thyroid cancer.[124]^48 Further, hsa-miR-485-5p is associated with lupus nephritis,[125]^49 pulmonary tuberculosis,[126]^50 and cancers such as lung cancer,[127]^51 hepatocellular carcinoma,[128]^52 oral tongue squamous cell carcinoma,[129]^53 and glioblastoma.[130]^54 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. Acknowledgments