Long non-coding RNAs ENST00000429730.1 and MSTRG.93125.4 are associated with metabolic activity in tuberculosis lesions of sputum-negative tuberculosis patients

Differentially expressed mRNAs and differentially expressed lncRNAs

386 mRNAs (141 down-regulated and 245 up-regulated mRNAs) were differentially expressed in PET-high samples when compared to PET-low samples (Figure 1A, 1B). Hierarchical clustering of these differentially expressed mRNAs indicated that PET-high samples were indeed distinct from PET-low samples (Figure 1C). Moreover, 44 differentially expressed lncRNAs (24 down-regulated and 20 up-regulated lncRNAs) were identified between PET-high and PET-low samples (Figure 1D, 1E). Hierarchical clustering of differentially expressed lncRNAs again revealed that PET-high and PET-low samples were clustered in two distinct groups (Figure 1F).

GO functional enrichment analysis indicated that differentially expressed mRNAs were significantly related to the fibrillar collagen (enrich factor = 23.22), platelet-derived growth factor binding (enrich factor = 22.87), and protein heterotrimerization (enrich factor = 20.13) (Figure 2A). Pathway enrichment analysis suggested these differentially expressed mRNAs were involved in the pathways of ribosome (enrich factor = 6.38), systemic lupus erythematosus (enrich factor = 5.87), and staphylococcus aureus infection (enrich factor = 4.99) (Figure 2B).

LncRNA-mRNA co-expression network

The correlations between differentially expressed lncRNAs and differentially expressed mRNAs were revealed to construct the lncRNA-mRNA co-expression network. The lncRNA-mRNA co-expression network was consisted of 1956 positive correlations and 6 negative correlations involving 33 differentially expressed lncRNAs and 340 differentially expressed mRNAs. Five up-regulated lncRNAs, including MSTRG.93125.4 (40 degrees), ENST00000429730.1 (35 degrees), MSTRG.93124.1 (31 degrees), ENST00000602711.1 (21 degrees), and MSTRG.177254.3 (17 degrees) correlated with more than 15 targeted mRNAs (Figure 2C). Regarding mRNAs, eight up-regulated mRNAs had a degree number larger than 35. They were C15orf48 (41 degrees), RASSF4 (38 degrees), CCL5 (37 degrees), CYBB (37 degrees), PGD (37 degrees), LAP3 (36 degrees), LAPTM5 (36 degrees), and LGI2 (36 degrees) (Figure 2C).

The differentially expressed mRNAs in lncRNA-mRNA co-expression network were enriched with the fibrillar collagen, platelet-derived growth factor binding, protein heterotrimerization, regulation of dendritic cell differentiation, and leukocyte migration involved in inflammatory response (Figure 2D). Pathway enrichment analysis for differentially expressed mRNAs in lncRNA-mRNA co-expression network indicated that the top ranked pathways were involved systemic lupus erythematosus, Staphylococcus aureus infection, phenylalanine metabolism, Leishmaniasis, and asthma (Figure 2E).

Target gene prediction of differentially expressed lncRNAs and functional analysis

One cis-target gene and 440 trans-target genes were selected from the lncRNAs differentially expressed in PET-high lung tissues (See Materials and Methods). The subset of these genes that correlated significantly with the lncRNAs (|R| > 0.85 and FDR < 0.0005) were predicted as target genes and chosen for further functional and pathway enrichment analysis. All of these lncRNA target genes were differentially expressed (Supplementary Table 4). These predicted target genes were mainly involved in the biological processes of chemokine-mediated signaling pathway, collagen catabolic process, neutrophil chemotaxis, neutrophil migration, and collagen metabolic process (Figure 3A). KEGG pathway enrichment analysis indicated the predicted target gene subset were involved in the pathways of systemic lupus erythematosus, rheumatoid arthritis, Staphylococcus aureus infection, chemokine signaling pathway, and cytokine-cytokine receptor interaction (Figure 3B).

