Comparative transcriptomic analysis of THP‐1‐derived macrophages infected with Mycobacterium tuberculosis H37Rv, H37Ra and BCG

Tuberculosis (TB) remains a worldwide healthcare concern, and the exploration of the host‐pathogen interaction is essential to develop therapeutic modalities and strategies to control Mycobacterium tuberculosis (M.tb). In this study, RNA sequencing (transcriptome sequencing) was employed to investigate the global transcriptome changes in the macrophages during the different strains of M.tb infection. THP‐1 cells derived from macrophages were exposed to the virulent M.tb strain H37Rv (Rv) or the avirulent M.tb strain H37Ra (Ra), and the M.tb BCG vaccine strain was used as a control. The cDNA libraries were prepared from M.tb‐infected macrophages and then sequenced. To assess the transcriptional differences between the expressed genes, the bioinformatics analysis was performed using a standard pipeline of quality control, reference mapping, differential expression analysis, protein‐protein interaction (PPI) networks, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Q‐PCR and Western blot assays were also performed to validate the data. Our findings indicated that, when compared to BCG or M.tb H37Ra infection, the transcriptome analysis identified 66 differentially expressed genes in the M.tb H37Rv‐infected macrophages, out of which 36 genes were up‐regulated, and 30 genes were down‐regulated. The up‐regulated genes were associated with immune response regulation, chemokine secretion, and leucocyte chemotaxis. In contrast, the down‐regulated genes were associated with amino acid biosynthetic and energy metabolism, connective tissue development and extracellular matrix organization. The Q‐PCR and Western blot assays confirmed increased expression of pro‐inflammatory factors, altered energy metabolic processes, enhanced activation of pro‐inflammatory signalling pathways and increased pyroptosis in H37Rv‐infected macrophage. Overall, our RNA sequencing‐based transcriptome study successfully identified a comprehensive, in‐depth gene expression/regulation profile in M.tb‐infected macrophages. The results demonstrated that virulent M.tb strain H37Rv infection triggers a more severe inflammatory immune response associated with increased tissue damage, which helps in understanding the host‐pathogen interaction dynamics and pathogenesis features in different strains of M.tb infection.

In particular, during these processes, various bacterial metabolites and by-products can have some effects; these metabolites and byproducts thus serve as a crosstalk bridge between bacteria and host cells, 4  To address whether differential host responses or variances in the responses to the different strains of M.tb lead to differential disease severity, our current study aimed to investigate the immune responses induced by different genotypes strains of M.tb in macrophages by transcriptomic analysis. This will help to understand the underlying mechanisms of the virulence of this economically important pathogen.
Our findings indicated that the gene expression profiles of H37RV, H37Ra and BCG infections were significantly different. In comparison with H37Ra and BCG infections, 66 unique genes were expressed in response to H37RV M.tb infection, with 36 genes up-regulated and 30 genes down-regulated. Increased expression of pro-inflammatory cytokines and chemokines was observed in H37RV-infected macrophages; these diverse gene expression profiles were associated with altered activation of inflammation-related signalling pathways, including STAT1, NF-kB and MAPK; additionally, we observed that, rather than inducing cell apoptosis, H37RV infection in macrophages induces a significantly higher level of pyroptosis.

| Bacterial culture condition
The axenic cultures of Mycobacterium tuberculosis H37Rv, H37Ra and BCG strains were employed in all experiments. M.tb strains were grown routinely at 37ºC in 10 ml of 7H9 (BD Middlebrook, Difco, Sparks), supplemented with 10% of oleic acid albumin dextrose catalase (OADC) (Becton Dickinson Microbiology Systems), 0.5% Glycerol (Sigma) and 0.05% of Tyloxapol (Sigma). Growth was monitored daily by OD absorbance at 600 nm for 20 days.

| Infection of the THP-1 macrophage cell line
The THP-1 cell line was maintained in medium 1640 supplemented with 10% (v/v) foetal calf serum, 2 mM glutamine, 1 mM strain H37Rv infection triggers a more severe inflammatory immune response associated with increased tissue damage, which helps in understanding the host-pathogen interaction dynamics and pathogenesis features in different strains of M.tb infection.

K E Y W O R D S
infection, macrophage, Mycobacterium tuberculosis, regulation, transcriptome sodium pyruvate and 0.05 mM 2-mercaptoethanol in a humidified 5% carbon dioxide atmosphere at 37°C. The THP-1 was differentiated for 24 h using a culture medium containing 40 ng/ml phorbol 12-myristate 13-acetate. The cells were then washed with fresh culture medium and incubated for 48 h. Rapidly grown M. tuberculosis H37Rv, H37Ra and BCG were pelleted (3,200 rpm, 10 min), washed twice with PBS (phosphate-buffered saline) and re-suspended in a THP-1 culture medium. Differentiated THP-1 cells were optimal conditions infected with Mycobacterium tuberculosis at a multiplicity of infection of 5 (MOI) for 24h. 9 The experiments were divided into four groups, including: control, H37Rv, H37Ra, and BCG-infected groups, respectively.

