Transcriptomic analysis of longitudinal Burkholderia pseudomallei infecting the cystic fibrosis lung

The melioidosis bacterium, Burkholderia pseudomallei, is increasingly being recognised as a pathogen in patients with cystic fibrosis (CF). We have recently catalogued genome-wide variation of paired, isogenic B. pseudomallei isolates from seven Australasian CF cases, which were collected between 4 and 55 months apart. Here, we extend this investigation by documenting the transcriptomic changes in B. pseudomallei in five cases. Following growth in an artificial CF sputum medium, four of the five paired isolates exhibited significant differential gene expression (DE) that affected between 32 and 792 genes. The greatest number of DE events was observed between the strains from patient CF9, consistent with the hypermutator status of the latter strain, which is deficient in the DNA mismatch repair protein MutS. Two patient isolates harboured duplications that concomitantly increased expression of the β-lactamase-encoding gene penA, and a 35 kb deletion in another abolished expression of 29 genes. Convergent expression profiles in the chronically-adapted isolates identified two significantly downregulated and 17 significantly upregulated loci, including the resistance-nodulation-division (RND) efflux pump BpeEF–OprC, the quorum-sensing hhqABCDE operon, and a cyanide- and pyocyanin-insensitive cytochrome bd quinol oxidase. These convergent pathoadaptations lead to increased expression of pathways that may suppress competing bacterial and fungal pathogens, and that enhance survival in oxygen-restricted environments, the latter of which may render conventional antibiotics less effective in vivo. Treating chronically adapted B. pseudomallei infections with antibiotics designed to target anaerobic infections, such as the nitroimidazole class of antibiotics, may significantly improve pathogen eradication attempts by exploiting this Achilles heel.


