Genomic and Phenotypic Evolution of Achromobacter xylosoxidans during Chronic Airway Infections of Patients with Cystic Fibrosis

ABSTRACT Bacterial pathogens evolve during chronic colonization of the human host by selection for pathoadaptive mutations. One of the emerging and understudied bacterial species causing chronic airway infections in patients with cystic fibrosis (CF) is Achromobacter xylosoxidans. It can establish chronic infections in patients with CF, but the genetic and phenotypic changes associated with adaptation during these infections are not completely understood. In this study, we analyzed the whole-genome sequences of 55 clinical A. xylosoxidans isolates longitudinally collected from the sputum of 6 patients with CF. Four genes encoding regulatory proteins and two intergenic regions showed convergent evolution, likely driven by positive selection for pathoadaptive mutations, across the different clones of A. xylosoxidans. Most of the evolved isolates had lower swimming motility and were resistant to multiple classes of antibiotics, while fewer of the evolved isolates had slower growth or higher biofilm production than the first isolates. Using a genome-wide association study method, we identified several putative genetic determinants of biofilm formation, motility and β-lactam resistance in this pathogen. With respect to antibiotic resistance, we discovered that a combination of mutations in pathoadaptive genes (phoQ and bigR) and two other genes encoding regulatory proteins (spoT and cpxA) were associated with increased resistance to meropenem and ceftazidime. Altogether, our results suggest that genetic changes within regulatory loci facilitate within-host adaptation of A. xylosoxidans and the emergence of adaptive phenotypes, such as antibiotic resistance or biofilm formation. IMPORTANCE A thorough understanding of bacterial pathogen adaptation is essential for the treatment of chronic bacterial infections. One unique challenge in the analysis and interpretation of genomics data is identifying the functional impact of mutations accumulated in the bacterial genome during colonization in the human host. Here, we investigated the genomic and phenotypic evolution of A. xylosoxidans in chronic airway infections of patients with CF and identified several mutations associated with the phenotypic evolution of this pathogen using genome-wide associations. Identification of phenotypes under positive selection and the associated mutations can enlighten the adaptive processes of this emerging pathogen in human infections and pave the way for novel therapeutic interventions.