| Macrophages RNA isolation and sequencing
The total RNA extraction was prepared from M.tb-infected THP-1 cells with RNA Miniprep kit reagent as described by the manufacturer Vazyme biotech. NGS analysis was performed on the BGISEQ-500 platform by BGI Genomic Services, generating per sample on average about 6.86G reads. After filtering the reads, using HISAT clean reads were mapped to the reference genome.
Averagely, 90.73% reads were mapped, and by using Bowtie2 (v2.2.5), clean reads were mapped to reference transcripts; then with RSEM (v1.2.12), the level of gene expression for each sample was measured. 10 We used Interferome 2.0 and EnrichR for annotation of transcripts; by CIMminer, clustered image maps (CIMs) were rendered. In the NCBI's Gene Expression Omnibus database, NGS data were deposited. The sequence outcomes for each transcript were obtained as the FPKM (fragment per kilobase of exons per million reads).

| RNA-seq data analysis
The raw data were defined as reads containing the low-quality reads, high content of unknown bases and the sequence of the adaptor.
To decrease data noise, these reads were removed before downstream analysis. The steps of filtering are as follow: The RNA-seq libraries' quality was first assessed using the Fast QC v0.11.5 software. Then, according to the following parameters, the reads were subjected to standard criteria of quality control (QC): (1) trimming and cleaning reads that aligned to primers and/or adaptors, (2) reads with over 50% of low-quality bases (quality value ≤5) in one read, and (3) reads with over 10% unknown bases (N bases). After filtering, the remaining reads are called 'clean reads' and stored as FASTQ format. The present study used bioinformatics methods to identify

| Real-Time PCR and Western blot
Total RNA extraction was prepared from M.tb-infected THP-1 cells. The housekeeping gene GAPDH was used as an internal control. THP-1 mRNA primer sequences for the target genes are listed in supplementary Table 6

| Statistical analysis
All data were expressed as the mean ±SD from at least three independent experiments. When the distribution was normal and there was homogeneous variance, we used the ANOVA to compare amongst three or more groups. P-values of 0.05 or less were considered statistically significant. All analyses were performed with Graph Pad Prism 7.0 software (version 6.07).

| Quality control of sequencing data
The quality distribution of bases reflects the accuracy of sequencing reads. We first assessed the RNA sequencing Quality Control parameters including, sequencing yield, Q20 and Q30 (% of bases greater than Q30 is a metric for the overall quality and accuracy of the base calling of the sequencing run), Total Raw reads, Total Clean reads and Total clean bases. Overall, we can see that the ratio of bases with low quality (quality <20) is low, indicating the good quality of sequencing (Supplementary Figure 1A). The base content of each position of reads is stable, without AT or GC separation. In RNA sequencing of the Illumina platform, the 6-bp random primers used in reverse transcription of cDNA will make the base composition of the first few positions has a preference, so the fluctuation of the first 6-bp base proportion in the figure is a normal phenomenon. The AT and GC content of all our samples obeyed this rule (Supplementary Figure 1B). The sequencing quality of all the samples is shown in Figure 1C. We can see that the Q30 of all bases was 86.5%-87.5%, and Q20 exceeded 95%, showing that the sequencing accuracy was very high (Supplementary Figure 1C). In all the reads of our samples, the proportion of clean reads exceeded 90%, showing a high quality of sequencing (Supplementary Figure 1D). We further performed Venn diagram analysis on 6 groups of differentially expressed genes (DEGs) (Supplementary Figure 2C). In this regard, induction of host innate immune response in macrophages can be confirmed based on the data that the numbers of DEGs were detected in each group. There were several specific or common DEGs that existed in each group, which is convenient for our next screening work.

| Transcriptomic profile of M.tb H37Rvinfected macrophages vs -uninfected Control
In total, 635 DEGs were found in the M.tb H37Rv-infected macrophages, of which 322 genes were up-regulated, and 313 genes were down-regulated in the heat map ( Figure 1A). Next, we per- H37Ra-infected macrophage, versus the mock infection control, to identify hub nodes. We analysed 515 up-and 526 down-regulated genes by PPI and then analysed them according to node and interaction intensity. As shown in Figure 4A, in M.tb H37Ra-infected macrophages, we identified HIST2H2BE and MT1 family proteins were up-regulated in PPI. In addition, the PPI analysis of down-regulated genes was also compared. According to Figure 4B, it was found that IL-10 (degree = 10), IFIT1 (degree = 26), NCAM1 (degree = 15), CXCL5 (degree = 12), CXCL13 (degree = 8), CTSG (degree = 18) and CCR1 (degree = 12) were the listed hub nodes within the downregulated DEGs.

