Blood and bronchoalveolar lavage fluid (BALF) was collected from healthy adult volunteers, in strict accordance with the US Code of Federal Regulations and approved Local Regulations (The Ohio State University Human Subjects IRB numbers 2008H0135 & 2008H0119 and Texas Biomedical Research Institute/UT-Health San Antonio Human Subjects IRB numbers HSC20170667H & HSC20170315H), and Good Clinical Practice as approved by the National Institutes of Health (NIAID/DMID branch), with written informed consent from all human subjects. Healthy adult volunteers (18–45 years old) were TST- or IGRA-negative and were recruited from both sexes (50:50 male:female ratio) with no discrimination based on race or ethnicity. The following comorbidities were excluded: drug and alcohol users, smokers, asthma, acute pneumonia, upper/lower respiratory tract infections, acute illness/chronic condition, heart disease, diabetes, obesity, obstructive pulmonary disease (COPD), renal failure, liver failure, hepatitis, thyroid disease, rheumatoid arthritis, immunosuppression or taking nonsteroidal anti-inflammatory agents, human immunodeficiency virus (HIV)/AIDS, cancer requiring chemotherapy, leukemia/lymphoma, seizure history, blood disorders, depression, lidocaine allergies (used during the BAL procedure), pregnancy, nontuberculous mycobacterial infection, and TB.

BALF was collected from healthy adult donors in 80 mL of sterile endotoxin-free human isotonic saline (0.9% NaCl), filtered through 0.2 μm filters, and concentrated 20-fold using Amicon Ultra Centrifugal Filter Units with a 10-kDa molecular mass cutoff membrane (Millipore Sigma, Burlington, MA, USA) to obtain the ALF physiological concentration within the lung (1 mg/mL of phospholipid), as we described (6–13, 25). We use the in vivo relevant levels of 0.5 to 1.5 mg/ml of surfactant phospholipids in the lung mucosa (26), measured by phospholipid assay (27). The 20-fold concentration process removes free surfactant lipids and free fatty acids of <10-kDa, leaving functional hydrolases and other proteins such as surfactant proteins A and D (SP-A/SP-D), complement, and antimicrobial peptides in the ALF fraction, as well as other molecules of >10-kDa such as cholesterol. Effective lipid removal (<10-kDa) from ALF was confirmed in previous studies by mass spectrometry (MS) of fatty acid methyl ester (FAMES) derivatives (6). ALF was then aliquoted and stored at −80°C until further use. Samples were maintained at 4°C on ice and in endotoxin-free sterile conditions during the whole procedure.

