antibiotics Article Weighted Gene Co-Expression Network Analysis Identifies Key Modules and Hub Genes Associated with Mycobacterial Infection of Human Macrophages Lu Lu 1,2, *,† , RanLei Wei 3,† , Sanjib Bhakta 4 , Simon J. Waddell 5 and Ester Boix 2, * 1 College of Animal Science and Technology, Sichuan Agricultural University, Chengdu 610000, China 2 Department of Biochemistry and Molecular Biology, Faculty of Biosciences, Universitat Autonoma de Barcelona, 08290 Cerdanyola del Vallès, Spain 3 Laboratory of Omics Technology and Bioinformatics, West China Hospital, Sichuan University, Chengdu 610017, China; weijinrl@163.com 4 Mycobacteria Research Laboratory, Department of Biological Sciences, Institute of Structural and Molecular Biology, Birkbeck, University of London, London WC1E 7HX, UK; s.bhakta@bbk.ac.uk 5 Global Health and Infection, Brighton and Sussex Medical School, University of Sussex, Brighton BN1 9PX, UK; S.Waddell@bsms.ac.uk * Correspondence: lu.lu@sicau.edu.cn (L.L.); Ester.Boix@uab.cat (E.B.) † These authors have contributed equally to this work. Abstract: Tuberculosis (TB) is still a leading cause of death worldwide. Treatments remain unsatis- factory due to an incomplete understanding of the underlying host–pathogen interactions during infection. In the present study, weighted gene co-expression network analysis (WGCNA) was con-  ducted to identify key macrophage modules and hub genes associated with mycobacterial infection.  WGCNA was performed combining our own transcriptomic results using Mycobacterium aurum- Citation: Lu, L.; Wei, R.; Bhakta, S.; infected human monocytic macrophages (THP1) with publicly accessible datasets obtained from Waddell, S.J.; Boix, E. Weighted Gene three types of macrophages infected with seven different mycobacterial strains in various one-to-one Co-Expression Network Analysis combinations. A hierarchical clustering tree of 11,533 genes was built from 198 samples, and 47 dis- Identifies Key Modules and Hub tinct modules were revealed. We identified a module, consisting of 226 genes, which represented the Genes Associated with Mycobacterial Infection of Human Macrophages. common response of host macrophages to different mycobacterial infections that showed significant Antibiotics 2021, 10, 97. enrichment in innate immune stimulation, bacterial pattern recognition, and leukocyte chemotaxis. https://doi.org/10.3390/ Moreover, by network analysis applied to the 74 genes with the best correlation with mycobacteria antibiotics10020097 infection, we identified the top 10 hub-connecting genes: NAMPT, IRAK2, SOCS3, PTGS2, CCL20, IL1B, ZC3H12A, ABTB2, GFPT2, and ELOVL7. Interestingly, apart from the well-known Toll-like Academic Editor: John E. Gustafson receptor and inflammation-associated genes, other genes may serve as novel TB diagnosis markers Received: 28 November 2020 and potential therapeutic targets. Accepted: 14 January 2021 Published: 20 January 2021 Keywords: mycobacterium; macrophage; WGCNA; transcriptome; network analysis Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affil- 1. Introduction iations. Tuberculosis (TB) is still a global threat and one of the leading infectious diseases according to the World Health Organisation, causing 1.5 million deaths and 10 million new cases in 2018 [1]. TB infection begins when the causal bacilli, Mycobacterium tuberculosis (Mtb), reach the alveolar air sacs of the lungs, where they invade and replicate within Copyright: © 2021 by the authors. macrophages [2,3]. Granulomas are formed when macrophages, T lymphocytes, B lympho- Licensee MDPI, Basel, Switzerland. cytes, and fibroblasts aggregate, surrounding the infected macrophages. The granuloma This article is an open access article prevents dissemination of the mycobacteria but provides a local environment for the distributed under the terms and conditions of the Creative Commons infection of immune cells [3]. The formation of granulomas may also suppress host im- Attribution (CC BY) license (https:// mune responses, as macrophages and dendritic cells in the granulomas are unable to creativecommons.org/licenses/by/ present antigen to lymphocytes [4]. In addition to active disease, Mtb may also cause 4.0/). latent infection [5]. Epidemiology studies indicate that 90% of those infected with TB are Antibiotics 2021, 10, 97. https://doi.org/10.3390/antibiotics10020097 https://www.mdpi.com/journal/antibiotics Antibiotics 2021, 10, 97 2 of 16 asymptomatic, with only 10% lifetime risk to develop active disease [6]. Unfortunately, the risk of developing active TB increases dramatically in HIV patients [7]. The emergence of a growing number of multi-drug resistant and extensively drug-resistant TB cases makes the global situation much worse and emphasises the need for the development of new anti-TB agents able to eradicate tuberculosis [8,9]. It is also important to note that different mycobacteria induce distinct host responses, and this may contribute to the highly variable outcomes of infection ranging from asymptomatic to severe extra-pulmonary disease [10]. The impact of Mtb strain variation on human disease is not well established, and Mtb genomic diversity remains to be unambiguously correlated to the intracellular/lung lesion phenotypic diversity [11]. Understanding this host–pathogen relationship is fundamental to the development of new transformative therapies to prevent TB. Transcriptomic profiling, the analysis of genome-wide RNA expression, is a com- mon and powerful approach to investigate host and pathogen processes in infection models. Transcriptomic data of macrophages infected with various mycobacteria, includ- ing M. tuberculosis H37Rv (Mtb-H37Rv) [12,13], M. tuberculosis GC1237 (Mtb-GC1237) [13], M. smegmatis [13], M. bovis BCG (BCG) [13–15], and M. abscessus (MAB) [16] exist in public databases; however, the comparative response behaviours of these human macrophages to these mycobacterial species has not been addressed. Moreover, M. smegmatis, BCG, M. abscessus as well as M. aurum, are frequently used as surrogate experimental models for Mtb, as Mtb-H37Rv is a slow-growing and pathogenetic species. In particular, M. aurum has been shown to be a good surrogate model that might be used in high-throughput screening for tuberculosis drug discovery [17]. However, we do not know to what extent the host macrophage response to M. aurum mimics the host response to H37Rv infection. Weighted gene co-expression network (WGCNA) analysis is a systems biology method for describing correlation patterns among genes in genome-wide RNA profiling experi- ments, from microarray samples or RNAseq samples, which has been successfully used to identify highly connected hub genes and significant modules [18]. WGCNA alleviates the multiple testing problem inherent in microarray data analysis. Instead of relating thousands of genes to a microarray sample trait, it focusses on the relationship between a few (typically less than 10) modules and the sample trait. It also enables datasets from different platforms to be compared together by focussing on groups of genes rather than on the expression of individual genes [18]. Moreover, cut-off criteria are not necessary for WGCNA, and hence, important information omitted by other statistical testing methods can be retrieved. To address the above-mentioned questions, we first captured the expres- sion profiles of macrophages infected with M. aurum and then employed an integrative network-based approach to identify host response co-expression networks in other my- cobacterial infection models: Mtb-H37Rv, BCG, MAB-Rough, MAB-Smooth, M. smegmatis, Mtb-GC1237, and heat-inactivated Mtb-H37Rv. The integration of differential gene expres- sion, data-driven correlation, and co-expression networks enabled the reconstruction of novel signalling maps underlying the host response to mycobacterial infection. To avoid analytical bias from the individual experiments, we selected as a starting analysis material the raw mRNA sequencing data available at NCBI Sequence Read Archive or at EMBL-EBI European Nucleotide Archive. By applying the WGCNA methodology to 198 transcriptome profiles, we identified a defined set of core genes consisting of 226 genes that were commonly regulated dur- ing host macrophage cell infection with the majority of mycobacterial species. These genes are involved in immune response stimulation by mycobacteria infection in human macrophages, which are mainly linked to Toll-like receptor and TNF signalling pathways. Of them, 74 genes were identified as hub genes representing the transcriptomic response of macrophages to mycobacterial infection. Moreover, we identified a set of 10 key hub genes that may serve to outline a common infection profile. Antibiotics 2021, 10, 97 3 of 16 2. Results 2.1. Construction of Weighted Co-Expression Network To learn about the immune response to infection by different mycobacteria, we inves- tigated the in vitro gene regulatory response of macrophages to infection using multiple Mtb strains and related mycobacterial species (Table 1). Table 1. Human macrophage gene expression profiles by RNAseq. Bioproject Host Cell Mycobacterium MOI POI/h Ensemble ID Libraries Reference 2 PRJNA295153 1 MAB 10 1 h, 4 h, 24 h 21657 27 [16] THP1 BMDM M1 PRJNA355844 BCG 0.18 24 h 25343 15 [19] BMDM M2 HAM PRJNA471095 Mtb-H37Rv 2 2 h, 24 h, 72 h 20774 36 [20] BMDM PRJNA495462 BMDM Mtb-H37Rv 0.5 18 h 39354 6 [21] PRJNA575195 THP1 M. aurum 10 24 h 58233 6 Mtb-H37Rv H37Rv + PRJNA279959 BMDM M. smegmatis 2 4 h, 18 h, 48 h 12728 108 [13] Mtb-GC1237 BCG 1Abbreviation used: THP1, THP1 activated macrophage; BMDM, blood monocyte-derived macrophage, including type 1 (indicated as BMDM-M1) and type 2 (indicated as BMDM-M2); HAM, human alveolar macrophage. Mtb- H37Rv, Mycobacterium tuberculosis H37Rv; MAB, Mycobacterium abscessus; BCG, M. bovis Bacillus Calmette–Guérin (BCG); M. aurum, Mycobacterium aurum; H37Rv+, heat-inactivated Mtb-H37Rv; M. smegmatis, Mycobacterium smegmatis; Mtb-GC1237, Mycobacterium tuberculosis GC1237; 2 including two cell lines, MAB-smooth and MAB- rough. 2 See Table S1 for a list of all accession numbers corresponding to each Bioproject. Specifically, we collected RNAseq data of human macrophage infected and/or stim- ulated with the following mycobacteria: Mtb H37Rv (pathogenic and slow growing), the common strain often used in laboratory experiments; Mtb GC1237, which is a clinical isolate of the virulent Beijing family; Mycobacterium smegmatis, which is fast-growing and non-pathogenic; M. bovis bacillus Calmette-Guerin (BCG), which is slow growing and an attenuated Mycobacterium bovis strain used for vaccinations; MAB-R and MAB-S, which are the rough and smooth morphotypes of the fast-growing M. abscessus; and heat-killed Mtb H37Rv. These previously published transcriptome studies were combined with our own data using M. aurum (PRJNA575195). The detailed features of mycobacteria are listed in Table 2. Table 2. Characteristics of studied mycobacteria. Strain Growth Rate Pathogenicity Morphology Note Mycobacterium tuberculosis Slow growth Tuberculosis Rough Laboratory strain Mtb H37Rv Mycobacterium tuberculosis Slow growth Tuberculosis Rough GC1237 M. bovis Bacillus An attenuated strain of M. Slow growth Non-pathogenic Rough Calmette–Guérin (BCG) bovis, used as a vaccine for TB M. smegmatis Fast growth Non-pathogenic Rough M. abscessus Rough (ABR) Fast growth Pathogenic Rough M. abscessus Smooth (ABS) Fast growth Pathogenic Smooth M. aurum Fast growth Non-pathogenic Smooth To construct the gene co-expression networks, all data were downloaded from the Gene Expression Omnibus (GEO) database, including our own unpublished experimental data, as detailed in Table 1 and Table S1. Raw data were preprocessed identically using R for background correction and normalisation. A total of 11,533 genes (Appendix A (Additional File 1)) for 198 samples were used to build a WGCNA. First, genes and conditions were clustered using the flashClust function (see Figure S1, where information on infection status, infection time, or infection type/different bacterial strains is indicated). Selection of the soft-thresholding power is an important step when constructing a WGCNA [22]. We performed the analysis of network topology for thresholding powers from 1 to 50 and identified the relatively balanced scale independence and mean connectivity of the WGCNA. As shown in Figure 1, power value 30, which was the lowest power for the data, as detailed in Table 1 and Table S1. Raw data were preprocessed identically using R for background correction and normalisation. A total of 11,533 genes (Appendix A (Ad- ditional File 1)) for 198 samples were used to build a WGCNA. First, genes and conditions were clustered using the flashClust function (see Figure S1, where information on infec- Antibiotics 2021, 10, 97 tion status, infection time, or infection type/different bacterial strains is indicated). Selec- 4 of 16 tion of the soft-thresholding power is an important step when constructing a WGCNA [22]. We performed the analysis of network topology for thresholding powers from 1 to 50 and identified the relatively balanced scale independence and mean connectivity of the WGCNA. As shown in Figure 1, power value 30, which was the lowest power for the scale-free topology fittopology scale-free index fitonindex 0.85, on was selected 0.85, was selectedto to produce a hierarchical produce a hierarchical clustering clustering tree tree (dendrogram) of 11,533 (dendrogram) genes. of 11,533 From genes. FromthetheWGCNA analysis, WGCNA analysis, we identified we identified 47mod- 47 distinct distinct modules withules with mergeCutHeight mergeCutHeight of 0.3 of 0.3 (Figure2). (Figure 2). Figure 1.Figure 1. Determination Determination of soft-thresholding of soft-thresholding power power for the weighted forco-expression gene the weighted geneanalysis network co-expression (WGCNA) network analysis.analysis (a) Analysis of the scale-free (WGCNA) fit index (a) analysis. for various Analysissoft-thresholding powersfit of the scale-free (β);index (b) Analysis of the mean for various connec- soft-thresholding tivity for various soft-thresholding powers. In all, 30 was the most fit power value. powers Antibiotics 2021, 10, x FOR (β); (b) Analysis of the mean connectivity for various soft-thresholding powers.5 In PEER REVIEW all, 30 was of 16 the most fit power value. Figure of Figure 2. Construction Construction 2. co-expression of co-expression modules modules by by WGCNA of macrophage WGCNAresponses transcriptional of macrophage transcriptional re- to mycobacterial infection. The cluster sponses dendrogram of genes. Each to mycobacterial branch in the infection. Thefigure represents cluster one gene, and dendrogram ofevery colour genes. below Each repre- in the figure branch sents one co-expression module. represents one gene, and every colour below represents one co-expression module. 2.2. Correlation between Modules and Identification of Key Modules In WGCNA, modules are defined as clusters of highly interconnected genes; genes within the same cluster have high correlation coefficients. The module eigengene (ME) is representative of the corresponding module’s gene expression profile correlated with a defined trait (infection status). We analysed the interaction relationships of the 47 modules and plotted the network heatmap (Figure S2). The results showed that each module was Antibiotics 2021, 10, 97 5 of 16 2.2. Correlation between Modules and Identification of Key Modules In WGCNA, modules are defined as clusters of highly interconnected genes; genes within the same cluster have high correlation coefficients. The module eigengene (ME) is representative of the corresponding module’s gene expression profile correlated with a defined trait (infection status). We analysed the interaction relationships of the 47 modules and plotted the network heatmap (Figure S2). The results showed that each module was an independent validation from the others, which indicated a high level of independence among the modules and relative independence of gene expression within each module. Moreover, we calculated eigengenes and clustered them according to their correlations. The results demonstrated that the 47 modules were mainly divided into two clusters (Figure S3). As shown in Figure 3, we found that several modules were significantly correlated with the infection status (significance value < 0.001), which was defined as presence or absence of mycobacterial infection. Among them, the grey60 module showed the highest correlation (r value) compared with the other modules. Therefore, grey60 was considered as the most relevant differentially expressed gene module across these diverse mycobacteria and variety of host macrophage cells exposed to infection. Figure 4 illustrates the correlation between module membership and gene significance in the grey60 module.6 of 16 Antibiotics 2021, 10, x FOR PEER REVIEW The results indicated that genes in the grey60 module likely play an important role in responding to mycobacterial infection in macrophages. Figure 3. Heatmap Figure of the 3. Heatmap correlation of the correlationbetween moduleeigengenes between module eigengenes andand mycobacterial mycobacterial infection infection status. status. Theofcolour The colour of each each cell cell atatthe the row–column intersection indicates the correlation coefficient between the gene module and mycobacterial infection. in- row–column intersection indicates the correlation coefficient between the gene module and mycobacterial fection. TheThe number number of genes of genes in each in each modulemodule is indicated is indicated on the on left.the left. The Theofvalue value of the correlation the correlation coefficientcoefficient and signifi- and significance cancebetween between each module and infection status are indicated. each module and infection status are indicated. Figure 3. Heatmap of the correlation between module eigengenes and mycobacterial infection status. The colour of each cell at Antibiotics the10, 2021, row–column 97 intersection indicates the correlation coefficient between the gene module and mycobacterial in-6 of 16 fection. The number of genes in each module is indicated on the left. The value of the correlation coefficient and signifi- cance between each module and infection status are indicated. Figure Figure 4. Scatter plot of module 4. Scatterin eigengenes plot theofgrey60 module eigengenes module in the that was grey60 module significantly that was associated withsignificantly macrophage associated re- with macrophage response to mycobacterial infection. sponse to mycobacterial infection. 2.3. Functional Enrichment and Identification of Hub Genes 2.3. Functional Enrichment and Identification of Hub Genes Following this observation, we inspected in detail the composition of the grey60 Following this observation, we inspected in detail the composition of the grey60 module. The module is correlated with mycobacterial infection status and contains 226 module. The module is correlated with mycobacterial infection status and contains 226 genes, the expression of each gene is listed in Appendix B (Additional File 2). Appling genes, the expression of each gene is listed in Appendix B (Additional File 2). Appling Kyoto Encyclopedia of Genes and Genome (KEGG) and Gene Ontology (GO) enrichment Kyoto Encyclopedia of Genes and Genome (KEGG) and Gene Ontology (GO) enrichment analysis, we found these 226 genes significantly enriched into pathways such as Cytokine- analysis, we found these 226 genes significantly enriched into pathways such as Cytokine- cytokine receptor interaction, mitogen-activated protein kinase (MAPK) signalling, toumor necrosis factor (TNF) signalling, phosphatidylinositol 3-kinase (PI3K)-protein kinase B (Akt) signalling, interleukin 17 (IL-17) signalling, and biological processes such as response to lipopolysaccharides, response to molecule of bacterial origin, and positive regulation of peptidyl–tyrosine phosphorylation, respectively (full list included in Appendix C (Addi- tional File 3)). A heatmap of the weighted relative expression of each gene in the grey60 module (Figure 5) revealed that many module genes were strongly expressed after mycobac- teria infection compared with the uninfected control groups. Among the three macrophage controls, blood monocyte-derived macrophage (BMDM) had the lowest expression of these genes, followed by the THP1 and human alveolar macrophage (HAM) cells. Among the different mycobacterial infections, Mtb GC1237 and M. smegmatis caused the highest upregulation of these genes and were followed by BCG, Mtb H37Rv, MABS, MABR, heat inactivated Mtb H37Rv, and M. aurum. macrophage controls, blood monocyte-derived macrophage (BMDM) had the lowe pression of these genes, followed by the THP1 and human alveolar macrophage (H cells. Among the different mycobacterial infections, Mtb GC1237 and M. smegmatis ca Antibiotics 2021, 10, 97 the highest upregulation of these genes and were followed by BCG, Mtb7 ofH37Rv, 16 M MABR, heat inactivated Mtb H37Rv, and M. aurum. Figure 5. Heatmap and eigengene expression profile for grey60 module, induced by macrophages after mycobacterial Figure 5. Heatmap and eigengene expression profile for grey60 module, induced by macrophages after mycobacteria infection. (a) The heat map shows the relative expression of each of the 226 genes; (b) Plot of eigengene values for the infection. (a)grey60 The module. heat map The shows the relative y-axis indicates expression the value of the moduleofeigengene; each of the thex-axis 226 indicates genes; (b) Plottype. sample of eigengene values for the RV, Mtb-H37Rv; grey60 module. BCG, The y-axis indicates M. bovis-BCG; MABR, M. the value of the abscessus-rough; MABS,module eigengene;Smeg, M.abscessus-smooth; the x-axis indicates M. smegmatis; sample type. GC, Mtb-GC1237; RV+, RV, heat Mtb-H37Rv BCG, M. bovis-BCG; inactivated MABR, Mtb-H37Rv;M.Aurum, abscessus-rough; MABS, M. aurum; Control, M.abscessus-smooth; uninfected Smeg, M. smegmatis; GC, Mtb-GC1237; RV+ control macrophages. heat inactivated Mtb-H37Rv; Aurum, M. aurum; Control, uninfected control macrophages. Out of the 226 genes from the grey60 module, 74 genes showed the best linear correla- tion with mycobacteria infection (Figure 4, see the full gene list in Appendix D (Additional Out File of 4)).the 226 Thus, genes those fromwere 74 genes theused grey60 module,for as candidates 74hub genes geneshowed the analy- identification best linear c lation sis with mycobacteria by applying infectionplatform the NetworkAnalyst (Figureto4,map seethe the full gene list protein-protein in Appendix interactions and D (A visualised them as networks in Cytoscape (Figure 6). The first 10 genes tional File 4)). Thus, those 74 genes were used as candidates for hub gene identific in the grey60 module according to their interaction degree, connected to most of the other genes within analysis by applying the NetworkAnalyst platform to map the protein-protein interac the module, were highlighted as main hub genes: nicotinamide phosphoribosyltrans- and visualised theminterleukin ferase (NAMPT), as networks in Cytoscape 1 receptor associated(Figure kinase 2 6). The first (IRAK2), 10 genes suppressor in the g of cy- module according tokine signallingto3 (SOCS3), their interaction degree, connected prostaglandin-endoperoxide to most synthase of theC-C 2 (PTGS2), other motifgenes w chemokine ligand 20 (CCL20), interleukin 1 beta (IL1B), zinc finger CCCH-type containing the module, were highlighted as main hub genes: nicotinamide phosphoribosyltran 12A (ZC3H12A), ankyrin repeat and BTB domain containing 2 (ABTB2), glutamine fructose ase (NAMPT), 6 phosphateinterleukin transaminase 12 (GFPT2), receptorandassociated kinase elongation of 2 (IRAK2), very long chain fattysuppressor acids protein of cyt signalling 3 (SOCS3), 7 (ELOVL7). We also prostaglandin-endoperoxide performed enrichment analysis to synthase explore the2biological (PTGS2), C-C motif ch processes (BP) and pathways in which the key module was involved using the 74 kine ligand 20 (CCL20), interleukin 1 beta (IL1B), zinc finger CCCH-type containing genes that showed better linear correlation within the grey60 module. GO enrichment of BP was conducted by (ZC3H12A), ankyrin repeat and BTB domain containing 2 (ABTB2), glutamine fruct R programme clusterProfiler, as detailed in Figure 7. The findings are consistent with the re- sults using all 226 genes in the grey60 module (Appendix C (Additional File 3)). The results highlight that genes in the grey60 module were enriched in innate immune stimulation responses, such as pathogen-associated molecular patterns (PAMPs), cytokine secretion, 7 (ELOVL7). We also performed enrichment analysis to explore the biological processes (BP) and pathways in which the key module was involved using the 74 genes that showed better linear correlation within the grey60 module. GO enrichment of BP was conducted by R programme clusterProfiler, as detailed in Figure 7. The findings are consistent with Antibiotics 2021, 10, 97 the results using all 226 genes in the grey60 module (Appendix C (Additional File 3)).8 of The16 results highlight that genes in the grey60 module were enriched in innate immune stimu- lation responses, such as pathogen-associated molecular patterns (PAMPs), cytokine se- cretion, and immune cell migration, all of which were positively correlated with myco- and immune cell migration, all of which were positively correlated with mycobacterial bacterial infection. infection. Figure 6. Figure Identification of 6. Identification of hub hub genes genes in in the the grey60 grey60 module module induced induced by by macrophages macrophages after after mycobacterial mycobacterial infection. infection. The The interaction network of of the the selected selected 74 74genes genesfrom fromthe the226-gene 226-genegrey60 grey60module builtbybyNetwork modulebuilt NetworkAnalyst. Analyst.Size Sizeofof the circle the is circle is relative relative to to number number of of thethe connections. connections. Antibiotics 2021, 10, x FOR PEER REVIEW 9 of 16 Antibiotics 2021, 10, 97 9 of 16 Figure 7. Biological process (GO) enrichment analysis of 74 genes in the grey60 module. The enriched biological processes Figure were 7. Biological coloured in red process (GO) enrichment or blue (according analysis to adjusted of 74 genes p-values) in the grey60 in comparison to themodule. TheGO connected enriched terms, biological processes which are shown in were coloured in red or blue (according to adjusted p-values) in comparison to the connected GO terms, which are shown grey. in grey. 3. 3. Discussion Discussion The The interaction interaction between between macrophage macrophage and and bacterium bacterium is is central to the immunopathol- ogy ogy of mycobacterial diseases; a deeper understanding of this of mycobacterial diseases; a deeper understanding of this interplay interplay will will identify identify new new treatment strategies. Here, we aimed to address the question of whether different my- cobacteria similarly modulate the host response of macrophages. Direct comparison of the Antibiotics 2021, 10, 97 10 of 16 differentially expressed genes (DEGs) of macrophages infected with different mycobacte- ria identified from different individual studies was not appropriate to meet this goal, as the analyses were performed differently by each group and, more fundamentally, DEG analysis neglects any key but minor changes in gene expression. For example, when we first explored the transcriptional change of THP1-derived macrophages challenged with M. aurum by DEG analysis, we only identified 27 differentially expressed genes in com- parison with uninfected control macrophages (Appendix E (Additional File 5)). The most up-regulated and down-regulated genes were TMEM189-UBE2V1 and RPS18, respectively. TMEM189-UBE2V1 is related to the IL-1 signalling pathway [23,24]; however, standard DEG analysis did not identify any further IL-1 related genes. To note, M. aurum shares a high similarity with Mtb at the genomic level [25], and most of the genes reported to be associated with drug resistance are common [26]. To compare between published macrophage transcriptome signatures and to over- come the DEG analysis drawback of neglecting minor changes in gene expression, we applied WGCNA to identify a core set of macrophage mycobacterial-responsive genes, contrasting the RNA profiles of macrophages responding to different mycobacterial infec- tions including seven mycobacterial exposures and three types of macrophages. These mycobacteria are known to differ in pathogenicity, growth rate, morphology, and cell wall lipid composition (Table 2). The selected mycobacterial species and macrophage types represent most of the working models that have been widely used for Mtb infection characterisation [3,15,27]. In our study, by WGCNA analysis of the 198 transcriptome samples, we can categorise 11,533 genes into 47 modules according to their expression pattern (Figure 2). Among the defined modules, we identified a module (grey60), consisting of 226 genes, which strongly correlates with mycobacterial adaptations to infection (Figure 4). The grey60 module is mostly enriched with genes from biological processes such as response to pathogen-associated molecules, cytokine secretion, and chemoattractant release (Figure 7). Specifically in the module, we find genes induced in response to infection such as IL1B, TNF, IL6, IL12B, NFKBIA, JUN, and MAP3K8, which are related to Toll-like receptor signalling (Appendix B (Additional File 2)). To note, IL1B, TNF, IL1A, IL6, IL23A, IL12B, PLK3, and IRAK2 were already reported as Mtb-responsive genes [27–32]. Next, by protein interaction network analysis, we identified a set of hub genes that represent the general response of macrophage infected with different mycobacteria, the top 10, connecting the genes in the module together, were NAMPT, IRAK2, SOCS3, PTGS2, CCL20, IL1B, ZC3H12A, ABTB2, GFPT2, and ELOVL7. Of them, IRAK2 and IL1B are well known to be related to Toll-like receptor signalling pathways [33]; SOCS3, PTGS2, CCL20, and IL1B are linked to TNF signalling pathways [34]; ZC3H12A and PTGS2 are involved in cellular inflammatory responses [35]; and GFPT2, NAMPT, ABTB2, and ELOVL7 are related to cell cycle dysregulation, tissue damage, and autophagy [36]. IFN-γ, IP-10, CRP, TNF-α, CCL4, IL1β, and TLR4 have been associated with high accuracy to TB [37,38]. Genes involved in NAD+ biosynthesis, for example, NAMPT, were strongly upregulated during Mtb infection [39]. The current data emphasise the previously suggested key role of the IFN pathway in the macrophage response to Mtb and successful outcome of infection [39,40]. Interestingly, some of these genes are involved in nucleic acid or lipid turnover and are also activated during other diseases, such as autoimmunity or cancer [36,40–45]. Our results back up very recent studies of Mtb modulation of host cell metabolism pathways [39]. Moreover, we note here the role of cell migration factors in TB pathogenesis, besides cell response to bacteria, as the genes in grey60 split into the two directions (Figure 7). In addition, interferon-γ release assay (IGRA), based on IFN-γ response, is widely used in TB diagnosis [46]. Except for the well-known TB markers such as IRAK2, SOCS3, PTGS2, CCL20, and IL1B, the other hub genes we identified here, NAMPT, ZC3H12A, ABTB2, GFPT2, and ELOVL7, may be useful as new biomarkers for TB diagnosis and as therapeutic targets for host-directed strategies. The present results based on the host metabolism can provide a complementary approach to other studies strictly focused on the mycobacterial Antibiotics 2021, 10, 97 11 of 16 response to available drugs in the market [47–51]. Hopefully, the dual sequencing of host–pathogen models [52] will help to unravel the key genes associated to both virulence and host response to infection and set the path to a better tailored drug design. We conclude here that different mycobacterial infection models can trigger a common core set of genes related to macrophage activation and migration, which is separate from macrophage stimulation with heat-killed Mtb. However, we cannot disregard the inherent differences between mycobacteria strains and macrophage type, as shown in Figure 6. Control uninfected BMDM had the lowest expression of this gene module, HAM had the highest level of expression of these genes, whereas THP1 showed a medium profile. The three studied macrophage types include two primary macrophage cells: blood monocyte- derived macrophages and alveolar macrophages, together with an immortalised monocyte- like cell line (THP1). Alveolar macrophages are natural mature macrophages resident in the lung and are among the first immune cells to encounter invading mycobacteria, while blood monocytes and THP1 cells need to be activated to become macrophages. We cannot disregard the key differences between immortalized THP1 cells and their considered physiological counterparts. For example, THP-1 cells, when compared with monocytes, are far less responsive to some pathogen recognition patterns [53]. Therefore, our results emphasise once more the importance of selecting the appropriate cell infection model to suit each study. Comparing the different mycobacteria strains, Mtb H37Rv, Mtb GC1237, M. bovis BCG, M. smegmatis, MABS, and MABR showed similar induction of genes in the grey60 module, while infection with M. aurum or stimulation with heat-inactivated Mtb H37Rv resulted in a lower magnitude of upregulation of these genes. The two species, M. aurum and M. smegmatis, have been regarded as ideal surrogates for mimicking tuberculosis in drug development, both demanding a shorter time of culture and less stringent biosafety laboratory conditions [25,54].The differential regulation of the gene module identified here by all mycobacteria highlights the potentiality of some alternative infection models to live Mtb infection. 4. Materials and Methods 4.1. RNA Sequencing of Human Monocyte-Derived Macrophage (THP-1) Infected with M. Aurum M. aurum (NCTC 10437) infection of human monocytic macrophage cell line THP-1 (NCTC #88081201) were performed as previously described; triplicate biological repli- cates were performed [55]. After 24 h infection, uninfected and infected macrophage cells were washed 3 times with prewarmed PBS and collected by cell scrapper for RNA extraction. Total RNA was extracted using mirVanaTM miRNA Isolation Kit as described by the manufacturer (Ambion, Life Technologies, AM1560). RNA purity was determined by spectrophotometry, and RNA integrity was analysed using Agilent 2100 Bioanalyzer and calculated as an RNA integrity number. Following RNA extraction, RNA sequenc- ing libraries were prepared according to protocols provided by Illumina. 50 bp-long single-end sequencing was carried out in an Illumina HiSeq2000 sequencer with a depth of >20 million reads per sample at the CRG genomics (Centre for Genomic Regulation, Barcelona) facility. The raw transcriptomic data were uploaded in NCBI Sequence Read Archive: PRJNA575195. After sequencing, quality assessment of reads was carried out us- ing FastQC [56] to assess the distribution of Phred quality scores and mean percentage GC content across each read. Adapters were trimmed using Trimmomatic [57]. Trimmed reads were aligned to the latest human genome assembly (version GRCh38) using Hisat2 [58]. Aligned reads were converted and sorted in the SAM file format using Samtools, providing the sorted BAM file. Stringtie was used to count the number of reads mapping to the GTF file of the annotated genome using BAM file as input [59]. Following, DEseq2 was used to reveal differently expressed genes (DEGs) using the read count as input file, and the results of all statistical tests were adjusted for multiple testing (0.05) [60]. Antibiotics 2021, 10, 97 12 of 16 4.2. Data Collection and Sample Processing for WGCNA Mycobacterial macrophage infection gene expression studies were reviewed and fil- tered based on the inclusion–exclusion criteria of the following: (1) human macrophage; (2) at least 2 h post infection time points; (3) minimum sample size of 3 control/3 exper- imental. After filtering, 6 genome-wide expression studies were selected, including our own transcriptome of M. aurum-infected macrophages from this study, PRJNA575195. As shown in Table 1, these datasets involved three different macrophage cell types (BMDM, THP1, HAM) infected with 8 different mycobacteria (Mtb-H37Rv, BCG, MAB-S, MAB-R, M. smegmatis, Mtb-GC1237, M. aurum, and heat-inactivated Mtb-H37Rv). All data were generated using a whole genome platform (RNAseq). In total, 198 samples were analysed, as summarised in Table 1 and detailed in full in Table S1. Ensembl IDs from RNAseq were converted into gene symbols and merged gene expression into two groups. For sequencing data, genes with less than 10 average counts per million were filtered out of the analysis. Then, the expression data were transformed using log2 transformation to normalise the data [61]. The well-established ComBat procedures in SVA packages were performed to reduce the potential study-specific batch effect [62]. Four samples were detected as outliers and filtered out using flashClust [63] (cluster dendrogram detailed in Figure S1). 4.3. Construction of Weighted Gene Co-Expressed Networks We performed WGCNA to identify the gene modules of interest from the assembled dataset using R package WGCNA [18]. The standard WGCNA procedure generates a squared adjacency matrix, between genes, based on their correlation [64]. In the present study, the absolute value of the Pearson correlation between the expression profiles of all candidate genes was determined for the 11532 most varying non-redundant transcripts, which was transformed into a connection strength measure by using a power function (connection strength (i,j) = |correlation(i,j)|ˆ β). The scale-free topology criterion was applied to determine the lowest soft threshold power guaranteeing a scale-free topology fit (R2 > 0.85) [18]. To group genes with coherent expression profiles into modules, we used the WGCNA R packages average lineage hierarchical clustering to measure the dissimilarity by topological overlapping. The co-expression modules were constructed using the automatic network construction function blockwiseModules with the following settings: power, 30; degree of similarity, 0.75; minModuleSize, 30; and TOMType, signed. All other parameters of the modules were set to default [18]. Gene significance (GS) was defined as the log10 transformation of P-value (GS = lgP) in the linear regression between gene expression and sample information. In addition, module significance (MS) was defined as the average GS for all the genes in a module. Modules with significance less than 0.0001 and an absolute value of correlation coefficient higher than 0.5 were recognised as significantly associated modules. 4.4. Functional Enrichment of Recurrence-Associated Modules The significantly associated modules were applied to GO and KEGG databases for function and pathway enrichment analyses [65]. P-Value corrected by Benjamin–Hochberg FDR less than 0.05 was the threshold used to discriminate significant terms. 4.5. Identification of Hub Genes in Key Modules Hub genes that are highly interconnected with module nodes were considered as functionally significant. A protein–protein interaction network was constructed using NetworkAnalyzer [66] for the significantly associated module. The intramodular connectivity was calculated to identify hub genes (qvalue < 0.01). 5. Conclusions In this work, we compared the transcriptomic responses of different types of macrophages infected with seven different mycobacterial strains: four that are pathogenic and survive in- tracellularly and three that are non-pathogenic surrogates. By weighted gene co-expression Antibiotics 2021, 10, 97 13 of 16 network analysis (WGCNA) of 198 sample profiles, we identified a specific 226-gene mod- ule enriched with genes related to bacterial unique pattern recognition, cytokine secretion, and leukocyte chemotaxis. We also defined by protein–protein network analysis a core set of genes that may be used as novel diagnosis markers and therapeutic targets. We also observed that the magnitude of expression of this immune cell response signature dif- fered between the mycobacteria evaluated. Results highlight the close connection between the host immune system and pathogenicity, and suggest new targets for host-directed therapies. Supplementary Materials: The following are available online at https://www.mdpi.com/2079-638 2/10/2/97/s1, Table S1. Data information. Figure S1. Dendrogram of all samples. The leaves of the tree correspond to each sample, clustering expression ratios from 11,533 genes. The first colour band underneath the tree indicates which samples appear to be outlying (the outliers, coloured in red, were excluded from the following analysis), the second colour band represents the infection time. Similarly, the remaining bands mark the bacterial trait (different mycobacteria), cell (different macrophage type), infection status (infected or control, infected coloured in red), and each combination of macrophage cell type and mycobacterial specie. RV, Mtb-H37Rv; BCG, M. bovis-BCG; MABR, M. abscessus-rough; MABS, M.abscessus-smooth; Smeg, M. smegmatis; GC, Mtb-GC1237; RV+, heat inactivated Mtb-H37Rv; Aurum, M. aurum. Figure S2. Interaction relationship analysis of co-expressed genes. Different colours of horizontal axis and vertical axis represent different gene modules. The brightness of yellow to red of the intersects represent the degree of connectivity between different modules. There was no significant difference in interactions among different modules, indicating a high-scale independence degree among these modules. Figure S3. Cluster dendrogram of module eigengenes. Author Contributions: L.L. and E.B. designed the experimental work; L.L. performed the laboratory experimental work; L.L. and R.W. analysed the data; L.L., E.B. and R.W. prepared the first draft and edited the manuscript; L.L., E.B., R.W., S.B. and S.J.W. revised the final paper version. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the Ministerio de Economía y Competitividad, (MINECO): SAF2015-66007P and PID2019-106123GB-I00 co-financed by FEDER funds and Fundació La Marató de TV3 (ref. 20180310) to E.B., and the National Key R & D Program of China (2020YFA0509500). S.B. is CIPLA Distinguished Fellow in Pharmaceutical Sciences. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The transcriptomic data of M. aurum infected macrophage presented in this study are openly available in NCBI Sequence Read Archive: PRJNA575195. The other data PRJNA295153, PRJNA355844, PRJNA471095, PRJNA495462, and PRJNA279959 were downloaded from NCBI SRA. Conflicts of Interest: The authors declare that they have no competing interests. Appendix A Additional file 1. Input dataset for WGCNA. The file contains the batch-corrected log2 counts per million for the 11,533 Ensembl genes analysed in this study for each of the 198 samples. Appendix B Additional file 2. Excel of relative expression of genes in grey60 module. Appendix C Additional file 3. KEGG and GO enrichment of genes in grey60 module. Appendix D Additional file 4. Gene list of 74 candidate hub genes in the grey60 module. Antibiotics 2021, 10, 97 14 of 16 Appendix E Additional file 5. Differentially expressed genes triggered by M. aurum infection in THP1 macrophages. References 1. WHO. Global Tuberculosis Report 2019; WHO: Geneva, Switzerland, 2019; ISBN 9789241565714. 2. Gupta, A.; Kaul, A.; Tsolaki, A.G.; Kishore, U.; Bhakta, S. Mycobacterium tuberculosis: Immune evasion, latency and reactivation. Immunobiology 2012, 217, 363–374. [CrossRef] [PubMed] 3. Ferluga, J.; Yasmin, H.; Al-Ahdal, M.N.; Bhakta, S.; Kishore, U. Natural and trained innate immunity against Mycobacterium tuberculosis. Immunobiology 2020, 225, 151951. [CrossRef] [PubMed] 4. Tufariello, J.M.; Chan, J.; Flynn, J.L. Latent tuberculosis: Mechanisms of host and bacillus that contribute to persistent infection. Lancet Infect. Dis. 2003, 3, 578–590. [CrossRef] 5. Dey, B.; Bishai, W.R. Crosstalk between Mycobacterium tuberculosis and the host cell. Semin. Immunol. 2014, 26, 486–496. [CrossRef] [PubMed] 6. Connell, D.W.; Berry, M.; Cooke, G.; Kon, O.M. Update on tuberculosis: TB in the early 21st century. Eur. Respir. Rev. 2011, 20, 71–84. [CrossRef] 7. Corbett, E.L.; Watt, C.J.; Walker, N.; Maher, D.; Williams, B.G.; Raviglione, M.C.; Dye, C. The growing burden of tuberculosis: Global trends and interactions with the HIV epidemic. Arch. Intern. Med. 2003, 163, 1009–1021. [CrossRef] 8. Lohrasbi, V.; Talebi, M.; Bialvaei, A.Z.; Fattorini, L.; Drancourt, M.; Heidary, M.; Darban-Sarokhalil, D. Trends in the discovery of new drugs for Mycobacterium tuberculosis therapy with a glance at resistance. Tuberculosis 2018, 109, 17–27. [CrossRef] 9. Lu, L.; Li, J.; Moussaoui, M.; Boix, E. Immune Modulation by Human Secreted RNases at the Extracellular Space. Front. Immunol. 2018, 9, 1012. [CrossRef] 10. Coscolla, M.; Gagneux, S. Does M. tuberculosis genomic diversity explain disease diversity? Drug Discov. Today Dis. Mech. 2010, 7, e43–e59. [CrossRef] 11. Correa-Macedo, W.; Cambri, G.; Schurr, E. The Interplay of Human and Mycobacterium Tuberculosis Genomic Variability. Front. Genet. 2019, 10, 865. [CrossRef] 12. Jhingan, G.D.; Kumari, S.; Jamwal, S.V.; Kalam, H.; Arora, D.; Jain, N.; Kumaar, L.K.; Samal, A.; Rao, K.V.S.; Kumar, D.; et al. Comparative Proteomic Analyses of Avirulent, Virulent, and Clinical Strains of Mycobacterium tuberculosis Identify Strain-specific Patterns. J. Biol. Chem. 2016, 291, 14257–14273. [CrossRef] [PubMed] 13. Blischak, J.D.; Tailleux, L.; Mitrano, A.; Barreiro, L.B.; Gilad, Y. Mycobacterial infection induces a specific human innate immune response. Sci. Rep. 2015, 5, 1–16. [CrossRef] [PubMed] 14. Nalpas, N.C.; de Park, S.; Magee, D.A.; Taraktsoglou, M.; Browne, J.A.; Conlon, K.M.; Rue-Albrecht, K.; Killick, K.E.; Hokamp, K.; Lohan, A.J.; et al. Whole-transcriptome, high-throughput RNA sequence analysis of the bovine macrophage response to Mycobacterium bovis infection in vitro. BMC Genom. 2013, 14, 230. [CrossRef] [PubMed] 15. Malone, K.M.; Rue-Albrecht, K.; Magee, D.A.; Conlon, K.; Schubert, O.T.; Nalpas, N.C.; Browne, J.A.; Smyth, A.; Gormley, E.; Aebersold, R.; et al. Comparative ’omics analyses differentiate Mycobacterium tuberculosis and Mycobacterium bovis and reveal distinct macrophage responses to infection with the human and bovine tubercle bacilli. Microb. Genom. 2018, 4, 54. [CrossRef] [PubMed] 16. Aulicino, A.; Dinan, A.M.; Miranda-CasoLuengo, A.A.; Browne, J.A.; Rue-Albrecht, K.; MacHugh, D.E.; Loftus, B.J. High- throughput transcriptomics reveals common and strain-specific responses of human macrophages to infection with Mycobac- terium abscessus Smooth and Rough variants. BMC Genom. 2015, 16, 1046. [CrossRef] [PubMed] 17. Gupta, A.; Bhakta, S. An integrated surrogate model for screening of drugs against mycobacterium tuberculosis. J. Antimicrob. Chemother. 2012, 67, 1380–1391. [CrossRef] 18. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 2008, 9, 559. [CrossRef] 19. Lardone, R.D.; Chan, A.A.; Lee, A.F.; Foshag, L.J.; Faries, M.B.; Sieling, P.A.; Lee, D.J. Mycobacterium bovis Bacillus Calmette– Guérin Alters melanoma microenvironment Favoring antitumor T cell responses and improving M2 Macrophage Function. Front. Immunol. 2017, 8, 1–14. [CrossRef] 20. Papp, A.C.; Azad, A.K.; Pietrzak, M.; Williams, A.; Handelman, S.K.; Igo, R.P.; Stein, C.M.; Hartmann, K.; Schlesinger, L.S.; Sadee, W. AmpliSeq transcriptome analysis of human alveolar and monocyte-derived macrophages over time in response to Mycobacterium tuberculosis infection. PLoS ONE 2018, 13, e0198221. [CrossRef] 21. Coya, J.M.; De Matteis, L.; Giraud-Gatineau, A.; Biton, A.; Serrano-Sevilla, I.; Danckaert, A.; Dillies, M.-A.; Gicquel, B.; De la Fuente, J.M.; Tailleux, L. Tri-mannose grafting of chitosan nanocarriers remodels the macrophage response to bacterial infection. J. Nanobiotechnology 2019, 17, 15. [CrossRef] 22. Zhang, X.; Feng, H.; Li, Z.; Li, D.; Liu, S.; Huang, H.; Li, M. Application of weighted gene co-expression network analysis to identify key modules and hub genes in oral squamous cell carcinoma tumorigenesis. Onco. Targets. Ther. 2018, 11, 6001–6021. [CrossRef] [PubMed] 23. Wang, F.; Huang, S.; Gao, H.; Zhou, Y.; Lai, C.; Li, Z.; Xian, W.; Qian, X.; Li, Z.; Huang, Y.; et al. Initial Whole Genome Sequencing and Analysis of the Host Genetic Contribution to COVID-19 Severity and Susceptibility. Cell Discov. 2020. [CrossRef] [PubMed] Antibiotics 2021, 10, 97 15 of 16 24. Stelzer, G.; Rosen, N.; Plaschkes, I.; Zimmerman, S.; Twik, M.; Fishilevich, S.; Stein, T.I.; Nudel, R.; Lieder, I.; Mazor, Y.; et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr. Protoc. Bioinforma. 2016, 54, 1.30.1–1.30.33. [CrossRef] 25. Leprae, M.; Phelan, J.; Maitra, A.; Mcnerney, R.; Nair, M.; Gupta, A.; Coll, F.; Pain, A.; Bhakta, S.; Clark, T.G. The draft genome of Mycobacterium aurum, a potential model organism for investigating drugs against Mycobacterium tuberculosis and Mycobacterium leprae. Int. J. Mycobacteriology 2015, 4, 207–216. [CrossRef] 26. Namouchi, A.; Cimino, M.; Favre-Rochex, S.; Charles, P.; Gicquel, B. Phenotypic and genomic comparison of Mycobacterium aurum and surrogate model species to Mycobacterium tuberculosis: Implications for drug discovery. BMC Genom. 2017, 18, 530. [CrossRef] [PubMed] 27. Hickman, S.P.; Chan, J.; Salgame, P. Mycobacterium tuberculosis Induces Differential Cytokine Production from Dendritic Cells and Macrophages with Divergent Effects on Naive T Cell Polarization. J. Immunol. 2002, 168, 4636–4642. [CrossRef] 28. Cooper, A.M.; Mayer-Barber, K.D.; Sher, A. Role of innate cytokines in mycobacterial infection. Mucosal Immunol. 2011, 4, 252–260. [CrossRef] 29. Van Crevel, R.; Ottenhoff, T.H.M.; van der Meer, J.W.M. Innate Immunity to Mycobacterium tuberculosis. Clin. Microbiol. Rev. 2002, 15, 294–309. [CrossRef] 30. Boro, M.; Singh, V.; Balaji, K.N. Mycobacterium tuberculosis-triggered Hippo pathway orchestrates CXCL1/2 expression to modulate host immune responses. Sci. Rep. 2016, 6, 1–14. [CrossRef] 31. Dunlap, M.D.; Howard, N.; Das, S.; Scott, N.; Ahmed, M.; Prince, O.; Rangel-Moreno, J.; Rosa, B.A.; Martin, J.; Kaushal, D.; et al. A novel role for C–C motif chemokine receptor 2 during infection with hypervirulent Mycobacterium tuberculosis. Mucosal Immunol. 2018, 11, 1727–1742. [CrossRef] 32. Volpe, E.; Cappelli, G.; Grassi, M.; Martino, A.; Serafino, A.; Colizzi, V.; Sanarico, N.; Mariani, F. Gene expression profiling of human macrophages at late time of infection with Mycobacterium tuberculosis. Immunology 2006, 118, 449–460. [CrossRef] [PubMed] 33. Cohen, P. The TLR and IL-1 signalling network at a glance. J. Cell Sci. 2014, 127, 2383–2390. [CrossRef] 34. Xu, G.; Wang, J.; Gao, G.F.; Liu, C.H. Insights into battles between Mycobacterium tuberculosis and macrophages. Protein Cell 2014, 5, 728–736. [CrossRef] [PubMed] 35. Ishrat, R. Stage specific classification of DEGs via statistical profiling and network analysis reveals potential biomarker associated with various stages of TB. bioRxiv 2018, 414110. [CrossRef] 36. Ariel, O.; Gendron, D.; Dudemaine, P.L.; Gévry, N.; Ibeagha-Awemu, E.M.; Bissonnette, N. Transcriptome Profiling of Bovine Macrophages Infected by Mycobacterium avium spp. paratuberculosis Depicts Foam Cell and Innate Immune Tolerance Phenotypes. Front. Immunol. 2020, 10, 1–27. [CrossRef] 37. Goletti, D.; Lindestam Arlehamn, C.S.; Scriba, T.J.; Anthony, R.; Cirillo, D.M.; Alonzi, T.; Denkinger, C.M.; Cobelens, F. Can we predict tuberculosis cure? What tools are available? Eur. Respir. J. 2018, 52, 1801089. [CrossRef] 38. Walzl, G.; Ronacher, K.; Hanekom, W.; Scriba, T.J.; Zumla, A. Immunological biomarkers of tuberculosis. Nat. Rev. Immunol. 2011, 11, 343–354. [CrossRef] 39. Vrieling, F.; Kostidis, S.; Spaink, H.P.; Haks, M.C.; Mayboroda, O.A.; Ottenhoff, T.H.M.M.; Joosten, S.A. Analyzing the impact of Mycobacterium tuberculosis infection on primary human macrophages by combined exploratory and targeted metabolomics. Sci. Rep. 2020, 10, 7085. [CrossRef] 40. Wu, K.; Dong, D.; Fang, H.; Levillain, F.; Jin, W.; Mei, J.; Gicquel, B.; Du, Y.; Wang, K.; Gao, Q.; et al. An Interferon-Related Signature in the Transcriptional Core Response of Human Macrophages to Mycobacterium tuberculosis Infection. PLoS ONE 2012, 7, e38367. [CrossRef] 41. Zhang, W.; Bouchard, G.; Yu, A.; Shafiq, M.; Jamali, M.; Shrager, J.B.; Ayers, K.; Bakr, S.; Gentles, A.J.; Diehn, M.; et al. GFPT2- expressing cancer-associated fibroblasts mediate metabolic reprogramming in human lung adenocarcinoma. Cancer Res. 2018, 78, 3445–3457. [CrossRef] 42. Zhang, H.; Zhang, N.; Liu, Y.; Su, P.; Liang, Y.; Li, Y.; Wang, X.; Chen, T.; Song, X.; Sang, Y.; et al. Epigenetic regulation of NAMPT by NAMPT-AS drives metastatic progression in triple-negative breast cancer. Cancer Res. 2019, 79, 3347–3359. [CrossRef] [PubMed] 43. Rienksma, R.A.; Suarez-Diez, M.; Mollenkopf, H.-J.J.; Dolganov, G.M.; Dorhoi, A.; Schoolnik, G.K.; Martins dos Santos, V.; Kaufmann, S.H.E.; Schaap, P.J.; Gengenbacher, M.; et al. Comprehensive insights into transcriptional adaptation of intracellular mycobacteria by microbe-enriched dual RNA sequencing. BMC Genom. 2015, 16, 1–15. [CrossRef] [PubMed] 44. Mao, R.; Yang, R.; Chen, X.; Harhaj, E.W.; Wang, X.; Fan, Y. Regnase-1, a rapid response ribonuclease regulating inflammation and stress responses. Cell. Mol. Immunol. 2017, 14, 412–422. [CrossRef] [PubMed] 45. Gleeson, L.E.; Sheedy, F.J.; Palsson-McDermott, E.M.; Triglia, D.; O’Leary, S.M.; O’Sullivan, M.P.; O’Neill, L.A.J.; Keane, J. Cutting Edge: Mycobacterium tuberculosis Induces Aerobic Glycolysis in Human Alveolar Macrophages That Is Required for Control of Intracellular Bacillary Replication. J. Immunol. 2016, 196, 2444–2449. [CrossRef] 46. Petruccioli, E.; Scriba, T.J.; Petrone, L.; Hatherill, M.; Cirillo, D.M.; Joosten, S.A.; Ottenhoff, T.H.; Denkinger, C.M.; Goletti, D. Correlates of tuberculosis risk: Predictive biomarkers for progression to active tuberculosis. Eur. Respir. J. 2016, 48, 1751–1763. [CrossRef] Antibiotics 2021, 10, 97 16 of 16 47. Waddell, S.J.; Stabler, R.A.; Laing, K.; Kremer, L.; Reynolds, R.C.; Besra, G.S. The use of microarray analysis to determine the gene expression profiles of Mycobacterium tuberculosis in response to anti-bacterial compounds. Tuberculosis 2004, 84, 263–274. [CrossRef] 48. Boshoff, H.I.M.; Myers, T.G.; Copp, B.R.; McNeil, M.R.; Wilson, M.A.; Barry, C.E. The Transcriptional Responses of Mycobacterium tuberculosis to Inhibitors of Metabolism. J. Biol. Chem. 2004, 279, 40174–40184. [CrossRef] 49. Altaf, M.; Miller, C.H.; Bellows, D.S.; O’Toole, R. Evaluation of the Mycobacterium smegmatis and BCG models for the discovery of Mycobacterium tuberculosis inhibitors. Tuberculosis 2010, 90, 333–337. [CrossRef] 50. Rampacci, E.; Stefanetti, V.; Passamonti, F.; Henao-Tamayo, M. Preclinical Models of Nontuberculous Mycobacteria Infection for Early Drug Discovery and Vaccine Research. Pathogens 2020, 9, 641. [CrossRef] 51. Obregón-Henao, A.; Arnett, K.A.; Henao-Tamayo, M.; Massoudi, L.; Creissen, E.; Andries, K.; Lenaerts, A.J.; Ordway, D.J. Susceptibility of mycobacterium abscessus to antimycobacterial drugs in preclinical models. Antimicrob. Agents Chemother. 2015, 59, 6904–6912. [CrossRef] 52. Macho Rendón, J.; Lang, B.; Ramos Llorens, M.; Gaetano Tartaglia, G.; Torrent Burgas, M. DualSeqDB: The host–pathogen dual RNA sequencing database for infection processes. Nucleic Acids Res. 2021, 49, D687–D693. [CrossRef] [PubMed] 53. Bosshart, H.; Heinzelmann, M. THP-1 cells as a model for human monocytes. Ann. Transl. Med. 2016, 4, 438. [CrossRef] [PubMed] 54. Gupta, A.; Bhakta, S.; Kundu, S.; Gupta, M.; Srivastava, B.S.; Srivastava, R. Fast-growing, non-infectious and intracellularly surviv- ing drug-resistant Mycobacterium aurum: A model for high-throughput antituberculosis drug screening. J. Antimicrob. Chemother. 2009, 64, 774–781. [CrossRef] [PubMed] 55. Lu, L.; Arranz-Trullén, J.; Prats-Ejarque, G.; Pulido, D.; Bhakta, S.; Boix, E. Human Antimicrobial RNases Inhibit Intracellular Bacterial Growth and Induce Autophagy in Mycobacteria-Infected Macrophages. Front. Immunol. 2019, 10, 1500. [CrossRef] 56. Babraham Bioinformatics-FastQC A Quality Control tool for High Throughput Sequence Data. Available online: http://www. bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed on 21 August 2018). 57. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120. [CrossRef] [PubMed] 58. Pertea, M.; Kim, D.; Pertea, G.M.; Leek, J.T.; Salzberg, S.L. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 2016, 11, 1650–1667. [CrossRef] 59. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.-C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 2015, 33, 290–295. [CrossRef] 60. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [CrossRef] 61. Zhai, X.; Xue, Q.; Liu, Q.; Guo, Y.; Chen, Z. Colon cancer recurrence-associated genes revealed by WGCNA co-expression network analysis. Mol. Med. Rep. 2017, 16, 6499–6505. [CrossRef] 62. Zhou, G.; Stevenson, M.M.; Geary, T.G.; Xia, J. Comprehensive Transcriptome Meta-analysis to Characterize Host Immune Responses in Helminth Infections. PLoS Negl. Trop. Dis. 2016, 10, e0004624. [CrossRef] 63. Horvath, S. Integrated Weighted Correlation Network Analysis of Mouse Liver Gene Expression Data. In Weighted Network Analysis; Springer: New York, NY, USA, 2011; pp. 321–351. 64. Botía, J.A.; Vandrovcova, J.; Forabosco, P.; Guelfi, S.; D’Sa, K.; Hardy, J.; Lewis, C.M.; Ryten, M.; Weale, M.E. An additional k-means clustering step improves the biological features of WGCNA gene co-expression networks. BMC Syst. Biol. 2017, 11, 47. [CrossRef] [PubMed] 65. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. clusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. Omi. A J. Integr. Biol. 2012, 16, 284–287. [CrossRef] [PubMed] 66. Su, G.; Morris, J.H.; Demchak, B.; Bader, G.D. Biological Network Exploration with Cytoscape 3. Curr. Protoc. Bioinforma. 2014. [CrossRef] [PubMed]