INTRODUCTION 43
The Gram-negative soil-dwelling bacterium Burkholderia pseudomallei causes melioidosis, an 44 opportunistic tropical infectious disease of humans and animals that has a high fatality rate 45 (Wiersinga et al. 2012). B. pseudomallei is found in many tropical and subtropical regions 46 globally, and has been unmasked in temperate and even arid environments following unusually 47 wet weather events (Yip et al. 2015;Chapple et al. 2016;Sarovich et al. 2016). Infection occurs 48 following percutaneous inoculation from contaminated soil or water, inhalation, or ingestion. 49 Melioidosis symptoms vary widely due to the ability for B. pseudomallei to infect almost any 50 organ, with pneumonia being the most common presentation (Leelarasamee and Bovornkitti 51 1989; Currie et al. 2010). Individuals most at risk of contracting melioidosis include diabetics, 52 those with hazardous alcohol consumption, and the immunosuppressed. There has been 53 increasing recognition that people with chronic lung diseases such as cystic fibrosis (CF) are 54 also at a heightened risk (Holland et al. 2002;O'Carroll et al. 2003). 55 CF is a heritable disorder of the CFTR gene, and defects in CFTR lead to exaggerated and 56 ineffective airway inflammation, an imbalance in salt regulation in the lungs and pancreas, and a 57 chronic overproduction of thick and sticky mucus in the airways and digestive system (Amaral 58 al. 2008). The most common and best-studied CF pathogen is P. aeruginosa, which can adapt 67 to the CF lung environment via various mechanisms. Convergent pathoadaptations in P. presentation is chronic carriage (76%), which is associated with accelerated lung function 79 decline (Geake et al. 2015). This prevalence contrasts with melioidosis in non-CF patients, 80 where chronic carriage is exceedingly rare, occurring in only 0.2% of cases (Price et al. 2015). 81 To better understand B. pseudomallei pathoadaptation in the CF lung, we recently investigated 82 the genome-wide evolution of isogenic B. pseudomallei strains isolated from seven Australasian 83 CF patients, which were collected between 4 and 55 months apart (Viberg et al. 2017). 84 Hallmarks of these infections included B. pseudomallei persistence despite multiple eradication 85 attempts, multidrug resistance, mutations in virulence, metabolism and cell wall components, 86 and the first-documented case of hypermutation in B. pseudomallei. In all except one case, 87 multiple single-nucleotide polymorphism (SNP) and insertion-deletion (indel) mutations were 88 identified, with a high rate of nonsynonymous mutations, many of which were predicted to affect 89 protein function (Viberg et al. 2017). 90

RNA-seq provides a detailed view of the transcriptional landscape in bacterial isolates grown 91
under different conditions or niches (Sharma et al. 2010), and is now a well-established method 92 for examining differential gene expression (DE) in bacterial pathogens (Creecy and Conway 93 2015). Here, we performed bacterial RNA-seq on five of the CF cases that we have recently

RESULTS 100
Differential gene expression among CF isogenic pairs. DE was observed in four of the five 101 CF pairs, with only the CF10 pair failing to yield significant transcriptional differences (i.e. no 102 genes with a ≥1.5 log2-fold change and a false discovery rate (FDR) of ≤0.01; Figure 1). We 103 have previously shown that no genetic variants separate the CF10 strains, which had the 104 shortest time between collection of only 10 months (Viberg et al. 2017). This lack of significant 105 transcriptional differences rules out epigenetic effects on gene expression between this pair, at 106 least under the tested growth conditions, and illustrates that RNA-seq is a robust methodology 107 that is not readily prone to false-positive results.  SNPs, indels, deletions or gene duplications) identified in CF8, CF11, CF6 and CF9, 124 respectively (Table 1). The elevated number of mutations seen in CF9 is due to a mutS 125 mutation in the latter strain, which confers a hypermutator phenotype, the first time 126 hypermutation has been described in B. pseudomallei (Viberg et al. 2017); this in turn 127 contributes to a high number of DE genes. However, when comparing the ratio of DE genes to 128 mutational events, CF11 had the highest proportion of DE genes (11.3), followed by CF6 (9.5), 129 CF9 (7.1) and CF8 (2.7). There was a significant skew in DE towards genes located on 130 chromosome II, which contains a lower proportion of housekeeping genes than chromosome I 131 (Holden et al. 2004 with DE genes in the latter isolates from CF9, CF6, CF11 and CF8, respectively, and when 167 compared with genes lost in MSHR6686, there was a 9, 15, 20 and 97% overlap, respectively 168 (Table S1). The proportion of downregulated genes varied across this dataset, ranging from 169 13% for CF11 to 100% for CF8, demonstrating that the effect on gene expression at these loci 170 is not unidirectional, with certain overlap loci in fact being upregulated in CF6, CF9 and CF11. 171 Of note, all 29 genes (BPSS1131-BPSS1159) that were lost due to a large deletion in the latter 172 with forty-six downregulated surface polysaccharide genes (Table S1). 188 The LPS loci wbiE (BPSL2676) and wbiD (BPSL2677) were downregulated by ~1.8-fold (3x) in 189 the latter CF9 isolate. In addition, the poorly characterized LPS biosynthesis-related membrane 190 protein loci BPSS1683-BPSS1685 were downregulated by ~5.9-fold (60x). This isolate also 191 exhibited downregulation of all CPS-I loci (except wcbC) by between 1.7-and 3.1-fold (3 to 8x). 192 In contrast, the latter isolate from CF11 showed upregulation of the CPS-I loci BPSL2793-193 BPSL2797 (wcbN-wcbM-gmhA-wcbL-wcbK) when compared with its initial isolate, with 194 increases ranging from 2.7 to 6.1-fold (7 to 69x). However, when compared with initial isolates 195 from CF6, CF8, CF9 and CF10, expression of BPSL2793-BPSL2797 in the latter CF11 strain 196 was in fact downregulated (3.4 to 4.5-fold; 11 to 23x). This observation was confirmed as 197 significant downregulation of these loci in the initial CF11 strain (by between 6.1-and 10.1-fold; 198 69 to 1,098x) when compared with all other initial strains, rather than significant upregulation of 199 these CPS-I loci in the latter CF11 strain. 200 Unlike CPS-I, expression of the CPS-II cluster (BPSS0417-BPSS0429) is induced when grown 201 in water, suggesting that this polysaccharide plays a role in environmental survival  Zenteno et al. 2010). One locus involved in CPS-II biosynthesis, BPSS0425, was 203 downregulated (1.8-fold; 4x) in CF6, and the entire cluster was downregulated in CF9 (range: 204 2.5-to 5.2-fold; 6 to 37x). Conversely, BPSS0417 and BPSS0418 were upregulated in CF11 205 (1.6-and 1.7-fold; 3x), respectively. However, as with CPS-I, both CF11 strains exhibited 206 significant downregulation of BPSS0417 and BPSS0418 when compared with initial strains from 207 CF6, CF8, CF9 and CF10 (2.0-and 3.1-fold; 4 to 9x). The genes encoding CPS-III (BPSS1825-208 BPSS1835) were significantly downregulated in the latter isolates from CF6 (2.4-to 3.1-fold; 5 209 to 9x) and CF9 (5.6-to 7.5-fold; 50 to 178x). Finally, two genes within the CPS-IV cluster 210 (BPSL2769-BPSL2785) were downregulated in CF9 (BPSL2782 and BPSL2785; 1.8-to 2.8-211 fold; 3 to 7x) but 11/17 loci from this cluster were upregulated in CF11 by 1.7-to 5.9-fold 212 (BPSL2769, BPSL2775-BPSL2784). However, unlike CPS-I and CPS-II, this pattern of 213 upregulation in both CF11 isolates was maintained for CPS-IV loci BPSL2769 and BPSL2775-214 BPSL2781 even when compared with the initial CF6, CF8, CF9 and CF10 isolates (2.0-to 4.0-215 fold; 4 to 16x). virulence genes lacking genetic mutations were found to be significantly downregulated in the 228 latter CF isolates. These loci included three Type IV pilus 7 (TFP7) loci (pilR, pilG and pilN; 1.7-229 fold; 3x), the lysozyme inhibitor BPSL1057 (3.4-fold; 11x), Burkholderia lethal factor 1 (3.2-fold; 230 10x), four T3SS-3 loci (bsaS, bsaP, bsaO and bsaN), and 16 flagellum loci in CF9 (average 2.5-231 fold; 6x), and a trimeric autotransporter adhesin (bpaC; BPSL1631) in CF11 (1.6-fold; 3x). Of 232 these, the four T3SS-3 loci are also missing in chronic P314 isolates. 233