RESULTS
Clinical isolates of A. xylosoxidans infecting patients with CF. A total of six patients with CF, chronically colonized by A. xylosoxidans, were identified between 2011 and 2019 at the CF center in Lund, Sweden (see Table S1A in the supplemental material). To investigate the genetic adaptation of A. xylosoxidans in the CF host airways and discover mutations responsible for altered phenotypes, we sequenced the genomes of 55 isolates collected longitudinally from these patients (median, 9 isolates per patient; range, 3 to 17) (Fig. 1). The sequenced isolates were collected over a median period of 5.0 years (range, 2.3 to 7.2). The patients' median age was 19.9 years at first isolation of A. xylosoxidans (range, 12.1 to 36.2). For all but one patient, the first isolate included in the study was collected within 6 months from the initial detection of A. xylosoxidans in the sputum. For patient CF1, no isolate was available within the initial 6.4 years of colonization. Colonization by A. xylosoxidans continued for all but two patients (CF4 and CF5) throughout the study period. CF4 cleared A. xylosoxidans after 2.3 years, while A. xylosoxidans colonization ended in CF5 after a double lung transplant in September 2018.
Hypermutator phenotype increases within-host genetic diversification of A. xylosoxidans. Initial analysis of multilocus sequence typing (MLST) using de novo assembled genomes showed that each of the 6 patients were infected by an individual unique clone throughout the study period (Table S1C). We therefore assigned the same nomenclature (CF1 to CF6) for clones from each patient. To gain insight into A. xylosoxidans within-host evolution in CF airways, we identified genetic variants that differentiated isolates of the same clone; i.e., we identified mutations that had accumulated during the period of infection in each of the six clones. In total, we identified 1,499 variants (median, 76; range, 12 to 931), including 837 nonsynonymous (NS) single-nucleotide polymorphism (SNPs), 165 insertions/deletions (indels), 124 intergenic SNPs, 110 intergenic indels, and 263 synonymous (S) SNPs (Table S1D). Ten of the 55 isolates had a large genetic distance (.100 SNPs) from the initial isolate of the same patient. One of the 10 isolates belonged to CF3, and the other nine isolates all belonged to CF1 ( Fig. 2A). The CF3-S12a isolate, found in the latest sputum sample from this patient, contained a nonsense mutation in the DNA mismatch repair system (MMR) gene mutL, which was absent in all other isolates of this clone and which could explain the elevated mutation number of this isolate. Moreover, a manual inspection of the MMR system genes in the 12 isolates of CF1 showed that all nine isolates with higher genetic distance had a 112-bp deletion in mutL. The deletion was found in all of the following genomes from the third isolate, but the fifth isolate, which also had the deletion, was below the threshold of large genetic distance (64 SNPs relative to the first isolate) ( Fig. 2A). As this deletion leads to a functional defect in MutL, it is likely responsible for the increased number of mutations in this clone. We also found higher levels of transitions to transversion (Ts/Tv) among SNPs of CF1 (Ts/Tv = 29) and CF3 (Ts/ Tv = 8) than among SNPs of other clones (Ts/Tv % 1) (Fig. 2B), which is in agreement with the expected effect of the MutL defects (23).
In the BEAST analysis of the nucleotide substitution rate of A. xylosoxidans clones infecting each patient, we found a higher estimated mean substitution rate of 19.8 SNPs/year for isolates of CF1 relative to the other clones (range, 0.9 to 5.9 SNPs/year) ( Table S1E). The higher substitution rate for CF1 was consistent with the hypermutator phenotype. In contrast, the estimated mean substitution rate for isolates of CF2 was 0.9 SNPs/year, which was lower than that of all other clones. The overall average substitution rate for all clones excluding the hypermutator CF1 was 2.7 SNPs/year (Fig. 2B), which is close to the previously reported 1.9 SNPs/year in CF adapted isolates of A. xylosoxidans genomes (24). In comparison, the mean substitution rates of CF adapted genomes of P. aeruginosa and B. multivorans were reported at 2.6 and 2.4 SNPs/year, respectively (22,25).
A. xylosoxidans clones have varying evolutionary dynamics. To identify evidence for either positive or negative selection of mutations in isolates of each clone, we calculated the relative rates of NS and S substitutions (dN/dS) (Fig. 2C). While mutations in isolates of CF1, CF5, and CF6 occurred mostly under neutral selection, dN/dS was higher for CF2 and CF4. However, we found a significant increase of dN/dS only in CF2 compared to other clones (Fisher's exact test, P = 0.04), possibly explainable by the lower total number of substitutions in CF4. As isolates of CF2 also had a lower substitution rate, it is clear that the majority of substitutions in this clone are NS. In contrast, isolates of CF3 exhibit a sign of negative selection with a dN/dS of 0.8, and the dN/dS ratio is significantly lower in this clone than the other clones (Fisher's exact test, P = 0.01). Interestingly, in-depth analysis of the selection signature in CF3 showed two significantly different phases of selection for the earlier eight isolates (dN/dS = 0.5) compared to the later nine isolates (dN/dS = 1.1) (Fig. S1) (Fisher's exact test, P = 0.01). However, one major caveat of this evolutionary dynamics analysis is that the calculated ratio of substitution rates (dN/dS) was performed for clones from a single population;  (Table S1B). All sequenced clinical isolates of A. xylosoxidans are represented by black dots. The colored lines show colonization of A. xylosoxidans before the first or after the last sequenced isolates. The years 2006 to 2010 are not shown, as no isolates from these years were sequenced.
A. xylosoxidans Genetic Adaptation in Patients with CF thus, the strength of positive selection is likely underestimated by the presence of polymorphisms that are not yet fixed in the population (26).
Independent clones of A. xylosoxidans show convergent evolution in pathoadaptive loci. To explore the function of genes enriched for mutations in CF-adapted A. xylosoxidans, we analyzed the distribution of NS mutations within 21 different clusters of orthologous groups (COG) of proteins found in A. xylosoxidans genomes (Table S1F). When NS mutations from all clones were analyzed (Fig. 3A), we found a significant enrichment of the group "Intracellular trafficking, secretion and vesicular transport" (U) among NS mutations compared with the presence of this group in the reference genome (Fisher's exact test, P = 0.001). Surprisingly, COG U was overrepresented by the contribution of CF1, since when NS mutations of this clone were excluded from the analysis, we found a significant enrichment of the group "Signal transduction mechanisms" (T) (Fisher's exact test, P = 0.001). We also observed some divergence in other groups targeted by NS mutations from each clone individually. For example, the group "Transcription" (K) was enriched in targets of CF2, CF5, and CF6 while being unchanged or even depleted in those of CF3 and CF1.
Despite these divergences in NS mutational targets, we still found consistent patterns of convergent evolution in certain genetic loci across independent clones. Using an intuitive approach based on length of the loci and the expected mutation number with respect to the normal distribution of mutations across the genome, we found four genes and two intergenic regions to be significantly enriched by parallel independent mutations across different clones (Fig. 3B). As these six loci are likely to be involved in pathogen host adaptation, we refer to them as pathoadaptive loci. All four genes were mutated only by NS mutations. The most frequently mutated gene across isolates of four clones was a transcriptional regulator gene with a product putatively identified as the biofilm-associated growth repressor BigR. As suggested by the annotation, this protein has been associated with regulation of biofilm growth in response to hydrogen sulfite in Xylella fastidiosa and Acinetobacter baumannii (27,28). Two of  (Table S3). (B) Summary data from each clone on the presence of deleterious mutations in DNA MMR gene mutL leading to the hypermutator phenotype, the fold change of Ts/Tv (transition/transversion), and the mean substitution rate estimated with BEAST. (C) Nonsynonymous-to-synonymous mutation ratio (dN/dS) for each clone. The red dotted line at a dN/dS of 1 represents neutral evolution. The asterisk indicates significance of change in number of N over S mutations in each clone compared to those found in the entire data set excluding that clone (Fisher's exact test, P , 0.05). the other pathoadaptive genes, encoding a LysR-family transcriptional regulator and the heat-inducible transcriptional repressor HrcA, were also transcriptional regulator genes. The fourth pathoadaptive gene encoded a putative HAMP domain histidine kinase PhoQ, which is associated with virulence, motility, and polymyxin resistance in Salmonella enterica serovar Typhimurium and P. aeruginosa (29)(30)(31). With the exception of hrcA, all other genes were found to be pathoadaptive in CF-adapted A. xylosoxidans from another study, confirming the relevance of these loci in genetic adaptation of A. xylosoxidans in CF (24). Moreover, the two intergenic loci had mutations upstream of genes encoding the septum formation protein Maf and a TonB siderophore receptor, which may be involved in metal uptake through siderophores similar to metal uptake in P. aeruginosa (32,33).
Evolution of adaptive phenotypes in independent clones of A. xylosoxidans. To understand the phenotypic evolution of A. xylosoxidans in CF, we tested the 55 sequenced isolates with some of the known adaptive phenotypes from other bacterial species found in the CF model (22,34).
In our study, the analysis of swimming motility showed that this phenotype is gradually lost over time following the initial colonization (Fig. 4). While evolved isolates of CF2, CF5, and CF6 had a rapid decline in motility over longer colonization time, decrease in motility was less pronounced for evolved isolates of CF1 and CF3. This may  (Table S1G). The presence of mutation in isolates of each clone is shown by the black square. US, upstream of intergenic loci. be because the first isolates of CF1 and CF3 already had lower swimming motility than first isolates from other patients. Moreover, unlike other clones, all three isolates of CF4 had an unchanged swimming motility. Looking at the possible mutations for the observed phenotype (Table S2), we found an indel in the gene encoding the flagellar protein FliS in isolate CF1-S4 that may be responsible for the diminished swimming motility of this isolate. We also identified NS mutations in the two genes encoding FimV and a fimbrial usher protein (AT699_RS04965) in isolate CF1-S9c, which had no swimming motility (Fig. S2).
Previous studies have shown that P. aeruginosa isolates converge toward a lower growth rate following adaptation in the CF lung (34). We also found some of the evolved isolates of A. xylosoxidans from CF2, CF4, CF5, and CF6 to be slower growing in the ABT minimal medium with glucose and Casamino Acids (Fig. 4). In contrast, most evolved isolates of CF1 and CF3 were relatively unchanged over time with respect to growth in minimal medium; likely because the first isolates of these two clones already had a lower growth rate than isolates of the other clones.
To dissect the genetic background behind development of these three phenotypes in A. xylosoxidans, we used a k-mer-based GWAS method (35). Utilizing this method, we identified three genetic components as being significantly associated with the motility and 29 components with the biofilm phenotypes (q , 0.05) (Table S1K to L). Some of the identified components were in genes with putative homologs that have previously been linked to these phenotypes in other species of bacteria, demonstrating their probable relevance for the observed phenotypes of A. xylosoxidans colonizing the CF host (Table 1).
Increased resistance to meropenem and ceftazidime is established by the direct or indirect regulatory effect of various mutations. All sequenced isolates of A. xylosoxidans from our collection were resistant to several classes of antibiotics. Based on the clinical data of antimicrobial susceptibility testing (Table S1M), 100% (n = 48) of the tested isolates were resistant to aztreonam, 98% (n = 54) were resistant to ciprofloxacin, 88% (n = 44) were resistant to gentamicin, 87% (n = 48) were resistant to tobramycin, 69% (n = 37) were resistant to colistin, 47% (n = 26) were resistant to ceftazidime, and 33% (n = 18) were resistant to meropenem. Furthermore, we found 49 isolates to be multidrug resistant (MDR), i.e., nonsusceptible to at least one agent of three or more classes of antibiotics. As most of the initial isolates were susceptible to ceftazidime and meropenem, we chose to analyze the subsequent isolates from each patient for  (Table S1H and I). The values are plotted against the number of years the isolate colonized the patient since the first isolate (Table S1B) to follow the change of the phenotype over colonization time in CF. The colored lines are the best-fit linear regression for isolates of each clone (colored symbols). resistance to these two antibiotics. Notably, there was a consistent increase in resistance to ceftazidime and meropenem in isolates from CF3 and CF5, which was probably selected in response to the regular treatment of these two patients with different b-lactam antibiotics like ceftazidime, meropenem, imipenem, and piperacillin-tazobactam (Fig. 5).
By using GWAS, we found a total of 28 genetic components to be linked to ceftazidime resistance and 5 components associated with meropenem resistance (Table S1N). Most of the components identified represented changes in the genome of isolates from CF3, which may be a result of the higher number of resistant and susceptible a The estimated effect of k-mer presence on the phenotype, the reference for the reported effect, and the annotation of each gene are reported. The full list available in Table S1K to L.  (Table S1M); the dotted lines represent the considered breakpoints of resistance, and the three bar colors represent different patterns of susceptibility to ceftazidime and/or meropenem. The middle panels represent the timeline of isolates from these two patients and the administered antibiotic combinations based on the clinical data. All b-lactams are in red. The question mark above MEM means that meropenem was used by patient CF3 via inhalation for an unknown period of time. The lower panels represent all genetic changes that may contribute to the increased resistance of ceftazidime and meropenem. Apart from the pathoadaptive gene hrcA, all other changes are in putative homologs of previously characterized components of b-lactam resistance in other bacterial species. The presence of the genetic change in each isolate from the topmost panel is denoted by colored squares.
A. xylosoxidans Genetic Adaptation in Patients with CF isolates from this patient. Moreover, a large number of a total of 33 identified components were in genes encoding signaling proteins (COG T; n = 5) or transcription proteins (COG K; n = 7) . Surprisingly, we also found three of the identified pathoadaptive genes (phoQ, hrcA, and bigR) to be present among the components associated with ceftazidime resistance, suggesting that these regulatory loci are vital for adaptive phenotypes such as antibiotic resistance (Fig. 5). The pathoadaptive gene encoding the putative histidine kinase PhoQ was associated with ceftazidime resistance by two independent genetic components. Mutations in homologs of this gene have an established role for resistance to polymyxin and aminoglycosides in P. aeruginosa (36), but they have also been selected in ceftazidime-treated evolved cultures of S. maltophilia (37). The pathoadaptive gene encoding the biofilm growth-associated repressor BigR has also been shown to have an indirect regulatory effect on antibiotic resistance. The putative homolog of this protein known as ArsR in Enterococcus faecium has been involved in positive regulation of a penicillin-binding protein (PBP), and its deletion has increased resistance to b-lactam antibiotics (38). No information could be obtained for the putative role of hrcA homologs in antibiotic resistance, but two other genetic components were found in regulatory encoding genes related to b-lactam resistance (Fig. 5). The first was a gene encoding a putative histidine kinase, CpxA, which has a well-established profile with respect to multidrug antibiotic resistance in several pathogenic bacteria. In this context, the CpxAR two-component response regulator is involved in regulation of the RND (resistance-nodulation-division)-type efflux pumps controlling resistance to b-lactam, aminoglycosides and certain antimicrobial peptides in Escherichia coli, P. aeruginosa, S. Typhimurium, Klebsiella pneumoniae, and Vibrio cholerae (39)(40)(41)(42)(43)(44). The other component was in the gene encoding the stringent response (p)ppGpp synthase/hydrolase SpoT. Among other functions, (p)ppGpp modulates the expression of PBP or efflux pumps and is linked to b-lactam resistance in other bacterial pathogens (45). While GWAS may be applicable for unbiased discovery of components linked to antibiotic resistance, it has a statistical limitation for detection of less frequent genetic changes. We therefore looked for other biologically relevant mutations to explain the evolution of b-lactam resistance and found mutations in two genes encoding putative direct effectors of b-lactam resistance in A. xylosoxidans (Fig. 5). In this context, the hypermutator isolate CF3-S12a, which lacked most of the other discussed mutations identified by GWAS, harbored two NS mutations in the gene encoding the penicillin-binding protein FtsI, which is associated with ceftazidime resistance in E. coli (46). Furthermore, in CF5, where a consistent development of resistance in both ceftazidime and meropenem occurred from the third isolate onward, a nonsense mutation was fixed in ampD, encoding N-acetyl-anhydromuramyl-L-alanine amidase. Inactivation of AmpD has been linked to constitutive overproduction of AmpC b-lactamase and to the rise of b-lactam resistance in most Gram-negative bacterial pathogens (47). Interestingly, we also found another fixed mutation in the gene encoding the pathoadaptive putative histidine kinase PhoQ in the third isolate and all later isolates.

DISCUSSION
Infections of A. xylosoxidans have recently become more prevalent in patients with CF (5,6,9), but there are few studies that characterize the genetic adaptation of Achromobacter spp. in CF infections (24,48,49). To complement these studies, we used a combination of functional and genomic approaches to reflect on the biology of A. xylosoxidans adaptation in the CF model. In this context, our data set included a larger number of genome-sequenced isolates from each patient to provide a higher resolution of the A. xylosoxidans within-host diversity and capture possible coexisting clades in the different compartments of the CF airway (50). Moreover, we focused exclusively on A. xylosoxidans, because the genetic and phenotypic diversity of different Achromobacter species may cause unwanted neglect of A. xylosoxidans specific traits (51,52). Finally, in addition to the genetic changes, we also investigated the phenotypic changes of A. xylosoxidans isolates and tried to link the genotype with the observed phenotypes. Having multiple isolates per each patient provided an advantage in the analysis of GWAS, where several biologically relevant genetic components were identified for biofilm formation, swimming motility, and b-lactam resistance.
The cystic fibrosis airway represents a complex biological system with multiple environmental stressors affecting adaptation of the colonizing microorganism (53). The differences in composition of infecting pathogens, antibiotic treatments, and immunological host factors impose various selection pressures on adaptation of A. xylosoxidans. This may account for some of the heterogeneity in the genomic characteristics of the six clones isolated from the patients in this study. The hypermutator phenotype was detected in 20% of the isolates (n = 11). Hypermutators have been found previously in A. xylosoxidans and other species isolated from CF infection, and it is assumed that the higher mutation rate has a positive effect on faster selection of adaptive phenotypes in this environment (14,25,49,54,55). It is also likely that this phenotype is fixed in CF1 as a result of its longer colonization in the patient, since it has been suggested that longer evolution history is associated with the hypermutator phenotype (56). Additionally, we observed different signatures of selection between the six clones of A. xylosoxidans that even shifted in one clone (CF3) at different phases of early and late colonization. This demonstrates that A. xylosoxidans has a dynamic adaptation and thrives by natural selection in different infection backgrounds. In this context, negative selection has previously been observed in CF isolates of P. aeruginosa following an initial period of rapid positive selection (57). However, despite these heterogeneities, we still found clear patterns of genetic adaptation of A. xylosoxidans during airway infections in CF.
Our analysis revealed that evolution of regulatory loci plays a major role in withinhost evolution of A. xylosoxidans during CF infection. Regulatory proteins often entail pleiotropic functions and are therefore expected to have substantial effects on the overall phenotype development of a bacterial pathogen during infection. We showed that genes encoding signal transduction proteins are the most frequent targets of NS mutations in most of the CF adapted A. xylosoxidans clones. Furthermore, genes encoding other regulatory proteins belonging to the transcription COG category (K) were also under positive selection in at least half of the clones. A possible way to interpret negative or neutral selection of the transcription genes in CF1 and CF3 may be that most of the adaptive NS mutations within these genes are already fixed in the first isolates of these clones as a result of earlier evolution in the CF airway. Previous studies have also shown that mutations in regulatory genes are fixed in the earlier adaptation of P. aeruginosa during CF infections (57). With regard to our patients, we know that CF1 was colonized with A. xylosoxidans for more than 6 years prior to the first isolate included in the study, but to our knowledge there were no previous isolates for CF3. It is possible that this patient had an undetected colonization for some time before the first positive culture, or that the patient was colonized by an evolved CF clone through patient-to-patient transmission (14). However, such a transmission is unlikely to come from patients within the same CF center, as all patients with chronic A. xylosoxidans infection in this study carried their own individual clones. One of the probable transcription genes that could have had preexisting mutations in CF1 is axyZ, which was also mutated in isolates of CF2 and CF5. This gene was shown to be under positive selection in other CF-adapted isolates of Achromobacter spp. and encodes a transcriptional regulator which represses the AxyXY-OprZ multidrug efflux system in Achromobacter spp. (24,58).
In addition to the enrichment in targets of NS mutations, we showed that all of the pathoadaptive genes in A. xylosoxidans encode regulatory proteins. This finding is also comparable to the enrichment of transcriptional regulators among targets of pathoadaptive genes in CF-adapted P. aeruginosa (18). Moreover, for the first time in studies of Achromobacter sp. CF adaptation, we identified two regulatory intergenic loci under positive selection across multiple independent clones. Considering the location of the mutations from these regions, they may regulate expression of the adaptive downstream genes. Previous functional studies emphasized the relevance of these mutations in regulation of adaptive bacterial phenotypes (59,60).
Finally, with regard to the b-lactam resistance phenotype, several of the components identified by the GWAS analysis were also in genes encoding regulatory proteins. As mentioned in Results, regulatory proteins with established profiles in resistance of other bacterial species, such as SpoT, BigR, CpxA, and PhoQ, are not only regulators of b-lactam resistance. On the contrary, most of these proteins operate within the broader pleiotropic resistance systems by, for example, regulating the transcription of multidrug efflux pumps. In a recent experimental evolution study, it was shown that mutations in genes involved in efflux pump expression were selected in P. aeruginosa in response to combination antibiotic therapy (60). As many of the patients included in our study also received a combination of several different classes of antibiotics, selection of regulatory proteins with pleiotropic resistance effects might be beneficial over selection of direct resistance factor genes. In support of this hypothesis, the last six isolates of CF3, which were resistant to both meropenem and ceftazidime, were also consistently resistant to colistin. Therefore, it is possible that the regulatory mutations in these isolates, among other possible mutations affecting colistin resistance, also had a positive effect on colistin resistance and were selected as a result of antibiotic combination therapy. Beyond the detection of mutations in regulatory protein-encoding genes, finding ampD or ftsI mutations in association with b-lactam resistance in A. xylosoxidans may seem more trivial. Nevertheless, it highlights the idea that, similar to that in other successful pathogens (61), the evolution of antibiotic resistance in A. xylosoxidans is multifactorial and the product of independent direct or indirect regulatory pathways.
One unique advantage of the k-mer-based GWAS method (DBGWAS) (35) is that it uses the assembled genomes of each isolate and includes variants of all possible accessory genes or mobile genetic elements specifically present in some isolates or clones. Changes within these features are missed by the reference genome-based variant calling. This enabled us to obtain several identified components that were absent in the reference genome. Therefore, this method is not limited to single nucleotide polymorphisms and short insertions or deletions, and more complex genetic changes can be captured. Moreover, we demonstrated that this method associates bacterial genomic with phenotypic evolution for a large collection of mutations and for phenotypes other than antibiotic resistance. However, while the GWAS analysis identified possible causal mutations for the antibiotic resistance, motility, and biofilm phenotypes, it is important to consider it as a preliminary and not a validated result. GWAS was underpowered by the limited number of genomes in our study, and some possible variants with an effect on function may be filtered out because of occurrence in fewer isolates; i.e., we expect that some phenotype-related mutations were not captured by the analysis. We therefore tried to overcome these limitations by focusing only on components of genes with putative established functions in other bacteria. We also performed a manual inspection to find additional mutations not identified by GWAS but with potentially similar putative function based on other studies. To ascertain the real function of these candidate mutations, the mutant alleles have to be constructed in environmental or CF strains of A. xylosoxidans, and the isogenic isolates should be tested for the relevant phenotypes. However, one constraint with the current studies of the pathogen-host adaptation is that the phenotypes are difficult to examine in vitro owing to the differences between the infection niche and the in vitro conditions. Another caveat for a functional approach is that adaptive regulatory mutations may function in epistatic ways to cause significant effects on adaptive phenotypes (62,63). Therefore, studying single mutations in one genetic background may not always lead to the optimal result. Future studies considering these challenges can further elucidate the function of the adaptive regulatory mutations identified in this work. Nevertheless, such approaches should also consider the difficulties of generating bacterial genetics in A. xylosoxidans clinical isolates due to lower growth rates and unavailability of antibiotic resistance markers.
In conclusion, our investigation of the within-host evolution of A. xylosoxidans identified putative genetic determinants for the adaptation of this pathogen in the CF lung (Fig. 6). The analysis identified several regulatory loci under positive selection, and we discuss the likely effect of mutations for the development of the CF adaptive phenotypes. With special focus on the evolution of multidrug antibiotic resistance in this pathogen, we showed that A. xylosoxidans develops further resistance to the few antibiotics to which it was susceptible before colonization. Notably, we found multiple plausible paths for emergence of b-lactam resistance in A. xylosoxidans through mutations in direct or indirect regulatory factors, which may facilitate resistance to other classes of antibiotics used in for CF infection treatment. The knowledge of phenotypic and genetic adaptation of A. xylosoxidans can be vital for advancement of potential therapeutic targets and diagnosis of its infections in patients with CF.

MATERIALS AND METHODS
Patient cohort, bacterial strains, and antibiotic susceptibility testing. Patients with CF registered at the CF center at Skåne University Hospital in Lund and with colonization of A. xylosoxidans in the lower airways were eligible for inclusion in the study. The colonization was defined as chronic in patients where A. xylosoxidans was found in at least half of the sputum cultures for a minimum of 1 year. Clinical isolates from the study participants were obtained from Clinical Microbiology, Labmedicin Skåne, Lund, Sweden. Species identification was done according to standard laboratory methods, including matrix-assisted laser desorption ionization-time of flight mass spectrometry (MALDI-TOF MS) (Bruker Daltonics, Bremen, Germany) (64) and confirmed by sequencing of the nrdA locus as previously described (65). Antibiotic susceptibility testing was performed as part of the clinical routine at the department of clinical microbiology according to EUCAST guidelines. Since no breakpoints for resistance of Achromobacter spp. have been defined by EUCAST, we considered breakpoints for nonfermenting Gram-negative bacteria as previously described (66). Clinical patient data were obtained from medical records. The study was approved by the Swedish Ethical Review Authority (reference number 2019-02825).
Genome sequencing, assembly, annotation, and clone type definition. Genomic DNA was prepared from overnight cultures of bacterial isolates grown on LB medium using a DNeasy blood and tissue kit (Qiagen) and sequenced on an Illumina NextSeq 500 platform generating 150-bp paired-end reads after  (67) with default parameters and k-mer sizes ranging from 21 to 127. The average assembly size was 213 contigs (range, 171 to 270). The assembled contigs were annotated using Prokka (68) using the NCTC10807 reference genome (NCBI accession number PRJEB6403) as the first annotation priority. To analyze the genetic distance between isolates of the six patients, multilocus sequence typing (MLST) analysis was performed on the genome of each lineage using the MLST database and as previously described (65). Genome alignment, variant calling, and phylogenetic trees. Genome alignment and within-host introduced variants were called with BacDist (https://github.com/MigleSur/BacDist) (24,69), which is based on snippy variant calling (70), and filtering of variants present in all samples, using NCTC10807 as a reference genome with default parameters. Hypermutator phenotype was checked by inspection of the BacDist variants and manual search of additional deletion or insertions in the A. xylosoxidans DNA mismatch repair system genes as previously described (49). The within-host-introduced variants (Table S2) were used for phylogenetic tree generation with RAxML (71) using the GTRCAT model.
Substitution rate estimation and signatures of selection. The substitution rate for bacterial isolates was estimated for all lineages using BEAUti (BEAST version 2.5.0) (72). The Markov chain Monte Carlo (MCMC) was run for 50,000,000 iterations with the sequence alignments from BacDist as the input and using the HKY substitution model, strict clock parameters, coalescent constant population tree prior, 1/X prior for population size, and gamma prior for the clock rate. Tracer software (version 1.7.1) (73) was used to control convergence in an effective sample size and parameter value traces. Multiple tests were performed for each sample to confirm reproducibility and convergence. The calculated clock rate (per site per year) was multiplied by the total size of the alignment to get the substitution rate per genome per year. Signature of selection (dN/dS) was calculated by division of the total number of nonsynonymous mutations by the total number of synonymous mutations, which had been multiplied by three, since the number of nonsynonymous sites is approximately three times higher than the number of synonymous sites in other similar bacterial genomes (57). Statistical significance for higher or lower dN/dS for any clone was calculated by two-tailed Fisher's exact using the number of S and NS mutations of that clone compared to that of all remaining clones.
Analysis of the COG groups. The reference genome NCTC10807 and the list of variants (Table S2) were annotated using EGGNOG-mapper version 1.0.3 (74) at the DIAMOND, eggNOG's bacterial database, with a 1 Â 10 28 E value cutoff and 0.8 minimum query sequence coverage settings. Within each of the 21 identified COG targets, the number of NS mutations was counted and divided by the sum to get the percentage for each COG among NS mutations. This value was divided by the percentage for each COG in the NCTC10807 reference genome to obtain the enrichments of each group. Statistical significance for enrichment of each group in the data set was calculated by the two-tailed Fisher's exact test and adjustment of the raw P value using the Bonferroni method.
Identification of genetic loci under positive selection. The A. xylosoxidans reference genome NCTC10807 with the annotated coding sequence regions was retrieved from the NCBI database. The total observed number of mutations in each coding or noncoding locus was recorded for each clone (Table S4). To identify loci under positive selection across different clones (convergent evolution), each locus with more than one mutation in the same clone was considered to have one clonal mutation. The number of clonal NS mutations in each clone (Table S1D) was distributed across the 6,214 genes based on the length of each gene divided by the total size of the coding region, where longer genes are expected to acquire more mutations than shorter ones (Table S5). The sum of the expected value of mutations in each gene across the 6 clones was compared with that of the observed number of mutations. Genes mutated in at least half of the clones (3 clones), with an observed versus expected enrichment of more than 10-fold, and a statistically significant difference in mutation density [Poisson, P (x;m) , 0.0001, where x is the observed mutation number and m is the expected mutation number], were defined as genes under positive selection for adaptation, i.e., pathoadaptive genes. Similarly, for identification of pathoadaptive noncoding loci mutations, the same method was applied but with the additional constraint that the clonal mutations had to be within close proximity, i.e., no more than 35 bp apart, similar to a previously described method (59).
Swimming motility, growth rate, and biofilm formation assays. ABTGC minimal medium (75) consisting of 2 g/liter (NH4) 2 SO 4 , 6 g/liter Na 2 HPO 4 , 3 g/liter KH 2 PO 4 , 3 g/liter NaCl, 1 mM MgCl 2 , 0.1 mM CaCl 2 , 0.01 mM FeCl 3 , 2.5 mg/liter thiamine supplemented with 0.5% (wt/vol) glucose, and Casamino Acids (Merck), was used during all phenotype assays and to make overnight cultures. To measure swimming motility of A. xylosoxidans, a similar-sized bacterial colony of each isolate was spotted in the middle of an ABTGC 0.3% (wt/vol) agar plate. The plates were incubated at 37°C for 24 h, and the diameters of migration from the spot were recorded. The growth experiments were performed in 96-well microtiter plates (Thermo Scientific). Briefly, the overnight cultures from each isolate were diluted to an OD 600 of 0.1, and 15 ml was added in a total volume of 150 ml in three wells of the plates. Bacterial growth (OD 600 ) was measured continuously for 22 h at 37°C with shaking using a SpectraMax multimode reader (Molecular Devices). The growth data were analyzed with Microsoft Excel (Microsoft), where the doubling time of each strain was calculated during exponential growth.
Biofilm formation was also analyzed in 96-well microtiter plates using crystal violet staining. Each overnight culture was diluted to an OD 600 of 0.5 and inoculated by addition of 18 ml to a total volume of 180 ml in each well. Medium alone (180 ml) was used as the negative control. The plate was incubated for 48 h at 37°C and 125 rpm shaking. Thereafter, the bacterial broth was discarded, the wells were washed three times with 200 ml of phosphate-buffered saline (PBS), and the biomass in the wells was fixed by addition of 200 ml methanol for 10 min. After removal of the methanol, the wells were allowed to dry for 2 h and then stained with 160 ml of 0.1% (wt/vol) crystal violet for 4 min, followed by 3 washes with 200 ml PBS to remove the excess free crystal violet. The crystal violet was solubilized with acetoneethanol (4:1 [vol/vol]) and quantified at OD 550 using SpectraMax multimode reader. Each experiment was performed in triplicate and repeated by three biological replicates. For each measurement, one-way analysis of variance (ANOVA) followed by Tukey's honestly significant difference (HSD) multiple-comparison test using GraphPad Prism 9 was used to calculate the significance of change in the mean from the earliest isolates of each patient. The adjusted P value of ,0.05 was used for significance. Isolates with mean OD 550 values greater than three times that of the negative control were considered biofilm producers (76).
Identification of genetic loci associated with motility, biofilm, growth, and antibiotic resistance. The significant genetic components associated with motility, growth, biofilm, ceftazidime, and meropenem resistance were identified using the DBGWAS version 0.5.4 software (34). The cutoff for ceftazidime and meropenem resistance was defined as described above. The other cutoff values were defined as motility of .10 mm, doubling time of .100 min, and biofilm OD 550 of .0.3. The de novo assembled contigs of the 55 isolates were used as input in the software. All available annotations of Achromobacter genes from the UniProt database (77) were used for unitig annotations (271,851 genes; retrieved 17 April 2020) for motility, doubling time, and biofilm GWAS analysis; all known bacterial resistance genes from the UniProt database (76) were used for unitig annotations (7,792 genes; retrieved 17 April 2020) in ceftazidime and meropenem resistance GWAS. All identified components with q value of ,0.05 were investigated further. In every component with more than one significant node, the node with the lowest q value was considered. Every k-mer present in only one isolate was discarded to avoid any possible phenotype bias.
Data availability. DNA sequence reads from the 55 A. xylosoxidans isolates are available in the European Nucleotide Archive under study accession number PRJEB43175.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.

ACKNOWLEDGMENTS
We thank Gisela Hovold for excellent technical assistance and Muna Al-jammal at Clinical Microbiology, Labmedicin Skåne, Lund, Sweden, for retrieval of the clinical isolates. We thank Åke Forsberg at the MIMS, University of Umeå, and Mattias Collin at the BMC, Lund University, for their feedback and support. We acknowledge the Center for Translational Genomics at Lund University for whole-genome sequencing of the isolates. The computations and data handling/SIMILAR were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at SNIC CENTRE, partially funded by the Swedish Research Council through grant agreement no. 2018-05973.