M.tb strains CDC1551, H₃₇R_v (ATCC# 25618), HN878 (BEI Resources, NR-13647), and W-7642 were cultured from frozen stocks in 7H11 agar plates (BD BBL), supplemented with oleic acid, albumin, dextrose and catalase (OADC) for 14 days at 37°C. Single bacterial suspensions in 0.9% NaCl were obtained as described (6–12, 25), and bacterial pellets (~ 1 x 10⁹ bacteria/ml) obtained after centrifugation (13,000 x g, 10 min) were further incubated in a 50 μl pool of ALFs (mixture of ALFs from n = 17 different donors in equal proportions), 50:50 male:female ratio, no discrimination in ethnicity and race) for 15 min or ~12 h at 37°C (no agitation). After exposure, ALF was removed and bacterial pellets washed with 0.9% NaCl. Pellets were further incubated in RNAProtect (Qiagen, Hilden, Germany) for 10 min at room temperature (RT), centrifuged at 13,000 × g and stored at −80°C until further use. For each of the strains, unexposed bacteria were processed in parallel and used as controls (bacterial pellets incubated for 15 min or ~12 h at 37°C without ALF). Two independent experiments were performed: two different exposures to the same pool of ALFs on different days, using single bacterial suspensions of newly fresh harvested bacteria grown in 7H11 agar for each strain studied.

Stored bacterial pellets were used for RNA extraction using the Quick-RNA Fungal/Bacterial Miniprep kit (Zymo Research, Irvine, CA, USA), following the manufacturer's protocol. Briefly, bacterial pellet was resuspended in lysis buffer and transferred to a ZR BashingBead Lysis tube containing 0.1 mm and 0.5 mm ceramic beads. A bead beating procedure to break the mycobacterial cell envelope was performed in a BioSpec Mini-Beadbeater-24 (six cycles of 45 s at 3,800 rpm, with 3 min intervals on a cooler rack). RNA in the supernatant was isolated using Zymo-spin columns, with an in-column DNAse I treatment, and eluted in nuclease-free water. A second DNAse treatment was performed on the isolated RNA using the TURBO DNAse reagent (Thermo Fisher Scientific, Waltham, MA, USA), and RNA was further purified using the RNA Clean & Concentrator kit (Zymo Research, Irvine, CA, USA). RNA final concentration was measured with a Qubit 4 Fluorometer using the HS RNA kit (Thermo Fisher Scientific, Waltham, MA, USA), and RIN numbers were obtained using the 4200 TapeStation System (Agilent, Santa Clara, CA, USA) to determine RNA quality before sequencing. A stranded whole transcriptome RNA library was constructed using the Zymo-Seq RiboFree Total RNA kit (Zymo Research) following the manufacturer's guidelines. RNA (500 ng) was used as input to construct the library, which included a rRNA depletion step of 30 min followed by 14 cycles of PCR. A 75-bp paired-end (PE) sequencing was conducted using a NextSeq 500 Mid Output Illumina platform, obtaining 5–8 million PE reads per sample. All samples (from n = 2 independent experiments) were sequenced in the same flow cell to avoid sequencing batch effects.

Raw reads were trimmed (phred score of 30) and aligned to the reference genome H₃₇R_v (Genbank NC_000962.3) (28) using STAR v2.5.3a with default values in the Partek Flow Genomic Analysis software (v9.0.20.0202, Partek Inc., Chesterfield, MO, USA). Further analyses were performed in the iDEP.94 software (http://bioinformatics.sdstate.edu/idep94/) (29): first, uploaded read count raw data were filtered (0.5 CPM, n = 1 library) to remove lowly expressed genes and transformed with edgeR as log₂ (CPM + c) using default settings (30), and used for clustering and principal component Analysis (PCA). Hierarchical clustering of normalized gene expression levels was calculated using default values (correlation distance, average linkage, cutoff Z score of 4), and shown as a heatmap. Then, significant differentially expressed genes (DEGs) between ALF-exposed M.tb and unexposed M.tb (control) for each of the strains and time points were identified using the DESeq2 method (31) in iDEP.94, with an FDR cutoff of 0.1 and a minimum fold change (FC) of 2 (corresponding to a log₂ FC equal or greater than an absolute value of 1). The Benjamini-Hochberg procedure was used to adjust the false discovery rate (FDR), and fold changes were calculated using an Empirical Bayes shrinkage approach (31) and the Wald test implemented in DeSeq2. Functional categories for DEGs were extracted from Mycobrowser (32). Gene enrichment analysis was performed in ShinyGO 0.77 using default settings (FDR cutoff of 0.05 and minimum pathway size of 2, Pathway database: all available gene sets) (33). A summary of the experimental conditions and workflow used for RNA sequencing and data analysis is shown in Supplementary Figure S1.

Some of the genes that were differentially expressed were selected for validation using RT-qPCR (Supplementary Table S1). Selection criteria was based on genes discussed more in detail in the manuscript, belonging to different pathways and/or categories, and that showed significant differences between strains and/or time points (echA20, KstR2, mymA, papA1, papA3, mmpL10, mmpL8, mprA, sigF, espC and esxH). Strains selected for RT-qPCR where the ones that showed the most differential expression for each of the genes. Briefly, 500 ng of RNA (same RNA that was used for RNA-seq) were used for the synthesis of cDNA using the RevertAid H Minus First Strand cDNA Synthesis kit (Thermo Fisher Scientific, Waltham, MA, USA) with random hexamer primers, following the manufacturer's instructions. Then, cDNA and specific primers were used in a 20 μl qPCR reaction with PowerUp SYBR Green Master Mix (Applied Biosystems, Waltham, MA, USA) run in an Applied Biosystems 7500 Real-Time PCR instrument with the following settings: reporter SYBR Green, no quencher, passive reference dye ROX, standard ramp speed, and continuous melt curve ramp increment, as described (25). Relative expression was calculated using the 2^−ΔΔCT method (34) with sigA as the housekeeping gene.

Statistical significance for bacterial uptake and intracellular growth rate was determined by multiple paired Student's t tests (ALF-exposed vs. unexposed M.tb). A 2-way ANOVA for multiple comparisons with an uncorrected Fisher's LSD test and main effects model was used for qPCR data (ΔCt ALF-exposed vs. ΔCt unexposed M.tb for each of the strains and time points).

RNA from four M.tb strains (H₃₇R_v, CDC1551, HN878, and W-7642) exposed to human ALF for 15-min and 12-h was extracted and sequenced. The total average number of reads/sample after removing low quality reads during the initial quality control (QC) was 6.15 million, ranging from ~ 4.1 to 11.3 million. The number of aligned reads to the M.tb reference genome was ~ 90% on average. Samples with <70% of aligned reads were analyzed using Kraken (35) to discard potential sample contamination (data not shown). All samples had a coverage of more than 97% of their genome after alignment, with GC content of 62%. A Pearson correlation matrix was computed, showing a correlation of more than r = 0.97 between sample replicates (Supplementary Figure S2A). Our quality post-alignment statistics and quality control metrics are found in Supplementary Table S2.

Raw aligned data was then transformed and used for sample clustering and Principal Component Analysis (PCA). PCA grouped the samples based on M.tb strain (CDC1551, H₃₇R_v, HN878, W-7642) and incubation time (15-min and 12-h), the main sources of variance in this dataset, indicating temporal and strain-specific changes in transcriptional profiles of M.tb both before and after ALF exposure (Figure 1). Strain CDC1551 was furthest from the others at both time points tested, suggesting a more distinct transcriptional profile regardless of ALF exposure. We note here that the exposure time (15-min vs. 12-h) seems to have a bigger effect on the strains' transcriptome than the ALF exposure itself, since ALF-exposed and unexposed conditions are closely grouped in all the strains and time points. This indicates that the act of going from agar plates to being in suspension is causing a great number of transcriptional changes in M.tb, as one would expect. For this reason, unexposed bacteria was processed in parallel for each of the time points and used as control/baseline, and differential expression analyses were focused on the effects of ALF exposure (ALF-exposed vs. unexposed M.tb) at each time point rather than comparing exposure times.

Hierarchical clustering of normalized gene expression classified the samples in four main clusters (Supplementary Figure S2B). Cluster 1 (HN878 and W-7642 at 15-min) was closest to cluster 2 (HN878 and W-7642 at 12-h), indicating that those two strains presented the most similar transcriptional profile at each of the time points tested. They were followed by H₃₇R_v at both time points (cluster 3) and CDC1551 (cluster 4), furthest from the others.

We performed DE analysis of ALF-exposed vs. unexposed M.tb for each of the strains at both time points. A total of 397 DEGs were identified, with DR-M.tb strain W-7642 being the one with the most transcriptional changes after exposure to ALF for 12-h (152 downregulated and 87 upregulated genes) compared to unexposed bacteria (Figure 2A). Exposure for 15-min to ALF resulted only in a few DEGs: 32 in CDC1551, 12 in H₃₇R_v, 9 in HN878, and 6 in W-7642. Hierarchical clustering using DEGs clearly classified our samples in two main clusters based on the ALF exposure time (15-min vs. 12-h) (Figure 2B), indicating a temporal adaptation to the lung mucosa, with a much greater number of DEGs at 12-h compared to 15-min.

Identified DEGs belonged to different functional categories, including but not restricted to: conserved hypotheticals, intermediary metabolism and respiration, cell wall and cell processes, virulence, detoxification and adaptation, lipid metabolism, PE/PPE, information pathways, and regulatory proteins (Figure 3, Supplementary Table S3). At 15-min post-ALF exposure, the few DEGs found were associated to either lipid or intermediate metabolism, respiration, regulatory proteins or cell wall processes, except for strain CDC1551 which had more DEGs compared to the other strains and included other categories such as information pathways and PE/PPE (Supplementary Figure S3). At 12-h after ALF exposure, DEGs belonged to additional categories like virulence, detoxification and host adaptation, as well as insertion sequences and mobile genetic elements (MGEs) (Supplementary Figure S3).

We next performed gene enrichment analysis for the 397 total identified DEGs, and found 17 significantly enriched GO terms, including response to environmental pH, fatty acid biosynthetic process and lipid transport and associated beta-ketoacyl synthase activity, acyl transferase activity, cholesterol metabolic process, and mycofactocin and heme activity (Figure 4). Other enriched GO terms that did not pass the FDR cutoff of 0.05 were sulfolipid metabolic process, protein secretion by type VII secretion system, PPE superfamily, lipid biosynthesis, response to temperature and hypoxia, sigma factor activity, and phenolic phthiocerol biosynthetic process, among others (Supplementary Table S4).

The cholesterol catabolic pathway is highly conserved in mycolic acid-containing actinomycetes, and is encoded by more than 100 genes in M.tb (37). There are three distinct sub-pathways or processes during cholesterol catabolism, according to the three parts of the molecule: side chain degradation, A & B ring degradation, and C & D ring degradation (38). Two major regulons, controlled by two TetR-type transcriptional regulators, have been associated with these processes: KstR1 regulates the transcription of side chain and A & B ring degradation genes, while KstR2 controls the transcription of C & D ring degradation genes (39). Other genes though to be regulated by cholesterol belong to the mammalian cell-entry Mce3R regulon, implicated in stress resistance (40), and the SigE regulon, essential for virulence (41).

Our results indicate that most KstR2 genes (12/15) were significantly upregulated (Log₂FC >1 and FDR <0.1) either at 15-min or 12-h post-ALF exposure or both time points in one or more of the M.tb strains tested, with the exception of strain W-7642, where none of the KstR2 genes reached significant upregulation at 15-min (Figure 5, Supplementary Table S3). This included the genes for C & D ring degradation (3 out of 5 genes were significantly upregulated in at least one of the strains/conditions tested). In contrast, genes associated with the KstR1 regulon, including ATP-dependent cholesterol import genes encoded by the Mce4 transport system (42) and side chain and A&B ring degradation genes, were not differentially expressed when compared to unexposed M.tb at any of the ALF exposure periods of time tested (Supplementary Table S3). Only ltp3 from the KstR1 regulon, encoding a keto acyl-CoA thiolase involved in the cholesterol side chain degradation (43), was significantly downregulated in strain CDC1551 at 12-h post-ALF exposure. Other significant DEGs associated with cholesterol side chain degradation were fadE30, fadE32, echA20, and fadA6 (all upregulated and controlled by KstR2), and fadB2 (downregulated in all four strains after 12-h of ALF exposure) (Supplementary Table S3). In addition, of the 23 genes associated with the Mce3R regulon, only two were significantly differentially expressed: yrbE3B, encoding a hypothetical membrane protein with phospholipid transport activity, was upregulated in strain HN878 at 12 h post-exposure, while mce3F was downregulated in H₃₇R_v and W-7642 at that same time point (Supplementary Table S3). None of the cholesterol-associated sigE regulon genes were differentially expressed in our dataset under the experimental conditions tested.

In addition to the cholesterol utilization pathways, most M.tb strains showed upregulation of genes required for maintaining an appropriate mycolic acid layer composition and permeability of the M.tb cell envelope at 15 min post-ALF exposure, achieving statistical significance in strains CDC1551 and W-7642 (Figure 6, Supplementary Table S3). This included shared genes Rv3083 (mymA), Rv3084 (lipR), and Rv3085 (sadH), as well as Rv3086 (adhD), Rv3087, and Rv3088 (tgs4), only significant in W-7642. In contrast, these genes were downregulated after 12 h of ALF exposure (Figure 7A). All these genes belong to the mymA operon, associated with maintaining the M.tb cell envelope under harsh conditions (44, 45).

At 12-h post-ALF exposure, we found several up and downregulated genes shared among two or more of the M.tb strains evaluated. Interestingly, both hypervirulent HN878 and drug resistant W-7642 showed significant upregulation of diacyltrehalose (DAT), penta-acyltrehalose (PAT) and sulfolipid-1 (SL-1) biosynthesis and transport genes pks4, papA3, mmpL10, papA1, and mmpL8 (Figure 7B). These glycolipids are major M.tb cell envelope components (46), and some of them are significantly reduced on the M.tb cell envelope surface after exposure to human ALF (6), indicating a potential compensatory upregulation of these genes to restore the M.tb cell surface. Indeed, pks1, involved in the synthesis of immunomodulatory phenolic glycolipids (PGLs) produced by some M.tb strains and associated with cell envelope impermeability, phagocytosis, defense against oxidative stress and biofilm formation (46–48), was upregulated in all M.tb strains in this study, although was only significant in CDC1551 and W-7642 (Figure 7B, Supplemental Table S3). We note here, though, that many clinical isolates (including CDC1551), as well as laboratory strain H₃₇R_v, contain a pks1-15 polymorphism associated with the inability of these strains to synthesize PGL (49).

Rv1685c, a TetR-like transcriptional regulator that controls an efflux pump with a potential role in the response to acid nitrosative stress (50), was also significantly expressed in HN878 and W-7642 (Figure 7B). Other upregulated genes at 12 h post-ALF exposure were associated to mycobactin assembly (mtbB, mtbC), involved in iron acquisition, growth in macrophages and virulence (51, 52), and stress response gene hspX in H₃₇R_v and HN878, as well as potential immunogenic and virulence-associated genes (Rv2190c, Rv2624c) in CDC1551 and H₃₇R_v.

We also found several shared downregulated genes after 12-h of ALF exposure (Figure 7A). Genes associated with resistance to aminoglycosides (acc, Rv0263c, Rv0264c) (54), lipid metabolism (icl1, fadB2) (55), and osmotic stress response and persistence (Rv0516c, mprA) (56, 57) were consistently downregulated in all four M.tb strains. Further, fatty acid (FA) transporter mmpL12 was down in all M.tb strains but CDC1551. We also observed downregulation of mycofactocin synthesis genes, which have a potential multifaceted role in alcohol metabolism of M.tb, hypoxia adaptation, and redox homeostasis (58), as well as several lipoproteins and toxin-antitoxin (TA) systems in most of the M.tb strains, among others (Figure 7A).

In addition to shared transcriptional responses among M.tb strains due to ALF exposure, we also observed unique expression signatures attributed to the M.tb genomic background at both time points tested, shown in Supplementary Table S5. At 15-min post-ALF exposure, most of the strain-specific DEGs were still related to the KstR2 regulon and maintaining the mycolic acid layer composition and permeability of the M.tb cell envelope. We found a gene encoding a probable transmembrane protein responsible for the transport of glutamine across the membrane (Rv0072) significantly increased in the hypervirulent HN878 (log₂FC >5). The highly transmissible CDC1551 showed the most changes at this time point, with 25 unique DEGs (Figure 6). It included downregulation of essential mycolic acid biosynthesis genes acpM, kasA, and kasB (59, 60), pks16 involved in biofilm formation (61), and cell envelope-associated genes such as rpfA (stimulates resuscitation of dormant cells), Rv0888 (induces formation of neutrophil extracellular traps or NETs) (62), and glnA1 (regulates nitrogen levels during poly-L-glutamine or PLG synthesis in the cell envelope) (63, 64). Interestingly, a global transcriptional regulator (Rv2359 or zur) that represses transcription of genes involved in zinc homeostasis was upregulated in CDC1551 (Figure 6).

At 12-h post-ALF exposure, CDC1551 downregulated several cell envelope genes, including Rv0559c and phoY, with suggested roles in drug resistance (65, 66). Sigma factor E (sigE), involved in persistence to anti-TB drugs and essential for virulence and phagosome maturation in macrophages, was downregulated (67, 68), as well as msrB, which encodes a repair enzyme for proteins that are inactivated by oxidation and reported as upregulated in the resuscitation of dormant M.tb (69). In addition, several PE/PPE genes, involved in virulence and the activation of proinflammatory responses, and stress-induced genes (e.g., hspR, clpB) were also downregulated. In contrast, only a few unique genes were upregulated in CDC1551, including Rv2958c and Rv2959c, involved in the glycosylation of cell envelope PGLs, and some fatty acid biosynthesis genes, which were initially downregulated at 15-min ALF post-exposure.

Interestingly, H₃₇R_v showed upregulation of several latency-associated dosR regulon genes after 12-h of ALF exposure (Rv2623 or TB31.7, Rv2625c, hrp1) (70), as well as pncB2, involved in NAD biosynthesis during host invasion that was found increased during hypoxia (71), indicating a potential role during persistence (Supplementary Table S5). Gene irtB, involved in iron homeostasis (72, 73), and a couple of genes associated with resistance to oxidative stress (lpdA, Rv1050) were also upregulated. Similar to CDC1551, we found two sigma factors with decreased expression (sigB and sigD), potentially involved in the regulation of the stationary phase of M.tb (74). Other unique downregulated genes included pup, with a role in proteasomal degradation, and TA-associated genes (Supplementary Table S5).

Of the 28 unique DEGs found in hypervirulent HN878 at 12-h post-ALF exposure, 24 showed significant increased expression (log₂FC >1 and FDR <0.1) (Supplementary Table S5). Within the cell envelope functional category, we found esxI upregulated, which encodes an ESAT-6 protein secreted by the ESX-5 system that is involved in the virulence of M.tb (75). We also found upregulated the lipoprotein lppJ of unknown function, and the cluster Rv1686c/Rv1687c, a putative efflux pump responsible for the transport of undetermined substrates (potentially drugs) across the membrane that is activated under several environmental stressors (76, 77). Further, unique DEGs from the intermediate metabolism and respiration category were all upregulated in HN878, including glbN that encodes an oxygen transport protein upregulated during macrophage infection (78), molybdenum cofactor biosynthesis protein MoaA1, strongly induced under hypoxia (79), and dxs2, involved in the mycobacterial ability to prevent acidification of the phagosome (80) (Supplementary Table S5). Contrary to CDC1551 and H₃₇R_v, in HN878 two PE/PPE genes showed increased expression. The only unique downregulated genes in HN878 were two transmembrane proteins of unknown function, an uncharacterized hypothetical protein, and fadE8, involved in lipid degradation (Supplementary Table S5).

Finally, MDR-M.tb W-7642 showed the most strain-specific DEGs (#155) after 12-h post-ALF exposure compared to unexposed bacteria, with 66 upregulated and 89 downregulated genes (Figure 7). We found upregulation of several genes associated with the ESX-3 system and other genes related to iron acquisition and utilization (esxH, eccD3, mycP3, Rv0285 or PE5, hemE, hemL, Rv2047c), essential for mycobacterial growth in vitro and virulence (81) (Supplementary Table S5). In contrast, ESX-1 genes espD, espC, and espA were downregulated. The ESX-1 system secretes host membrane-targeting proteins such as ESAT-6 that interrupt the innate immune response against M.tb (82, 83). Other ESAT-6-like proteins were found downregulated (EsxJ, EsxP, EsxD) in W-7642, with the exception of EsxL that was upregulated. Overall, genes generally induced during starvation, persistence, and/or during the stationary phase of M.tb showed decreased expression (e.g., rsbW, vapC44, mce1R, lipU, Rv3241c) (84), as well as sigma factors sigF and sigI (Supplementary Table S5). However, we observed increased expression of M.tb cell envelope biosynthesis genes and lipid metabolism (Rv2974c, aftD, murE, murF, Rv3807c, accD6, Rv3031, kasB, acpS, fas, ppsC, ppsD, pks15), amino acid metabolism (leuC, leuD, argD, trpA), and DNA repair genes (ruvB, Rv3201c, Rv3433c, radA). Interestingly, W-7642 was the only strain with decreased expression of ethA (85), but increased expression of Rv2989 (86), both associated with drug resistance. Other unique DEGs in strain W-7642 are shown in Supplementary Table S5.

RNA-seq results for selected genes of interest were validated by qPCR, showing similar fold change values (Supplementary Figure S4).