Profiling the mRNA and miRNA in Peripheral Blood Mononuclear Cells in Subjects with Active Tuberculosis

Abstract

Purpose

To identify candidate hub genes and miRNAs associated with active tuberculosis (ATB) and reveal the potential molecular mechanisms of disease progression.

Patients and Methods

The expression of mRNA and miRNA was evaluated in peripheral blood mononuclear cells (PBMC) from 4 ATB patients and 4 healthy donors (HD) using high throughput sequencing (HTS) and bioinformatics analysis. Moreover, differentially expressed miRNAs were validated with 35 ATB patients and 35 HDs using reverse transcription quantitative real-time PCR (RT-qPCR).

Results

A total of 2658 significantly differentially expressed genes (DEG) including 1415 up-regulated genes and 1243 down-regulated genes were identified in the ATB group compared with HDs, and the DEGs enriched in immune-related pathways, especially in TNF signaling pathway, cytokine–cytokine receptor interaction, mitogen-activated protein kinase (MAPK) signaling pathways and tuberculosis. Additionally, 10 hub genes were acquired according to protein–protein interaction (PPI) analysis of DEGs. Moreover, 26 differentially expressed miRNAs were found in ATB group compared with HDs. Furthermore, RT-qPCR results showed that hsa-miR-23a-5p (P=0.0106), hsa-miR-183-5p (P=0.0027), hsa-miR-193a-5p (P=0.0021) and hsa-miR-941(P=0.0001) were significantly increased in the ATB patients compared with HD group, and the hsa-miR-16-1-3p was significantly decreased (P=0.0032).

Conclusion

Our research provided a characteristic profile of mRNAs and miRNAs expressed in ATB subjects, and 10 hub genes related with ATB were found, which will contribute to explore the role of miRNAs and hub genes in the pathogenesis of ATB, and improve the ability of differential diagnosis and treatment for the disease.

Keywords:

• mRNA
• miRNA
• active tuberculosis
• high throughput sequencing
• RT-qPCR

Introduction

Tuberculosis (TB) is a worldwide disease mainly caused by Mycobacterium tuberculosis (Mtb), which is one of the top ten causes of death in the world.¹ According to the report of the World Health Organization (WHO) in 2019, there are 10 million people falling ill with TB each year, and the morbidity is more than 1.33‰, the annual death toll for Mtb infection is up to 1.5 million,² which lead to serious health problems and socioeconomic burdens. Therefore, it is particularly important to conduct early diagnosis and treatment for TB. In recent decades, some progress has been made in the diagnosis of Mtb infection; however, there are few reliable, simple, and effective ways to definitively diagnose tuberculosis with the exception of T-SPOT.TB assay.³ Although the combined usage, such as rifampicin (RFP), isoniazid (INH), pyrazinamide (PZA), and other medicine could effectively treat TB, resistant tuberculosis has appeared and gradually threatened TB control efforts.⁴

Additionally, T-SPOT.TB assay could determine whether the subjects were currently infected with Mtb with high specificity and sensitivity, and the test results are not influenced by BCG vaccination and host immune function.⁵ Our previous research found that characteristic T cell receptor beta variable (TCRBV) in the peripheral blood T cell repertoire could be used in the differential diagnosis of the subject with active tuberculosis (ATB) or latent TB infection (LTBI).⁶ However, the molecular mechanism of the occurrence and development of Mtb infection has not been fully understood, besides the diagnosis and differential diagnosis of ATB.

MicroRNA (miRNA) is widely found in various eukaryotes such as animal, plant, parasite, and so on,⁷ miRNA is small, evolutionarily conserved, noncoding RNA with length (~22 nts).⁸ Studies have shown that microRNA regulation of life activities, including regulating the release of cytokines, immune cell signal transduction, cell apoptosis and autophagy and so on.⁹ It is a regulator of post-transcriptional regulation of gene expression in many eukaryotes.¹⁰ Studies have found that some miRNAs can be more stable in the serum of tuberculosis patients,¹¹ and some of them will be up-regulated or down-regulated, which is an important surrogate (biomarker) of the diagnosis for TB, or target of effective treatment for TB. Kathirvel et al indicated that changes in circulating miRNA expression levels could be involved in host immune response disorders to tuberculosis.¹² Although some miRNAs have been found to be associated with diagnostic biomarker of ATB, there are few reports about the role of miRNA–mRNA interaction in the pathogenesis of ATB, and about the miRNA related signaling pathways involved in.