Several regulators with decreased transcription in the CF isolates are absent in B. mallei. 234
We hypothesized that downregulation of transcriptional regulators, particularly those absent in 235 B. mallei, would be identified in the latter CF isolates due to niche adaptation. As predicted, two 236 of the four pairs exhibited significant downregulation of transcriptional regulators. The first of 237 these, the Fis family regulatory protein YfhA (encoded by BPSL0350), was downregulated by 238 ~2.7-fold (~7x) in both CF6 and CF9. In E. coli, Fis is a global regulator that is induced under 239 nutrient-rich conditions and plays a role in the regulation of myriad processes including the 240 initiation of DNA replication, ribosomal RNA transcription activation and capsule expression 241 The latter isolate from CF6, which encodes a T314fs mutation in bpeT and is highly resistant 259 towards TMP/SMX (MIC ≥32 µg/mL), showed 5.2 to 6.8-fold upregulation of bpeEF-oprC (38 to 260 111x; Table S1). This isolate also harbours an R20fs mutation in ptr1 (BPSS0039; folM), which 261 encodes a pteridine reductase that is involved in TMP/SMX resistance (Podnecky et al. 2017); 262 this frameshift truncates Ptr1 from 267 to 91 residues. Both strains isolated from CF11 are also 263 highly resistant to TMP/SMX (MIC ≥32 µg/mL) and encode the BpeS missense variants V40I increased by 4.5-fold (22x) in the latter CF6 strain. Six proximal genes (BPSS0945; BPSS0948-281 BPSS0952) were also upregulated by 3.1-to 4.7-fold (9 to 26x; Table S1). One of these, 282 BPSS0945, is a peptidase and a putative virulence factor that may play a role in multinucleated 283 giant cell formation (Singh et al. 2013). 284 A gene duplication event encompassing penA has also been documented in the CF11 isolates. 285 The initial strain showed an elevated MIC towards CAZ (12 µg/mL), corresponding with a ~10x 286 duplication of a 36.7kb region that includes penA, whereas the latter strain had a 2x duplication 287 of this region and a wild-type CAZ MIC (2 µg/mL) (Viberg et al. 2017). As expected, penA was 288 downregulated by 2.1-fold (4x) in the latter isolate due to five times fewer copies of this gene. 289 Downregulation of other genes within the 36.7kb locus ranged from 1.4 to 3.3-fold (3 to 10x; 290 Table S1).  Both CF11 strains encode a large deletion in amrR (amrR ΔV62-H223 ). This mutation results in a 301 2.6-to 3.2-fold (6 to 9x) upregulation of amrAB-oprA in these isolates. In addition, a previously 302 undocumented S130L mutation in BPSL3085 likely contributes to the decreased susceptibility 303 observed towards doxycycline. 304 Evidence of convergent DE between early and latter CF isolates. Finally, we performed a 305 comparison of expression profiles from all CF cases to identify a signal of convergent gene 306 expression (pathoadaptation) across early vs latter isolates. To yield the most robust and 307 relevant analysis, we excluded the latter isolate from CF10 due to a lack of DE in this strain, and 308 the initial isolate from CF11, which was retrieved >3 years after infection and had already 309 undergone substantial genetic and transcriptional modifications. Exclusion of both strains was 310 supported by a lack of convergent signal when they were included in the analysis (data not 311 shown). Using these parameters, 17 genes were found to be significantly upregulated, and two 312 were significantly downregulated (Table S2). Five (26%) loci encode for hypothetical proteins 313 with no known function, of which four were upregulated. Among the convergently upregulated 314 genes with known function was the RND efflux pump BpeEF-OprC (4.8-to 6.1-fold; 28 to 69x), 315 the CydAB cytochrome bd quinol oxidase (5.5-to 5.9-fold; 45 to 60x), and the quorum sensing The causative agent of the tropical disease melioidosis, B. pseudomallei, is an uncommon 321 pathogen in CF, with fewer than 30 cases documented worldwide to date (Geake et al. 2015). 322 We have recently performed comparative genomic analysis of isogenic strains collected 323 between 4 and 55 months apart from the airways of seven of these cases (Viberg et al. 2017). 324 Here, we sought to further characterize these chronic cases by examining the transcriptomes of 325 five paired B. pseudomallei isolates retrieved between 10 and 55 months apart. Isolates were 326 cultured in an artificial CF sputum medium (Fung et al. 2010) to mimic their original in vivo 327