We then visualized the top 30 KEGG pathways, the predicted target genes, as well as the differentially expressed lncRNAs together in a lncRNA-target-pathway network map (Figure 3C). 5 lncRNAs (ENST00000429730.1, MSTRG.93125.4, MSTRG.93124.1, ENST00000602711.1, and ENST00000555688.1) and 9 target genes (HLA-DMB, RASSF4, CCL5, CXCL9, C15orf48, CCR5, CYBB, LGI2, and PGD) had a degree of number more than 25. These extensively connected genes were involved in pathways that were enriched with connections: cytokine-cytokine receptor interaction (involving CCL5, CXCL9 and CCR5), chemokine signaling pathway (involving CCL5, CXCL9 and CCR5), and cell adhesion molecules (involving HLA-DMB and HLA-DPB1). RASSF4 regulates PIP2 production, a critical step for phagocytosis and phagosome formation. C15orf48 encodes a novel subunit of the mitochondrial electron transport chain (32045714). Leucine-rich repeat glioma inactivated gene member encoded by LGI2 plays an important role in remodeling of synapses between nerve cells (23713523). PGD encodes 6-phosphogluconate dehydrogenase in the pentose phosphate pathway, which converts 6-phosphogluconate into ribulose 5-P, a metabolite product that supports inflammatory macrophage responses (25904920). Other highly-connected pathways in the lncRNA-target-pathway network map included complement and coagulation cascades (involving C1QB and C1QA), lysosome (involving CD68, CTSD, and CTSK), phagosome (involving CYBB and CYBA). Four highly-connected target genes CCL5, CXCL9, CYBB, and HLA-DMB were predicted target genes of two highly-connected lncRNAs ENST00000429730.1 and MSTRG.93124.1 with 54 and 53 predicted target genes, respectively.

Expression validation and prediction efficiencies of lncRNA indicators

Expression levels of 7 mRNAs (C1QB, CD68, CCL5, CCL19, MMP7, HLA-DMB, and CYBB) and 4 lncRNAs (ENST00000429730.1, MSTRG.93125.4, MSTRG.93124.1, and ENST00000602711.1) that had relatively larger number of degree in lncRNA-target-pathway network were measured in the lung PET-high and PET-low tissue samples collected from 11 patients. The expression levels of C1QB, CD68, CCL5, CCL19, MMP7, HLA-DMB, and CYBB were validated to be significantly up-regulated in lung PET-high tissue samples (Figure 4A). Obviously higher expressions of two lncRNAs ENST00000429730.1 and MSTRG.93125.4 were found in lung PET-high tissue samples (Figure 4A). ROC curve analysis of ENST00000429730.1 showed that the AUC was 0.750 with a sensitivity of 58.33% and a specificity of 100.00% (p = 0.024) (Figure 4B), and the AUC of MSTRG.93125.4 was 0.813 with a sensitivity of 58.33% and a specificity of 100.00% (p = 0.001) (Figure 4C).

Bioinformatics analysis of ENST00000429730.1 and MSTRG.93125.4

Genomic location of ENST00000429730.1 was predicted to be on chromosome 2q33.3 and assumed to be unable to encode genes (Figure 5A). The optimal secondary structure for ENST00000429730.1 had several hairpin loops with a minimum free energy (MFE) of -110.50 kcal/mol (Figure 5B). There were strong interactions of ENST00000429730.1 with proteins of ELAV1, TRA2B, PTBP1, SRSF9, and SRS10. The interaction between ENST00000429730.1 and ELAV1 has an interaction propensity of 68, and with the discriminative power of 97% (Figure 5C).

Genomic location of MSTRG.93125.4 was predicted to be on chromosome 14KI270846v1_alt and unable to encode genes (Figure 6A). The optimal secondary structure for MSTRG.93125.4 had the MFE of -359.20 kcal/mol (Figure 6B). For MSTRG.93125.4, strong interactions were found with proteins of SFPQ, PCBP1, LN28B, SRSF2, and PCBP2. The interaction between MSTRG.93125.4 and SFPQ has an interaction propensity of 180, with the discriminative power of 100% (Figure 6C). Prediction analysis of potential transcriptional factors showed that MSTRG.93125.4 can combine with 67 transcriptional factors (such as Kid3, AML1, myogenin, MyoD, and AP-4), and ENST00000429730.1 can combine with 82 transcriptional factors (such as AML1, TBX5, E2F, IPF1, and E2F-1) (Figure 7).