| Transcriptomic profile of BCG-infected macrophages vs -uninfected Control
Our results revealed that compared with the control group, a total of  Figure 5D).  Figure 6, the selected network showed the signalling between cytokines and chemokines affecting T-cell profiling and apoptosis in immune cells. According to Figure 6A, it was found that the increased expression of ESR1 (degree = 9), NTN1 (degree = 11) and TGFA (degree = 5); the decreased expression of IL17R (degree = 4), PTAFR (degree = 15), MELK (degree = 12), GBP5 NLRP3 (degree = 13), PYGER3 (degree = 19) and IFIT family proteins were the main features of the PPI network ( Figure 6B).

| Validation of DEGs and signalling pathway activation in H37RV, H37Ra and BCGinfected macrophage
To verify the transcriptome analysis data, we further analysed the DEGs and signalling transduction pathway activation in H37Rv,

H37Ra and BCG-infected macrophages by Q-PCR and Western blot.
As shown in Figure 9, in Q-PCR analysis, there was a significant increase in several pro-inflammatory cytokines, including IL-1β and TNFα in H37Rv-infected macrophages, compared with H37Ra and BCG infection (Figure 9). Interestingly, we found that the expression of immunosuppressive factors, such as IL-10, was significantly downregulated upon H37Rv infection (Figure 9). In terms of chemokine expression, CCL3, CCL4 and CXCL8 levels were significantly upregulated in response to H37Rv infection, whereas CCL2 levels were significantly down-regulated ( Figure 9). Other immune regulatory molecules, such as CD14 and CSF2, were significantly up-regulated, while others, such as TLR4, IRF9 and CD36, were significantly downregulated after infection with H37Rv ( Figure 9).
In the transcriptome analysis, through GO and KEGG analysis of DEGs, we noticed that the signal pathway and cell function of cellular inflammatory response were significantly activated after infection with M.tb, especially in H37Rv. Therefore, we further validated the inflammatory-related signalling pathway activation by measuring the pathway-related protein expression by Western blot. As shown in Figure 10, the expression levels of p-P65, p-ERK and p-JNK, which are associated with the activation of NF-κB and MAPK related pathways, as well as the activation of STAT1, were significantly higher in H37TV infected cells compared with H37Ra or BCG infection, indicating that the virulent strain of MTB infection is likely to trigger a higher level of inflammation, which is consistent with our Q-PCR and the transcriptome analysis data ( Figure 10).
Meanwhile, we investigated the mode of cell death following infection by measuring the expression of caspase-1, cleaved-caspase-1, Caspase-3 and Cleaved-Caspase-3. We discovered a significant increase in cleared-caspase-1 expression in H37RV-infected macrophages when compared to H37Ra or BCG infection; however, the expression of Cleared-Caspase-3 was not detectable in all cases, indicating that pyroptosis, rather than apoptosis, occurred after virulent MTB infection ( Figure 10).

| DISCUSS ION
The Macrophage polarization is increasingly known as an influential pathogenesis factor in multiple types of infectious and inflammatory diseases. Classical M1 and alternative M2 polarization are the two extremes of macrophage phenotype transformation. 18 Pro-inflammatory M1 macrophages show anti-infection activity by promoting T helper (Th) 1 response. M2 macrophages promote Th2 responses and contribute to tissue repair, which are supposed to be anti-inflammatory. 19 In accordance with previous studies, 16 we confirmed that H37Rv infection in THP-1-derived macrophages induced increased TNF-a expression and indicated the M1 polarization potential; compared with BCG or H37Ra infection, H37Rv infection selectively induced IL-23 rather than IL-12 in THP-1-derived F I G U R E 7 The heat map, Venn diagram, GO and KEGG classification of the DEGs screening. A, Shows the number of DEGs in three different datasets and the crossing area indicates the cross-DEGs in different datasets. Including NC vs H37Rv (orange), NC vs H37Ra (green) and NC vs BCG (red). B, The heat map shows the relative transcript levels of the 133 differential genes in THP-1 cells infected with Mycobacterium tuberculosis H37Rv, HA36Ra and BCG strains. C, Common 133 differential genes GO biological process enrichment analysis of significant DEGs. D, Common 133 differential genes KEGG pathway enrichment analysis of significant DEGs macrophages. This result indicated that H37Rv is more likely to  In Leishmania-infected macrophages to establish an antiinflammatory response, the expression and phosphorylated status of Ago2 protein regulated miRNA-targeted pro-inflammatory cytokines. 21 As for M.tb infection, single nucleotide polymorphisms (SNPs) of AGO2 have also been found linked with tuberculosis susceptibility in the Chinese Uygur population. 28 Here, we detected increased expression of the AGO2 gene in H37Rv-infected macrophages, further supporting its anti-M.tb function. The ATP10A gene belongs to the family of P-type cation transport ATPases. Its functions include the translocation of phosphatidylcholine, and it is involved in plasma membrane dynamics. Its function in M.tb infection is not fully known. However, it has been identified as a novel biomarker for the discrimination of TB from latent TB infection individuals (LTBI). 23

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in GEO reference number GSE162729.