In the present study, we used high throughput sequencing (HTS) technique to evaluate the expression profiles of mRNA and miRNA in peripheral blood mononuclear cells (PBMC) from Chinese ATB patients, and analyze the correlation between the miRNA and ATB infection status. Moreover, we constructed the miRNA–mRNA co-regulatory network for ATB patients by combing the miRNAs and mRNAs that expressed significant difference compared with that of healthy donors. Therefore, the results will contribute to study the role of miRNAs and mRNA in the pathogenesis of ATB, and improve the diagnosis and treatment for ATB.

Materials and Methods

Sample Collection

Thirty-nine subjects with active tuberculosis (ATB) and 39 healthy donors (HD) were enrolled in the Infectious Disease Department (IDD) and Health Examination Center (HEC) of the First Affiliated Hospital, Zhejiang University School of Medicine. ATB was identified by sputum or effusion smear and T-SPOT.TB assay, and confirmed on the basis of clinical syndromes and radiological findings, the more detailed diagnostic criteria could be found in previous reports.¹³ Additionally, subjects were excluded if they were co-infected with other infectious diseases, such as cytomegalovirus (CMV) infection, or human immunodeficiency virus (HIV), more exclusion criteria were shown in our previous report.⁶ Healthy donors were individuals who were excluded from LTBI by T-SPOT and other clinical diagnostic indicators. Moreover, individuals with cancer, diabetes, immune-compromised conditions, and allergic diseases were excluded. The additional characteristics of the enrolled subjects at the time of the study are shown in Table 1.

Table 1 Clinical Characteristics of the Participants

ID	T1	T2	T3	T4	H1	H2	H3	H4
Gender	Female	Male	Female	Female	Male	Female	Female	Female
Age	23	24	56	45	47	34	28	27
Cough for>2 weeks	+	+	+	+	-	-	-	-
Fever for>2 weeks	+	-	-	+	-	-	-	-
T-SPOT	+	+	+	+	-	-	-	-
GeneXpert	+	+	+	+	NA	NA	NA	NA
Radiographic features	+	+	+	+	NA	NA	NA	NA
HIV	-	-	-	-	-	-	-	-
HBV	-	-	-	-	-	-	-	-
Cancer	-	-	-	-	-	-	-	-
ADA	18.8	17.3	11.3	11.4	NA	NA	NA	NA
Alb	40.7	34.2	40.3	40.6	49.1	43.1	75.4	53.3
TP	73.7	55.3	67.9	70.8	73.1	72.6	53.3	75.4
RBC	4.14	4.65	3.94	3.65	5.11	4.71	4.33	4.33
HGB	120	116	116	104	153	125	132	132
WBC	7.4	6.2	2.9	6	4.4	7.9	5.4	5.4

Notes: “+” represent “positive”; “-”represent “negative”.

Abbreviations: HIV, human immunodeficiency virus; HBV, hepatitis B virus; NA, not applicable; ADA, adenosine deaminase; Alb, albumin; TP, total protein; RBC, red blood cell; HGB, hemoglobin; WBC, white blood cell.

This study complies with the declaration of Helsinki, and the protocols involving humans were approved by the ethics committee of the First Affiliated Hospital, Zhejiang University School of Medicine. All participants, including ATB subjects and HDs, agreed to participate in the study and provided written informed consent.

PBMC Preparation

About 5mL of peripheral vein whole blood were obtained from each subject, and the peripheral blood mononuclear cells (PBMCs) were isolated using density gradient centrifugation with Ficoll solution according to operating instructions. Simply, the peripheral vein whole blood was diluted with PBS (×1) in 1: 1 ratio, and the mixture was transferred into a 15mL centrifugal tube pre-filled with 5mL Ficoll. The PBMC were carefully collected in the middle layer (flocculent layer) of the tube following centrifugation for 20 min (500× g). Subsequently, the PBMC were washed and centrifuged with PBS for 6 minutes at 350× g, repeated this step once, and the precipitate was collected into 1.5mL EP tube containing 1.0 mL of TRIzol reagent.

RNA Isolation and Sequencing