Patients and samples

Five patients who diagnosed as sputum-negative pulmonary tuberculosis in our hospital from June 06 to December 2018 were enrolled for RNA-sequencing analysis. There were one male and four females, aged 20-38 years, with an average age of 31.2 ± 7.4 years (Supplementary Table 1). Lung tissue samples with high metabolic activity (PET-high) and low metabolic activity (PET-low) demonstrated by (FDG-PET/CT) were collected from these five participants for the subsequent sequencing (Figure 8). Moreover, lung PET-high and PET-low tissue samples collected from 11 enrolled subjects were used for experimental validation of expression changes by quantitative real-time PCR (qRT-PCR). There were 7 male and 4 females, aged 20-57 years, with an average age of 38.3 ± 11.3 years (Supplementary Table 2).

All the subjects were confirmed as pulmonary tuberculosis by pathological examination with sputum smear and culture negative preoperatively, and accompanied by tuberculous pleurisy or other extrapulmonary tuberculosis. The included participants were demonstrated to have residual metabolic activity by FDG-PET/CT, and in non-human immunodeficiency virus (non-HIV) infection status. Patients were excluded when they met one of the following criteria: sputum positive for mycobacterium tuberculosis complex by Xpert or culture (bacteriologically confirmed tuberculosis disease) before surgery; in the presence of tuberculosis symptoms; signs or symptoms of acute illness; with the diagnosis of malignancy; any clinical condition requiring systemic steroid or other immunosuppressive medication in the preceding six months; uncontrolled diabetes mellitus, pregnant or breastfeeding; anemia (hemoglobin < 7g/dL); significant smoking history (> 30 pack-years). All patients signed informed consent to participate in the study. This study was approved by the Ethics Committee of the Shanghai Public Health Clinical Center.

RNA-sequencing by Illumina HiSeq

Total RNA of each lung sample was extracted by using TRIzol Reagent (Invitrogen, US) and further purified using the RNeasy Mini Kit (Qiagen, Germany) according to the manufacturer’s instructions. Total RNA of each sample was quantified and qualified by Agilent 2100 Bioanalyzer (Agilent Technologies, USA), NanoDrop (Thermo Fisher Scientific Inc., USA) and 1% agarose gel. Total RNA (1 μg) with the value of RNA integrity number (RIN) > 7 was used for library preparation of the following next generation sequencing, with the NEBNext Ultra II Directional RNA Library Prep. Kit for Illumina. The rRNA was depleted from total RNA using Ribo-Zero™ rRNA removal Kit (Illumina, USA), and cDNA libraries were generated following manufacture’s protocols. Then, libraries with different indices were multiplexed and loaded on Illumina HiSeq instrument (Illumina, USA) for 2 × 150 paired-end sequencing according to manufacturer’s instructions in the Medical Laboratory of Nantong ZhongKe.

RNA-sequencing data analysis

Generated raw reads were saved as fastq format. Raw reads were processed by Trimmomatic (version 0.30) to filter dirty reads that with sequence adaptors, primers, fragments, or length of bases less than 20. Afterwards, high quality clean data were aligned to reference genome via software Hisat2 (version 2.0.1). Comparing with the PET-low samples, the mRNAs differentially expressed in PET-high samples were uncovered by the DESeq Bioconductor package. The cut-off criteria were set as |log₂ fold change (FC) | > 1 and q value < 0.05, after adjusted by Benjamini and Hochberg’s approach. Hierarchical clustering analysis based on centered Pearson correlation was performed for all the differentially expressed mRNAs via the pheatmap (version 1.0.8) [35, 36]. Meanwhile, differentially expressed lncRNAs between PET-high and PET-low samples were identified by the DESeq Bioconductor package, with the thresholds of |log₂ FC| > 0.5 and false discovery rate (FDR) < 0.05.

In order to identify functions and specific signaling pathways of differentially expressed genes, Gene Ontology (GO) and KEGG pathway analyses were performed. The thresholds were set as gene number ≥ 5 and enrichment test p-value ≤ 0.01 [37, 38].

Construction of lncRNA-mRNA co-expression network