environment. 328
Under these conditions, DE was detected in four of the five cases and ranged from 32 to 792 329 genes, with the hypermutator strain from CF9 contributing the greatest number of DE loci (Table  330   S1). Interestingly, when compared with the number of genetic changes occurring in each isolate 331 pair, the latter isolate from CF11 had a higher proportion of DE loci to mutations (11.3) than CF9 332 (7.1), demonstrating that hypermutation does not necessarily lead to a similarly high number of 333 transcriptional differences. The one case with no DE, CF10, exhibited no genetic changes (i.e. 334 SNPs, small indels, copy-number variants, or large deletions) and had the shortest time 335 between isolate collection at 10 months (Viberg et al. 2017). All other cases encoded genetic 336 differences between pairs. The DE genes fell into several functional categories (Table S1) gmhA-wcbL-wcbK in both CF11 isolates. In both cases, it is likely that CPS-I production was 380 either substantially reduced or abolished. Although CPS-III is not required for virulence, it is 381 noteworthy that this locus was also downregulated in CF9 and CF11, as it has been previously 382 shown that CPS-III expression is tied to that of CPS-I genes (Ooi et al. 2013). In addition to TMP/SMX resistance, the initial isolate from CF11 is resistant to CAZ (12 µg/mL) 422 and the latter isolate from CF6 is highly resistant to CAZ (≥256 µg/mL). Our prior genomic study 423 showed that CAZ resistance in the initial CF11 strain was due to a 10x duplication of a 36.7kb 424 region encompassing the β-lactamase gene, penA, the first time that gene duplication has been 425 shown to confer CAZ resistance in B. pseudomallei (Viberg et al. 2017). In contrast, a 2x 426 duplication of this region in the latter strain did not raise the CAZ MIC above wild-type levels 427 transcriptional changes provides a potential means to predict pathogen behavior and evolution 442 across multiple CF cases in a relatively straightforward manner. Such predictability could 443 conceivably be exploited to improve the diagnosis or treatment of intractable CF infections, or 444 ideally, to prevent them from progressing in the first place. Therefore, a major objective of this 445 study was to identify evidence of convergence in B. pseudomallei gene expression during its 446 transition to a chronic infection. Despite the small number of CF melioidosis patients available 447 for this study, a signal of convergent pathoadaptation was identified between the initial and latter 448 isolates, with 19 significantly DE loci identified, 17 of which were upregulated (Table S2). This 449 convergence is noteworthy given the large size of the B. pseudomallei genome and the many 450 redundant pathways that could lead to similar adaptive phenotypes, a phenomenon that is well-451 recognized in P. aeruginosa (Marvig et al. 2015). One advantage of identifying convergence 452 using transcriptomics rather than genomic data is that it can reveal the transcriptional 453 consequence of multiple genetic mutations; for example, we have observed that multiple high-dose antibiotic therapy (Viberg et al. 2017). It was therefore not surprising to identify the 462 convergent upregulation of bpeEF-oprC (4.8-to 6.1-fold; 28 to 69x), which was significantly 463 upregulated in two of the four patients with DE (CF6 and CF11), and which led to TMP/SMX 464 resistance as discussed above. The second convergently upregulated locus was the cydAB 465 operon (BPSL0501 and BPSL0502), which encodes for cytochrome bd quinol oxidase (5.5-to 466 5.9-fold; 45 to 60x); this locus was significantly DE in CF6 and CF9. CydAB is an aerobic 467 terminal oxidase that oxidizes ubiquinol-8 and reduces oxygen to water under oxygen-limiting 468 conditions. This enzyme is better able to scavenge oxygen under microaerobic conditions 469 compared with cytochrome o oxidase, which otherwise predominates as the terminal respiratory 470 enzyme in electron transport-associated energy production (Cotter et al. 1997). Voggu and 471 colleagues demonstrated that the cydAB loci encoded by non-pathogenic Staphylococcus 472 species were better able to resist P. aeruginosa antagonism in the CF lung compared with the  Of the two convergently downregulated loci, only one, BPSS0351, has an assigned function, 515 although little is known about the role of this gene and its product in B. pseudomallei. This gene 516 encodes CueR (3.4-fold; 11x), a MerR family copper response regulator that is highly sensitive 517 to the presence of copper (Cu) and which regulates the transcription of genes that protect 518 against toxic metal ion concentrations (Brown et al. 2003;Singh et al. 2004). Cu has a long 519 history as an effective antimicrobial agent due its ability to generate reactive oxygen species, 520 with Cu accumulation in the mammalian host purported to act as an innate immune defense 521 mechanism to restrict pathogen growth (Samanovic et al. 2012). Thus, downregulation of cueR 522 in the latter CF isolates may represent a mechanism for mitigating Cu toxicity in the host, 523 similarly to E. coli (Singh et al. 2004). However, there are contradictory reports as to whether Cu pers. comm.). Final concentrations of the components were: 10 g/L mucin, 1.39 g/L salmon 579 sperm DNA (Sigma-Aldrich), 5 g/L NaCl, 2.2 g/L KCl, 0.22 g/L CaCl2, 5 g/L casein acid 580 hydrosylate (Sigma-Aldrich), 10 g/L bovine serum albumin (Roche Diagnostics, Castle Hill, 581 NSW, Australia), 0.005% diethylenetriaminepentaacetic acid, and 0.5% of egg yolk emulsion. 582 Each batch was tested for sterility prior to use by plating 100 µL onto Luria Bertani (LB) agar 583 (Oxoid) and incubating aerobically for 24 h. pH was tested using an aliquot of the medium to 584 ensure it was within the desired range (pH ~6.5-7). The medium was stored at 4 o C for no 585 longer than four weeks prior to use. 586 Viability counts. Two sets of viability counts were performed for this study. The first was 587 conducted to determine the number of colony-forming units (cfu) at OD590, which enabled us to 588 standardize the starting number of cells inoculated into the artificial sputum medium. The 589 second was conducted to verify the final concentration of cells across all CF isolates, which 590 enabled us to determine the number of cells for nucleic acid extraction to ensure that 591 approximately equal cell amounts were processed for each pair. The CF isolates were 592 subcultured from glycerol stocks onto LB agar at 37 o C for 24 h. Cells were suspended into 593 phosphate-buffered saline (PBS) followed by spectrophotometric measurement at OD590=1.0 in 594 a WPA CO 8000 cell density meter (Biochrom Ltd, Blackburn, VIC, Australia). Tenfold dilutions 595 and plating of cultures onto LB agar was carried out, followed by enumeration at 24 h. Viable 596 counts demonstrated that all CF isolates exhibited similar cell density when normalized to an 597 OD590=1.0 (range 1.3x10 8 to 4.9x10 8 ). Based on these counts, the starting amount of culture for 598 the CF isolates was standardized to 10 5 cfu for all subsequent experiments. 599 Growth curves in the artificial sputum medium. To minimize laboratory passage, each 600 culture was again subcultured from the original glycerol stocks onto LB agar at 37 o C for 24 h, 601 followed by a replication of OD590 measurements as determined previously. Based on the 602 viability count data, samples were then diluted to 10 6 cfu/mL in PBS. One hundred µL of this 603 suspension (~10 5 cfu) was used to inoculate 1.9 mL of the sputum medium, which was 604 aliquoted into 14 mL Nunc round-bottom culture tubes (Thermo Fisher Scientific, Scoresby, VIC, 605 Australia). Due to biosafety concerns, cells were grown in closed-capped tubes, and were 606  hands, we found this DNase I treatment to be insufficient for removing all contaminating DNA. 632 Extractions for RNA-seq were therefore further treated with TURBO DNA-free kit (Ambion, 633 Scoresby, VIC, Australia). For each sample, 35 µL of extracted RNA was incubated with 6 U 634 TURBO DNase at 37 o C for 32 min. The remaining RNA was not treated with this second round 635 of DNase; instead, this sample was used as template for PCR contamination screening, as 636 described below. All samples were transferred to clean RNase/DNase-free tubes for 637 downstream processing. 638 RNA quality control. To verify the removal of DNA from the total RNA extractions, two 639 contamination screens were performed. The first was used to determine the removal of salmon isolates vs. all latter CF isolates, with the latter CF10 isolate excluded due to a lack of DE in this 685 strain and the initial CF11 isolate excluded due to >3 years of infection prior to its isolation 686 (Viberg et al. 2017). For all analyses, DE was defined as a log2 fold change of ≥1.5 and a false 687 discovery rate of ≤0.01. To improve visualization of DE loci in the volcano plot of the initial and 688 latter comparison (Figure 1), highly expressed DE loci in only a single strain were omitted. 689

DATA ACCESS 690
The RNA-seq data generated in this study are available on the Sequence Read Archive 691 database under BioProject number PRJNA398168 and submission numbers SRR6031143 to 692

DISCLOSURE DECLARATION 712
The authors declare no conflicts of interest. 713