Total RNA from PBMCs was extracted using the TRIzol reagent (Life Technologies, Carlsbad, CA, USA) according to the manufacturer’s protocol. The extracted total RNA was quantified using a Nanodrop 2000 spectrophotometer (Thermo Fisher, America). The RNA integrity was measured by visualization of the 28S and 18S RNA transcripts on a 1.2% agarose, and the quality of RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Total RNA was reversely transcribed into cDNA, and sequencing library was constructed by Shanghai Majorbio Co. (Shanghai, China) according to the manufacturer’s instructions. Illumina TruSeq Small RNA Library Prep kit and Illumina TruseqTM RNA sample prep kit were used to create small RNA library and mRNA library, respectively. Additionally, the eight sequencing libraries (named T1, T2, T3, T4, H1, H2, H3 and H4) were constructed and sequenced, and the mRNA sequencing was performed on an Illumina Novaseq 6000 machine, the miRNA sequencing was performed on Illumina HiSeq X Ten machine. The miRBase database was used to analyze known miRNAs of ATB and healthy donor groups, and miRDeep2 was used to predict new miRNA. The more detailed protocol was shown in a previous report.¹⁴

Identification of DEGs and DEMs

The count data of mRNA and miRNA sequencing were subjected to principal component analysis (PCA) in R software (version 3.6.1, http://www.R-project.org) to detect outliers. Significantly differentially expressed genes (DEGs) and differentially expressed miRNA (DEMs) were screened based on the DESeq2 package in R software. The raw P value was adjusted by the Benjamini and Hochberg method (BH), and the |log2FC| ≥ 1 and adjusted P value <0.05 were regarded as the cutoff criteria and significantly.

Functional Gene Enrichment

To further understand the function of identified DEGs, Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed. The Goatools software (https://github.com/tanghaibao/goatools) was used to perform GO enrichment analysis, and the R program to perform KEGG pathway enrichment analysis. P value<0.05 was considered statistically significant.

Construction of PPI Network and Module Analysis

The STRING database (version: 11.0, https://string-db.org) is a biological database designed to analyze the protein–protein interaction (PPI) network. The brief steps were the DEGs were mapped to STRING to evaluate the interactive relationships, and the cutoff criterion was defined as more than 0.4 (medium confidence). Finally, Cytoscape (version: 3.7.1, www.cytoscape.org) was used for network construction, and Hub genes were ranked by the CytoHubba plugin. The sub-network modules in the PPI network were obtained by the plug-in unit (Molecular Complex Detection, MCODE). An MCODE score >4 and a number of nodes >5 were used as the criteria to define a module. The default parameters were set as node score cutoff=0.2, degree cutoff=2, k-core=2 and max depth=100.

Construction of miRNA–mRNA Network

To explore the potential functions of the significantly differentially expressed miRNAs (DEMs), the Miranda (http://www.miranda.org/) and RNAhybrid (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid/) algorithms were used to predict their target genes. We selected significant DEMs and their corresponding target genes, which were also identically expressed in our test of mRNA expression (DEGs), to construct the miRNA–mRNA interactive network by Cytoscape. A group of up-regulated miRNAs and corresponding down-regulated target mRNAs and another group of down-regulated miRNAs and corresponding up-regulated target mRNAs were determined. Finally, we constructed the regulatory network based on the negative correlation of miRNA and its target gene.

RT-qPCR for miRNA Verification

The RNA of PBMC from 35 ATB patients and 35 HDs were used to verify the DEMs with RT-qPCR. The simple operation is as follows, cDNA was synthesized using the PrimeScriptTM RT reagent Kit with gDNA Eraser (Takara, Dalian, China). RT-qPCR analysis was performed using SYBR Premix Ex Taq (Takara) on the Applied Biosystems 7500 (Applied Biosystems, CA, USA). RT-qPCR reactions were performed on the Applied Biosystems 7500 (Applied Biosystems, Foster City, CA). The expression of miRNAs was determined using the Bulge-LoopTM miRNA RT-qPCR primer kit (RiBoBio, Guangzhou, China). U6 was used as reference gene. 2-ΔΔCt method was used to quantify qPCR data.

Statistical Analysis

Statistical analysis was performed with GraphPad Prism 7 software (GraphPad Software, La Jolla, CA, USA). The nonparametric Mann–Whitney test was used to evaluate the differences between groups. The miRNA expression value in the subject with healthy donors was normalized to 1. P value<0.05 was considered statistically significant.

Results

Sequencing Information and Quality Control

There are two sets of high throughput sequencing (HTS) data (mRNA, and miRNA) obtained at the present study. Sequencing readings of each sample and their profile of mapping to the genome were summarized in Table S1 and S2. Moreover, the results of principal component analysis (PCA) for the sequencing results are shown in Figure S1.

Identification of DEGs

According to the above results, the 8 samples were eligible for differential expression analysis. Differentially expressed genes (DEGs) were screened with fold change (|log2FoldChange| > 1) and adjusted P (P value < 0.05) determined by DESeq2 analysis. In comparison with the HD group, we identified 1415 up-regulated DEGs (uDEGs) and 1243 down-regulated DEGs (dDEGs) in the ATB individuals. The volcano plot of DEGs is shown in Figure 1. Additionally, the gene expression of ATB and HD samples were significantly distinguished in the heat-map clustering of DEGs (Figure S2).

Figure 1 Volcano plot of differentially expressed genes (DEGs). Red, significantly up-regulated gene; green, significantly down-regulated gene; gray, non-significantly differently expressed gene.

Gene Ontology and KEGG Analysis of DEGs

Top 20 Gene Ontology including biological process, cellular component and molecular function were analyzed. Go terms with adjusted P value <0.05 were selected and ranked by adjusted P value. The GO analysis results showed that up-regulated DEGs were predominantly enriched in the following biological processes (BP): negative regulation of leukocyte activation, cellular response to extracellular stimulus and negative regulation of secretion. The down-regulated DEGs were mainly involved in DNA deamination and small molecule metabolic process in biological processes. As for cellular component (CC), the down-regulated DEGs were enriched in intracellular organelle and intracellular membrane-bounded organelle, but there were no up-regulated DEGs enriched in this subgroup. When regarding molecular function (MF), the down-regulated DEGs were enriched in transcription regulator activity, DNA-binding transcription factor activity and DNA binding. Similarly, no up-regulated DEGs were enriched in this subgroup either (Figure 2A and B).

Figure 2 Bubble diagram of top 20 enriched GO terms and KEGG pathways for significantly differentially expressed genes (DEGs) between ATB and HD. (A) Top 20 ranked GO terms of up-regulated DEGs. (B) Top 20 ranked GO terms of down-regulated DEGs. (C) Top 20 pathways of up-regulated DEGs. (D) Top 20 pathways of down-regulated DEGs. Rich factor: the ratio of genes enriched in the GO term (KEGG) to the number of annotation genes, the enrichment degree was stronger with a bigger Rich factor. The size of dots indicates the number of genes in the GO term (KEGG).

Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The results of KEGG pathway analysis showed that the up-regulated DEGs were significantly enriched in MAPK signaling pathway, cytokine–cytokine receptor interaction, TNF signaling pathway, TB (tuberculosis) and peroxisome proliferator-activated receptor (PPAR) signaling pathway. Most of the down-regulated DEGs were related to cell adhesion molecules (CAMs), antigen processing and presentation, and T cell receptor signaling pathway (Figure 2C and D).

PPI Generation and Module Analysis

Top 1000 up-regulated DEGs and top 1000 down-regulated were added to the PPI network, including 1429 nodes and 7134 edges. Hub genes were ranked by cytoHubba of Cytoscape, and found CUL3, KLHL2, KLHL22, SPSB2, KBTBD8, FBXO7, KCTD7, FBXO40, FBXL16 and FBXL8 as hub genes. Cytoscape MCODE was applied to screen the modules of the PPI network. Besides, top 2 primary modules from the PPI network were selected for functional enrichment analysis. Module 1 consists of 20 nodes and 190 edges, and Module 2 comprises 16 nodes and 96 edges (Figure 3). Enrichment analyses for Module 1 demonstrated that the pathways were markedly enriched in ubiquitin-mediated proteolysis and Hedgehog signaling pathway. The most significantly enriched BP in module 1 included protein poly-ubiquitination, ubiquitin-dependent protein catabolic process and modification-dependent macromolecule catabolic process (Figure S3). Additionally, the KEGG pathways in Module 2 were mainly associated with Chemokine signaling pathway, cytokine–cytokine receptor interaction and Neuroactive ligand–receptor interaction. The most significantly enriched BPs in Module 2 was response to lipopolysaccharide, response to molecule of bacterial origin and regulation of signaling (Figure S4).

Figure 3 Top 2 modules of PPI sub-network by MCODE in Cytoscape software. (A) module 1; (B) module 2.

Significant Difference of miRNA Expression

There are 9 up-regulated miRNAs and 17 down-regulated miRNAs in ATB compared with HD group with differential criteria (|log2FC|>1, adjusted P value<0.05). Among the 26 DEMs, 22 miRNAs are known, and 4 miRNAs are novel (Table 2). Additionally, 428 target genes were identified for 9 up-regulated miRNAs, and 55 target genes were identified for 8 down-regulated miRNAs (no target gene found for another 9 down-regulated miRNAs) by Miranda and RNAhybrid algorithms.

Table 2 Differentially Expressed miRNAs (DEM) in ATB Compared with HDs

miRNA	Type	Log2FC(ATB/HD)	P-Adjust	Regulate
hsa-miR-941	Known	1.12891203	0.005289709	Up
hsa-miR-760	Known	1.358401543	0.048009942	Up
hsa-miR-23a-5p	Known	1.21548475	0.002130772	Up
hsa-miR-193a-5p	Known	1.049993343	0.003776702	Up
hsa-miR-183-5p	Known	1.743313601	0.028929263	Up
hsa-miR-12,136	Known	1.751294075	0.009083745	Up
6_14541	Novel	2.11648939	0.023954136	Up
3_8196	Novel	1.005929855	0.048009942	Up
2_6381	Novel	3.26489339	0.03996549	Up
hsa-miR-1277-3p	Known	−2.206447856	8.89431E-05	Down
hsa-miR-33a-5p	Known	−2.328391239	0.002130772	Down
hsa-miR-1277-5p	Known	−2.10788954	0.00253324	Down
hsa-miR-19a-3p	Known	−1.27165401	0.005289709	Down
hsa-miR-1307-5p	Known	−1.998188963	0.005289709	Down
hsa-miR-652-5p	Known	−1.165243492	0.005989237	Down
hsa-miR-19b-3p	Known	−1.223579915	0.007794562	Down
hsa-miR-29c-3p	Known	−1.214234207	0.016080995	Down
hsa-miR-33a-3p	Known	−1.05087557	0.017908658	Down
hsa-miR-590-5p	Known	−1.727717534	0.023835661	Down
hsa-miR-590-3p	Known	−1.01612367	0.023954136	Down
hsa-miR-101-3p	Known	−1.248748508	0.039055415	Down
hsa-miR-16-1-3p	Known	−1.954883466	0.042332706	Down
hsa-miR-1537-3p	Known	−2.892938848	0.042332706	Down
hsa-miR-340-5p	Known	−1.12770789	0.048009942	Down
hsa-miR-545-5p	Known	−1.575959365	0.048009942	Down
1_592	Novel	−2.313761221	0.002236797	Down

miRNAs–mRNAs Network

The intersection genes between DEGs (mRNA) and the target gene of DEMs (miRNA) were determined by the Venn diagram, and found 4 up-regulated genes, and 19 down-regulated genes (Figure 4). Subsequently, the DEMs and intersection genes were used to construct the miRNA-mRNA network, and there were 5 up-regulated DEMs (miRNA) corresponding to 19 down-regulated DEG (mRNAs), and 1 down-regulated DEMs (miRNA) corresponding to 4 up-regulated DEG (mRNAs) (Figure 5).

Figure 4 Intersection genes between DEGs and target gene of DEMs. DEGs, significantly differentially expressed genes (mRNA) between ATB and HD; DEMs, significantly differentially expressed genes (miRNA) between ATB and HD.

Figure 5 MiRNA-mRNA regulatory network was analyzed using Cytoscape software. Diamond-shaped represents miRNAs, and circular nodes represent target genes. Red, up-regulated miRNA or mRNA; green, down-regulated miRNA or mRNA in active tuberculosis subjects.

Validation of DEMs by RT-qPCR

To confirm the DEMs resulted from HTS and Bioinformatics analysis, 5 DEMs (hsa-miR-23a-5p, hsa-miR-193a-5p, hsa-miR-941, hsa-miR-183-5p, hsa-miR-16-1-3p) were selected and validated by RT-qPCR analysis with samples of 35 ATB patients and 35 healthy donors. The results were consistent with the prediction, the expression of 4 out of the 5 miRNAs were significantly increased in the ATB patients compared to healthy donors; however, the another one (hsa-miR-16-1-3p) was significantly decreased (P=0.0026) (Figure 6).

Figure 6 We analyzed the expression of five miRNAs (hsa-miR-23a-5p, hsa-miR-183-5p, hsa-miR-193a-5p and hsa-miR-941 and hsa-miR-16-1-3p) in 35 active tuberculosis patients versus 35 healthy donors. U6 was used as reference gene. 2-ΔΔCt method was used to quantify RT-qPCR data. Mann–Whitney test was used to evaluate the differences between groups. P value<0.05 was considered statistically significant.

Abbreviation: RT-qPCR, reverse transcription quantitative real-time PCR.

Discussion

Mycobacterium tuberculosis (Mtb) infection is a serious health problem on worldwide among the young, elderly, and immunocompromised population.¹⁵ The transcriptional profile of PBMC in the subjects with Mtb infection are complicated and could be due to multiple factors, such as age and genetic background.¹⁶ In the present study, we examined the transcriptomes of ATB patients and healthy individuals (HDs) and uncovered specific molecular features and pathways associated with ATB.

Mtb infection results in the complex intracellular bacterial infection, and the control of Mtb infection requires innate and adaptive immune responses.^17,18 The cell-mediated immune response is considered to play a vital role in controlling Mtb.¹⁷ When immune cells are activated by Mtb-specific antigens, they can release proinflammatory cytokines, such as TNF-α, IFN-γ and IL-17. Interestingly, in the study, the KEGG pathway analysis showed that DEGs were enriched in TNF signaling pathway, cytokine–cytokine receptor interaction, IL-17 and MAPK signaling pathway. In previous research, TNF-α was found to initiate innate cytokine response and phagocyte activation in Mtb infection.¹⁹ In addition to controlling Mtb infection, TNF-α can also cause serious tissue damage.²⁰ IL-17 is essential for survival following Mycobacterium tuberculosis infection.²¹ Recently, Choi demonstrated that the expansion of IFN-γ/IL-17-producing cells in the lung is of great significance for the protective effect in the tuberculosis subunit vaccine.²² Su et al found that Rv1016c protein inhibits the Th17 cell differentiation normally polarized by BCG-infected DCs and decreased the production of IL-17. As for the function of the MAPK signaling pathway, some studies have suggested that MMPs induce lung tissue remodeling and regulate granuloma formation after Mtb infection.^23,24 Various MMP inhibitors have been used to study the immunomodulatory role of MMPs in Mtb infection and showed promising effect.^25,26 For example, Xu et al reported that inhibition of matrix metalloproteinase (MMP) enhances the in vivo potency of isoniazid and rifampicin.²⁷ In summary, our results revealed a series of processes and pathways related to Mtb infection, which provide a new basis for further explore the pathogenesis of ATB.

According to PPI analysis of DEGs, CUL3, KLHL2, KLHL22, SPSB2, KBTBD8, FBXO7, KCTD7, FBXO40, FBXL16 and FBXL8 were identified as hub genes and most of them are related to the process of ubiquitination. Ubiquitination is very important in regulating eukaryotic cell processes, including innate immune system responses and cell cycle progression.^28,29 Some findings also reveal the close connection between ubiquitination and the biological process of tuberculosis. Cheng et al indicated that the host cytosolic RNA sensing pathway stimulates ICAM-1 expression through the inhibition of CRL4^COP1/DET1, an E3 ubiquitin ligase, and promotes T Lymphocyte-mediated mycobacterial killing in macrophages.³⁰ Jia et al found that Gal9 can be recruited by impaired lysosomes to regulate ubiquitination and contribute to autophagic control during Mtb infection.³¹ Additionally, Chai et al demonstrated that Mtb surface protein Rv1468c could recruit ubiquitin to trigger host xenophagy.³² Furthermore, miR-325-3p was found to target the ubiquitin ligase LNX1, and leads to abnormal accumulation of NEK6, thus inhibiting the cell apoptosis of Mycobacterium tuberculosis.³³ Our results are consistent with other’s reports. In the study, we recognized some important genes related to ubiquitin which may regulate the biological process in Mtb infection.

Module analysis of the PPI network indicated that ATB infection was related to ubiquitin-mediated proteolysis, hedgehog signaling pathway, cytokine–cytokine receptor interaction, chemokine signaling pathway and neuroactive ligand–receptor interaction. Among them, ubiquitin-mediated proteolysis and cytokine–cytokine receptor interaction pathway have been discussed above. In addition, the chemokine signaling pathway, neuroactive ligand–receptor interaction and hedgehog signaling pathway are important pathways that should be concerned about. Chemokines play a central role in cell recruitment to the site of Mycobacterium tuberculosis infection and within the granuloma.³⁴ Interestingly, a study reported that Mtb Heat-Shock Protein 16.3 promotes macrophages to M2 phenotype via CCRL2/CX3CR1.³⁵ In an in vitro model of central nervous system tuberculosis, the Hedgehog pathway was down-regulated and impaired blood-brain barrier function.³⁶ Neuroactive ligand–receptor interaction has been reported to be associated with various cancers including lung cancer by pathway analysis.³⁷

The important function of microRNA is to regulate the expression of target genes.³⁸ Each miRNA can have multiple target genes and affect the synthesis of multiple proteins. Several miRNAs can also jointly regulate the same target gene to promote or inhibit the process of translation. There are two main modes of action: one is that the target gene is not completely complementary to inhibit the expression of the target gene without affecting its stability; the other is that it is completely complementary with the target gene, leading to the degradation of the target gene.³⁹ Therefore, the analysis of candidate genes and pathways related to ATB can provide a cognitive basis for the occurrence of diseases. Moreover, miRNAs are key in both adaptive and innate immunity by regulating the differentiation of various immune cell subsets as well as their immunological functions.⁴⁰ In the present study, we identified 26 miRNAs significantly different between ATB and HD samples. Function of differentially expressed miRNAs (DEMs) are still not well understood, such as miR-12,316 and miR-1537. However, some miRNAs, such as hsa-miR-23a-5p, have been connected to specific signaling pathways.⁴¹ In response to Mtb infection, miR-23a-5p regulates Mtb survival and autophagy via TLR2/MyD88/NF-κB pathway by targeting TLR2.⁴² Consistent with our sequencing results, hsa-miR-101 was significantly decreased in the serum of pulmonary TB patients compared with healthy controls.⁴³ High miR-183 level was found to be correlated with macrophage activity in the serum of TB patients.⁴⁴ In addition, Yi et al reported that hsa-mir-19b-2 was significantly down-regulated in TB group compared with controls.⁴⁵ The expression of has-mir-16 in serum was significantly down-regulated in the ATB group compared with the BCG-inoculated group.⁴⁶ In this study we validated 5 significantly different miRNA (hsa-miR-23a-5p, hsa-miR-193a-5p, hsa-miR-941, hsa-miR-183-5p and hsa-miR-16-1-3p) using RT-qPCR analysis. To our knowledge, hsa-miR-193a-5p and hsa-miR-941were first reported to be significantly elevated in ATB compared with that of HD control. Furthermore, we constructed the miRNA–mRNA regulatory network to explore the function of DEMs. The result indicated that hsa-miR-23a-5p/AIFM2, hsa-miR-193a-5p/ZSCAN20, hsa-miR-941/OLFM2 regulatory axes pay a key role in ATB onset and progression, although, which should be further confirmed in the future with more cases.

Conclusion

In summary, the characteristic profiles of mRNA and miRNA in PBMCs were determined in patients with ATB infection, and some predicted targets of characteristic miRNA are identical to our determined genes, although this should be confirmed in the future with more cases. Our results could provide a theoretical basis for studying the role of miRNAs in the pathogenesis of ATB and indicate the different mRNAs and miRNAs contribute to improve the diagnosis and treatment for ATB infection and distinguish the different clinical status of Mtb infection.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Acknowledgments

We sincerely thank all of the subjects participating in the study.

Disclosure

The authors report no conflicts of interest in this work.