The correlations of differentially expressed lncRNAs with differentially expressed mRNAs were evaluated by calculating the Pearson correlation coefficient R values. Fisher's asymptotic test implemented in the WGCNA library of R was adopted to estimate p-values for each correlation pair. The p-value was further adjusted into FDR by Bonferroni multiple test correction implemented in the multtest package of R [39]. The lncRNA-mRNA correlations with the |R| > 0.9 and FDR < 0.005 were regarded as co-expressed and used to construct lncRNA-mRNA co-expression network. The co-expression network was visualized by Cytoscape software. GO functional and KEGG pathway enrichment analysis were conducted for the differentially expressed mRNAs in the lncRNA-mRNA co-expression network, with the criterion of p < 0.05.

LncRNA target prediction and regulatory network construction

Target genes of differentially expressed lncRNAs were predicted via cis- or trans-regulatory effects. LncRNAs and potential target genes were paired and visualized using University of California Santa Cruz (UCSC) genome browser (http://genome-asia.ucsc.edu/index.html). Genes transcribed within a 10 kbp window upstream or downstream of lncRNAs and co-expressed with the lncRNA genes were considered as cis target genes using Bedtools (version 2.17.0). The algorithm for regulation in trans is based on mRNA sequence complementarity and then RNA duplex energy prediction. First, BLAST (version 2.2.28+) software was used sequence complementarity. RNAplex was then used to select trans-acting target gene with threshold parameter of duplex energy–e set as 20.

KEGG pathway analysis was performed for the predicted cis- or trans-target genes among differentially expressed genes to identify important signaling pathways. The correlations of differentially expressed lncRNAs with predicted cis- or trans-target genes were evaluated by calculating the Pearson correlation coefficient R values. The lncRNA-target gene correlations with the |R| > 0.85 and FDR < 0.0005, and the top 30 KEGG pathways were used to construct lncRNA-target pathway network.

Validation of lncRNA and mRNA expressions

The expressions of lncRNAs and mRNAs that had high number of degree in lncRNA-target-pathway network were further confirmed by qRT-PCR. Lung PET-high and PET-low tissues were collected from 11 patients. Total RNA was isolated from each sample using TRIzol reagent and treated with DNase I (Invitrogen) according to the manufacturer’s protocol. RNA was reverse transcribed into cDNA using the PrimeScript RT reagent Kit (Takara, Shiga, Japan) and amplified using the SYBR Green Kit (Takara, Shiga, Japan) on an ABI Prism 7300 sequence detection system (Applied Biosystems, Foster City, USA). The reaction mixture contained cDNA, 200 nM of each primer and 10 μl of 2 × SYBR Green PCR Master Mix (Takara, Shiga, Japan). The primer sequences of lncRNAs and mRNAs were listed in Supplementary Table 3. Reactions were incubated at 95° C for 30s, followed by 40 cycles of 95° C for 10 s and 60° C for 30 s. Reactions were run in triplicate for analysis. At the end of the PCR cycles, melting curve analysis was performed to verify product specificity. Expression levels of lncRNAs and mRNAs were normalized to housekeeping gene GAPDH and the 2^-ΔΔCt method was applied to calculate the relative expression levels.

Prediction efficiencies and bioinformatics analysis of lncRNA-binding proteins

The lncRNAs validated by qRT-PCR were further remained as indicators for lesion with high metabolic activity. The prediction efficiencies of lncRNAs were evaluated by receiver operating characteristic (ROC) curve analysis via pROC package (version 1.14.0) to calculate the area under the curve (AUC) [40].

Genomic locations of lncRNA-binding proteins were predicted by using UCSC Genome Browser (http://genome-asia.ucsc.edu/index.html). Secondary structures were identified via RNAfold minimum free energy estimations based on RNAfold webserver (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). TRANSFAC (http://www.gene-regulation.com/index2.html) was used to predict the potential transcription factors (TFs) of lncRNAs. Moreover, catRAPID analysis (http://service.tartaglialab.com/) was conducted to predict the potential interacting proteins of lncRNAs.

Ethical approval

The Permission of the Hospital Ethics Committee was obtained by the Shanghai public health clinical center's ethics committee before operation for each patient. (Approval number:2019-S009-02).