Antagonistic Pleiotropy in the Bifunctional Surface Protein FadL (OmpP1) during Adaptation of Haemophilus influenzae to Chronic Lung Infection Associated with Chronic Obstructive Pulmonary Disease

Nontypeable Haemophilus influenzae is an important pathogen in patients with chronic obstructive pulmonary disease (COPD). To elucidate the bacterial pathways undergoing in vivo evolutionary adaptation, we compared bacterial genomes collected over time from 13 COPD patients and identified recurrent genetic changes arising in independent bacterial lineages colonizing different patients. Besides finding changes in phase-variable genes, we found recurrent loss-of-function mutations in the ompP1 (fadL) gene. We show that loss of OmpP1/FadL function reduces this bacterium’s ability to infect cells via the hCEACAM1 epithelial receptor but also increases its resistance to bactericidal fatty acids enriched within the COPD lung, suggesting a case of antagonistic pleiotropy that restricts ΔfadL strains’ niche. These results show how H. influenzae adapts to host-generated inflammatory mediators in the COPD airways.

pathways undergoing in vivo evolutionary adaptation, we compared bacterial genomes collected over time from 13 COPD patients and identified recurrent genetic changes arising in independent bacterial lineages colonizing different patients. Besides finding changes in phase-variable genes, we found recurrent loss-of-function mutations in the ompP1 (fadL) gene. We show that loss of OmpP1/FadL function reduces this bacterium's ability to infect cells via the hCEACAM1 epithelial receptor but also increases its resistance to bactericidal fatty acids enriched within the COPD lung, suggesting a case of antagonistic pleiotropy that restricts ΔfadL strains' niche. These results show how H. influenzae adapts to host-generated inflammatory mediators in the COPD airways. KEYWORDS CEACAM1, FadL/OmpP1, Haemophilus influenzae, adaptive evolution, antagonistic pleiotropic, chronic obstructive pulmonary disease, comparative genomics, convergent evolution, free fatty acids, in situ evolution B acterial pathogens evolve during long-term chronic infections (1), and identifying the genetic changes that arise in bacteria over time within hosts is a powerful means of deciphering host selective pressures acting on bacteria that colonize new niches and the specific adaptations that allow bacteria to persist (2)(3)(4). Dissecting the genetics of bacterial adaptation is also key to understanding and predicting the emergence of antimicrobial resistance and immune evasion traits within bacterial populations, as well as for identifying new interventional targets.
Colonization of the human lower airways by opportunistic pathogens of environmental origin leads to the irreversible progression of major respiratory diseases, such as cystic fibrosis (CF), and bacterial evolution during long-term lung infection has been extensively analyzed in this context (2,4,5). In particular, genomic analyses of serially collected isolates of Pseudomonas aeruginosa and Burkholderia species have provided in vivo evidence of adaptation during long-term infection of CF lungs (3,4,(6)(7)(8)(9)(10)(11)(12)(13). In contrast, we only scarcely understand how different members of the human microbiome adapt to new niches within the human body and become pathogens, and our knowledge is limited largely to Staphylococcus aureus evolutionary dynamics during CF lung infection and transitions from carriage to invasive disease (14)(15)(16)(17)(18).
Chronic obstructive pulmonary disease (COPD), the fourth leading cause of death globally, whose primary risk factor is smoking, is an irreversible airflow obstruction accompanied by emphysema, fibrosis, neutrophil airway infiltration, mucus hypersecretion, inflammation, and long-term lower airway colonization by opportunistic pathogens (19). Though nontypeable Haemophilus influenzae (NTHi) is typically a benign commensal of the human upper airways, it is also an opportunistic pathogen commonly isolated from the lower airways of COPD patients during both exacerbations and clinically stable periods, and it is responsible for 20 to 30% of all acute COPD exacerbations. NTHi persists within COPD patient lungs for months to years, contributing to chronic airway inflammation that results in worsening of symptoms and speeds disease progression (19).
In this study, we addressed molecular genetic changes underlying bacterial adaptation to the COPD airway by investigating the genetic variants arising from within-patient evolution of NTHi in an independent cohort of patients. We identify recurrent polymorphisms in several genes and provide experimental evidence for the biological significance of such convergent variation in the ompP1 (also called fadL) gene. Overall, our results provide novel insights into the evolution of a human-restricted respiratory commensal to be a persistent pathogen of the lower airway.

RESULTS
Longitudinal collection of NTHi genomes from COPD patients. To identify natural genetic variation in bacterial genomes associated with adaptation to the COPD lung, we serially collected NTHi isolates from a set of COPD patients over time and sequenced and comparatively analyzed their genomes. Sputum samples were collected from COPD patients during exacerbations requiring hospitalization or at primary health care routine visits; NTHi colonies were phenotypically identified and isolated. Thirteen patient series were selected on the basis of the presence of four or more longitudinally sampled NTHi isolates, as well as the presence of more than one independent isolate of the same pulsed-field gel electrophoresis (PFGE) type, suggesting persistent longterm infection by that strain. Sampled subjects were Spanish male smokers or exsmokers; clinical information on each subject is summarized in Table 1. The bacterial strain collection consists of 92 isolates, with 4 to 18 isolates per patient, collected over a period ranging from 1 to 9 years (Fig. 1a). Details on each isolate are summarized in Data Set S1a and b in the supplemental material.
The genomes of all isolates were subjected to Illumina-based shotgun wholegenome sequencing, and a subset of 17 isolates were also sequenced using the long-read Pacific Biosciences RSII technology (PacBio) to serve as complete reference sequences. Genome assembly, gene annotation, taxonomic classification, and in silico multilocus sequencing typing (MLST) were performed. Summary metrics are shown in Data Set S1c; results of DNA modification motif analysis for 10 PacBio assemblies are shown in Data Set S1d. Assemblies from the two sequencing methods showed high overall agreement ( Fig. S1; see also Materials and Methods).
Gene clustering and phylogenomic analysis. To provide context to these new COPD genomes for comparative genomic analyses, we combined our genome collection with 391 publicly available H. influenzae assemblies (of which 210 were from a recent large study of NTHi in COPD patients by Pettigrew et al. [36], 28 of which were of strains collected from COPD patients or patients with other lower airway infections, and the remaining 153 were non-COPD patient isolates), along with 35 representative genomes from sister species H. haemolyticus and H. parainfluenzae. Pan-genome analysis was conducted by clustering together homologous protein-coding genes using the pan-genome analysis pipeline Roary (3), resulting in 6,453 total homologous gene For molecular typing and clustering of closely related isolates into clonal types (CTs), besides performing PFGE and MLST analyses, we applied the goeBurst algorithm in PHYLOViZ to allele assignments that we made from 309 core protein-coding genes shared by all genomes, including the sister species. For each single-copy core gene, unique nucleotide sequences were defined as distinct alleles, and strings of allele identities (IDs) for each isolate were used as an MLST-like input. Strains with fewer than 15 allelic differences across all core genes were clustered together, and each resulting connected component was assigned an arbitrary CT number (Materials and Methods;  Table 1; see also Data Set S1a and b in the supplemental material). Axes indicate patient identifier (ID) numbers and the time of sampling for each isolate. Symbols indicate whether the isolate was collected outside an acute exacerbation (triangles) or not (circles). Colors indicate clonal type (CT) based on 309 single-copy core protein-coding genes. Black is reserved for CTs collected only once. Filled symbols indicate full-length fadL; open symbols indicate a truncated fadL allele. (b) Minimum-spanning tree (MST) of CTs across all genomes in the study. Nodes represent individual CTs (collapsed from the genotype-based nodes seen in Fig. S3), and edges indicate links between closely related CTs. Node size indicates the number of strains across the data set belonging to that CT, and node color indicates to which data set strains in each CT belong. Primary colors indicate that CTs contain a strain only from this work (red), from the work of Pettigrew et al. (36) (blue), or from other works (yellow). Mixtures of the primary colors indicate Ͼ1 strain from two data sets, with black indicating CTs containing strains from all three sets of strains. Fig. S3). So, CTs are defined as closely related groups of strains but allow some recombination. A clear-cut correlation between CTs and MLSTs was observed, since strains grouped by CT the same way as they grouped by MLST (Data Set S1b), but the former offers higher resolution for distinguishing relationships among CTs.
Overall, 205 CTs were detected. H. influenzae accounted for 175 of the CTs, including 40 of which appeared in our COPD collection. Twenty of the COPD CTs contained more than one isolate from within our collection, and so the 72 COPD isolates belonging to these 20 CTs formed the basis for the search for recurrent genomic changes (described below). Figure 1b shows a minimum-spanning tree (MST) of all CTs, indicating whether strains belonging to that CT were from this study, the Pettigrew et al. study (36), another study, or some combination. These results indicate that many CTs can be found in both Spanish and American COPD patients, as well as in altogether distinct subjects. Nearly half of the H. influenzae CTs across the data set had only a single isolate (n ϭ 87), suggesting that many CTs remain unobserved.
Consistently with results from clonal typing, the species-level phylogenetic tree (built from alignments of all protein-coding homologs with at most one copy per strain) illustrated several findings ( Fig. 2; unrooted version with outgroups in Fig. S4a; fully labeled rectangular version in Fig. S4b). (i) Individual patients were colonized by numerous genetically diverse strains over time, though diversity at any given time point remains unclear (24,32,(37)(38)(39)(40). High polyclonality and a high average degree of divergence among CTs contrasts with studies of environmentally acquired opportunistic pathogens, which have typically supported colonization by single CTs followed by subsequent diversification (3,6,41). (ii) Groups of closely related isolates were collected both within and among patient series, and they sometimes included isolates from public databases with diverse clinical and geographical origins. (iii) Several instances of probable interpatient transmission were observed (CT14, -44, -48, -73, -95, and -137).
Most notably, for CT48, 17 isolates were found across four patients, all of whom live within a 400-m radius, effectively neighbors. Furthermore, patients 4 and 6 were hospitalized on the same floor of the hospital in April 2011, patients 4 and 7 were hospitalized on the same floor in December 2012 and November 2013, and patients 6 and 7 were in the emergency unit at the same time in December 2012. No such relationship was found for patients 7 and 8 or 6 and 8. These observations, along with SNV differences among this subset of strains, suggest likely transmission links among these patients, as shown in the time course in Fig. 3 (Data Set S1a and b), though we cannot exclude the possibility that there may be indirect links through medical professionals, mutual friends, or relatives. Alternatively, since the relative abundance of CTs circulating in the population remains unknown, CT48 may simply be common in the local population.
Recurrent genetic changes among clonal types, especially in genes encoding membrane-associated functions. Polymorphisms among isolates of the same CT, especially when collected from the same patient, point to recent mutations or recombination events that may contribute to the bacterium's adapting to the COPD lung environment. Because a single isolate was collected at each sampling time point from what may be polyclonal populations, we did not treat earlier isolates as having parental genotypes and later ones as having derived genotypes (41,42). Nevertheless, independent polymorphisms at the same locus in different CTs and different patients might be indicative of adaptive evolution (or, alternatively, hotspots of mutation or recombination).
To rank genes by evidence of recurrent changes in distinct CTs, we first identified genetic variants (SNVs, short indels, and gene possession differences) that occurred within each of the 20 multi-isolate CTs. This showed that intraclonal genetic variation was characterized by relatively few genetic differences affecting a small number of loci, either as isolated point mutations or clusters of closely spaced genetic variants, characteristic of natural-transformation events (1,34,43,44) (Fig. 3; Data Set S3). We next classified variants based on their impact on protein coding and ranked homologous protein-coding gene clusters by the number of CTs affected by intra-CT polymor-phisms. Of 6,453 total gene clusters, 299 were polymorphic for short variants in one or more of the 20 multi-isolate CTs. Of these, only 15 gene clusters had intra-CT polymorphisms in three or more CTs (Table 2) and only 5 clusters had high-impact (frameshift, nonsense, or nonconservative missense) short variants in three or more CTs. The top hits in Table 2 affected by high-impact variants encoded membrane proteins, enzymes determining lipooligosaccharide (LOS) structure, and a DNA methyltransferase subunit. Several of these are known to be phase variable due to SSRs, whose number determines if the reading frame is open (lic2A, oafA, losA, hsdM3) (45)(46)(47). The frequent frameshift variants arising within CTs in these genes are due to copy number changes in these repeats due to slip-strand mispairing during replication (Data Set S3). This observation nicely corroborates the recent one by Pettigrew et al. with regard to an independent set of COPD isolates (36). Several hits in Table 2 were also found in the Pettigrew et al. COPD study (36), including DNA methylase, LOS-related glycosyltrans-  (36), 181 publicly available genomes, and 92 genomes from this study. The patient ID is indicated by node color, whereas CT is indicated by the color in the outer ring, as in Fig. 1a. Only COPD isolates sequenced for this study are labeled; Fig. S4a and b, respectively, show an unrooted phylogeny, including H. haemolyticus and H. parainfluenzae outgroups, and a fully labeled rectangular tree that includes bootstrap results.

FIG 3
Genetic variation within clonal type 48 (CT48). (a) Whole-genome view of all short variants detected for the 17 isolates in CT48, collected from 4 patients, against the PacBio reference assembly for isolate P621 (the x axis shows the P621 genome coordinate; the left-hand side of the y axis shows results for patients 4, 6, 7, and 8; blue hatch marks indicate SNVs and short indels). Blocks are due to numerous closely spaced variants, which resemble natural transformation events seen in the laboratory. (b) Zoomed-in view at 920 to 940 kb, showing that recombination tracts in this region likely had independent donor sources. (c) Heatmap of pairwise variant counts among CT48 isolates, showing 4 distinct haplotypes, mostly due to recombination tracts distinguishing the isolates. P618 in haplotype A is more similar to haplotypes B and C than they are to each other and, as seen above, appears to be a hybrid of the two. (d) Time course of all isolates belonging to CT48, encompassing strains collected from four patients. ferase, acetyltransferase, and a TonB-dependent receptor. Also among top hits were several genes affected by dozens of low-impact intra-CT variants (conservative missense and synonymous variants) affecting the same contiguous gene cluster (e.g., hgpB, hgpC, and tbp2), which were due to overlapping recombination tracts in different CTs ( Fig. 3) and potentially indicate loci under diversifying selection (48,49).
Independent loss-of-function mutations in ompP1 or fadL arise within multiple lineages and patients. Aside from changes in phase-variable genes, the most striking observation of the parallel evolution analysis above was that of independent predicted loss-of-function changes in the ompP1, or fadL, gene ( Fig. 1a; Fig. 4; Table 2; Data Set S2; Data Set S3). Seven of the 20 multi-isolate CTs were affected by polymorphisms in fadL, and this included nine independent null alleles induced by frameshift or nonsense mutations. Since only intra-CT polymorphisms were considered above and this was sometimes affected by annotation problems at truncated alleles, we extended the analysis by extracting the fadL allele from each of our new genome assemblies using tblastx (50). Multiple alignment of the resulting 92 amino acid sequences identified 34 distinct fadL alleles encoding a full-length protein, as well as 16 distinct alleles with mutations predicted to encode truncated proteins (Data Set S1b and e). Amino acid variation among the 34 full-length FadL protein variants is represented in Fig. 4a and shows high variability in predicted surface-exposed loops. Predicted full-length and truncated variants are illustrated in Fig. 4b and show clear evidence for multiple independent null mutations at different positions in the gene, including both frameshift and nonsense mutations. Truncation alleles were found in 24 isolates collected from 11 out of the 13 patients and spanning 14 CTs. Where possible, we identified  full-length alleles that were the likely ancestors of specific truncated variants, when present in sets of isolates from the same patient and CT with no fadL variation outside the putative null mutation (groupings in Fig. 4b and Data Set S1e).
We further provide direct evidence that the FadL protein is not expressed in strains carrying predicted null alleles by SDS-PAGE of bacterial whole-cell lysates from representative isolates. A protein band whose weight was consistent with the FadL predicted molecular weight was present in strains P589, P631, and P608 but absent from P590, P595, P662, P614, P622, and P607 (Fig. 4c). Band identity as FadL was confirmed for strain P589 by peptide mass fingerprinting matrix-assisted laser desorption-ionization tandem time of flight (MALDI-TOF/TOF) mass spectrometry. Collectively, these results , and two had more than one (CT44 had -100). (b) Multiple independent fadL truncation alleles were observed and are schematically shown to indicate surface-exposed loops (in gray); they are paired with full-length alleles, when possible. (c) Putative null alleles do not express FadL protein. Whole-protein extracts from cultures were analyzed by SDS-PAGE and Coomassie brilliant blue staining. A protein band at FadL's predicted molecular weight (arrows) was present in strains containing full-length fadL alleles (P589, P631, and P608; variants 2, 3, and 8, respectively) but absent in strains containing truncated alleles (P590, P595, P662, P614, P622, and P607; variants T1, T4, T6, T7, T8, and T10, respectively).
NTHi Genomic Adaptation to COPD Chronic Lung Infection ® strongly suggest that fadL loss-of-function mutations are selected for during NTHi's adaptation to the COPD lung.
FadL variation affects hCEACAM1-dependent NTHi interactions with host cells. FadL in NTHi (FadL NTHi ) strains has been reported to be a bacterial ligand for the human carcinoembryonic antigen-related cell adhesion molecule 1 (hCEACAM1), facilitating NTHi entry into epithelial cells (51,52). Therefore, we assessed the impact of naturally occurring variation in fadL genes on the interplay between NTHi and host cells. By comparing bacterial invasion of HeLa cells (which do not normally express hCEACAM1) to that of a stably transfected HeLa derivative cell line expressing hCEACAM1 (HeLa-BGP cells) (53), we first tested whether bacterial cell entry depends upon the interaction between FadL and hCEACAM1. With isogenic fadL knockout mutations in reference strains NTHi375 and RdKW20 (54,55), gentamicin protection assays showed strong defects in bacterial entry into HeLa cells in the absence of either FadL or hCEACAM1 in both strain backgrounds (Fig. 5a). As expected, COPD clinical isolates also showed FadL-dependent invasion of hCEACAM1-expressing cells. Pairs of clonal isolates carrying ancestral full-length and derived null alleles of FadL (Data Set S1e) were also tested. In all cases, invasion of hCEACAM1-expressing cells was substantially decreased in the absence of full-length FadL (Fig. 5b). Differences in hCEACAM1-dependent invasion rates were observed among isolates carrying distinct full-length FadL alleles, which could not be explained solely by differences in fadL gene expression levels; for example, isolate P608 showed the highest invasion rate, but its fadL gene expression was comparable to those of the other strains tested (Fig. 5c).
We next examined the interaction of FadL NTHi with hCEACAM1 by using computational molecular modeling, since a three-dimensional (3D) perspective of the molecular recognition events is lacking. To evaluate the potential effect of the natural variation present among the FadL surface-exposed loops, homology models of FadL variants from RdKW20, NTHi375, and P608 (selected as a representative COPD isolate with high rates of HeLa-BGP cell invasion) were computed using FadL in Escherichia coli (FadL E. coli ) and refined by molecular dynamics (MD) simulations. All three computed 3D structures were similar to each other and to that of FadL E. coli , though the mobile variable surface-exposed loops had distinct conformations among structures, and some differences were seen within the ␤-barrel (Fig. 6a, Fig. S5a to c, and Text S1). Docking calculations (56) between the three FadL NTHi models and hCEACAM1 were also performed and the resulting protein-protein complexes submitted to MD simulations, rendering high stability along the simulation. Protein-protein interaction analyses of the three FadL NTHi -hCEACAM1 complexes identified several relevant polar, hydrophobic, and CH-interactions that remained stable along the MD simulation, which varied depending on the allele used (Fig. 6b, Fig. S5d, and Text S1). Among the residues involved in the interactions, we found several hCEACAM1 residues (Ser32, Tyr34, Val39, Gln44, Gln89, and Ile91) that were previously reported as necessary for the interaction with NTHi by site-directed mutational analyses (57,58). Our structural modeling clearly identifies these residues as directly involved in the interaction with FadL, thereby accounting for the interaction between FadL NTHi and hCEACAM1 and further suggesting that natural variation in FadL surface-exposed loops may affect NTHi's ability to interact with epithelial cell surfaces. Our molecular modeling calculations also allowed us to identify several FadL residues located in loops 3, 4, and 5 as likely to be involved in hCEACAM1 binding. High conservation among FadL full-length allelic variants possibly excludes many of these from an association with variability in epithelial invasion. In particular, Gly193 P608 in L3, Ile286 RdKW20 and Lys290 RdKW20 in L4, and Asp332 P608 , Ala334 NTHi375 , and Lys339 RdKW20 in L5 are highly conserved. However, other identified residues are polymorphic, in particular, Val196 P608 , Asn200 P608 , Lys213 NTHi375 , Leu214 P608 , Lys216 P608 , Lys217 P608 , Trp220 P608 , Phe222 RdKW20 , Lys225 RdKW20 , Ser276 NTHi375 , and Gly279 P608 . These variants may account for the observed variation in terms of epithelial invasion by the tested strains.
FadL loss-of-function confers resistance to the bactericidal effect of arachidonic acid. FadL NTHi is also a predicted transporter of exogenous long-chain fatty acids HeLa-BGP (a stably transfected derivative expressing hCEACAM1) cells were used to quantify epithelial invasion by WT and ΔfadL derivative NTHi375 and RdKW20. WT strains invaded HeLa-BGP cells at significantly higher rates than they invaded HeLa cells (for NTHi375 and RdKW20, P was Ͻ0.0001), in contrast to their respective ΔfadL mutants. Inactivation of fadL decreased the invasion of HeLa-BGP cells compared to that of their respective isogenic WT strains (P Ͻ 0.0001). (b) HeLa and HeLa-BGP cells were also used to quantify invasion by COPD patient isolates with the putative full-length FadL ancestor and derived truncated allele pairs (strains belonging to CT9, -48, -40, -44, and -76; FadL full-length variants 2, 5, 8, 13, and 24; FadL truncated variants T1, T7, T8, T10, T11, T14, and T15). Strains with full-length alleles invaded HeLa-BGP cells at significantly higher rates than they invaded HeLa cells (for P589, P600, P608, and P666, P was Ͻ0.0005; for P634, P was Ͻ0.005), in contrast to strains with truncated alleles. Moreover, fadL truncation alleles significantly decreased bacterial entry into HeLa-BGP cells compared to that of their full-length ancestors (P was Ͻ0.0005 when P590 was compared to P589; the P value was the same for comparisons of P614 and P600, P622 and P600, P607 and P608, P617 and P634, P664 and P666, and P665 and P666). (c) Expression of the fadL gene during exponential growth in sBHI as measured by quantitative reverse transcription-PCR. Data are shown for strains containing full-length FadL variants assayed in panels a and b, namely, NTHi375, RdKW20, P589, P600, P608, P634, and P666. In all cases, means Ϯ standard deviations (SD) are shown. ΔCt, change in threshold cycle.
(LCFAs). In E. coli and Sinorhizobium meliloti, the fadL gene products are required for the ability to use free fatty acids as a sole carbon source (59)(60)(61), and key amino acid residues required for fatty acid binding and lateral diffusion have been elucidated for FadL E. coli (59,61). To investigate the ability of FadL NTHi to bind LCFAs, we performed and FadL P608 (in lime green). The FadL E. coli crystal structure is shown on the left (purple). The modeled structures correspond to a 14-stranded ␤-barrel with seven extracellular loops (L1 to L7), of which L1, L3, L4, and L7 are mobile and exposed to the outside. (b) Protein-protein docking between the amino-terminal domain of hCEACAM1 (green) and the modeled structures of FadL NTHi375 (blue), FadL RdKW20 (cyan), and FadL P608 (lime green). These interactions remained stable along molecular dynamics (MD) simulations of each hCEACAM1-FadL complex (Text S1).
docking calculations of three LCFAs in the FadL NTHi375 , FadL RdKW20 , and FadL P608 homology models. To take into account the flexibility of FadL, additional normal mode analysis (NMA) of the three FadL NTHi 3D models provided distinct conformations representing motions and conformational changes. Tested LCFAs included arachidonic acid (AA), oleic acid (OA), and lauryl dimethylamine-N-oxide (LDA) as a control (since it was previously crystallized with FadL E. coli ). All the calculations predicted several binding poses for the three fatty acids at the entrance to and within the FadL NTHi ␤-barrel ( Fig. 7 and Fig. S6). All three fatty acids were predicted to interact with FadL NTHi at sites equivalent to those for LDA in FadL E. coli , but other distinct binding sites were also identified inside the ␤-barrel (Text S1). The docking calculations provide reasonable binding poses to FadL NTHi for the studied LCFAs and support its role as a fatty acid transporter. However, further use of LCFAs by NTHi as a sole carbon source was deemed unlikely due to the absence of a complete ␤-oxidation pathway in NTHi (26,54). Thus, as expected, no bacterial growth was observed in a defined minimal medium free of fatty acids when supplemented with AA as a carbon source (62). , and FadL P608 (lime green) (c). The details of the interactions are depicted as sticks. Calculations predicted two favorable binding modes of AA in both structures (FadL NTHi375 and FadL RdKW20 ); the first one was at the entrance of the ␤-barrel (site A) in a hydrophobic groove located between L3 and L4. Interactions with polar residues, specifically with Lys at L5 (Lys336 in FadL NTHi375 and Lys426 in FadL RdKW20 ), lipophilic interactions, and CH-interactions with aliphatic and aromatic residues were identified; and a second binding mode at the depth of the ␤-barrel (site B), establishing polar contact with Lys142 present in L2. These binding poses were also found for FadL P608 together with an additional binding site (Text S1). Details of the interactions are represented as sticks.
NTHi Genomic Adaptation to COPD Chronic Lung Infection ® Free fatty acids act as natural detergents, and a bactericidal effect on H. influenzae has been described for AA (63). We next examined the viability of wild-type (WT) and isogenic fadL mutant strains when they were incubated in defined minimal medium with AA. WT strains had reduced viability in a dose-dependent manner; by contrast, both fadL knockout mutants were more resistant to AA than their respective WT strains (Fig. 8a). The same effect was seen for OA, though only at higher concentrations (data not shown). Resistance to AA was also assayed for two representative pairs of COPD isolates with predicted ancestral full-length and derived null alleles of FadL (Data Set S1e). In both cases, loss of FadL function was associated with increased bacterial resistance to AA (Fig. 8b). Although one or more of the few additional genetic variants that distinguish the two pairs of isolates may instead be responsible for this increased AA resistance, fadL loss of function was the only common one in both isolate pairs ( Table 3). Given that AA metabolites are key players in COPD-related airway inflammation (64), these molecules may in turn provide a selective pressure within COPD patient lungs such that NTHi fadL loss-of-function mutations are adaptive within this niche.
FadL inactivation has no effect in a murine model of acute lung infection. Human-restricted bacterial pathogens are known to selectively interact with hCEACAM1 (65). NTHi murine pulmonary infection may uncouple FadL functional analysis as an LCFA transporter and a hCEACAM1 ligand, since the latter is absent from mice. So, CD1 mice were infected with WT and fadL mutant NTHi375 and RdKW20 strains, and bacterial loads were quantified in lung and bronchoalveolar lavage fluid (BALF) samples. At 24 h postinfection (hpi), bacterial counts were comparable for both WT and mutant strain pairs. At 48 hpi, no differences were observed for WT and mutant RdKW20 strains; in contrast, NTHi375ΔfadL rendered lower bacterial counts than its isogenic WT strain (P Ͻ 0.005 in lung, P Ͻ 0.01 in BALF) (Fig. 8c and d). Together, these results showed that fadL gene inactivation gives no substantial benefit in this infection model. We acknowledge that this model system may be suboptimal in this case, since NTHi fails to establish a long-term infection, showing only variation in times to clearance, which occurs in Ͻ3 days for the WT strains.
Truncated fadL alleles appear to arise in NTHi within COPD patient lower airways and are rare in isolates from other body sites. Because all the subjects in this study had COPD for years prior to being sampled for NTHi, we could not conclude that null alleles of fadL arose within COPD patient lungs, despite 6 of 20 multi-isolate CTs being polymorphic for null alleles, since these might also arise in the normal commensal population of NTHi in the nasopharynx. To test whether loss of fadL function was specifically associated with lower airway infections, we first examined the distribution of full-length and truncated fadL alleles across 181 publicly available genomes (excluding those from the Pettigrew et al. study [36]) derived from isolates with diverse clinical origins (26,66,67) (Fig. S4). Of 153 isolates collected outside the context of lower respiratory tract infection (including nasopharyngeal carriage, middle ear carriage in pediatric patients with otitis media, invasive bacteremia, and meningitis isolates), only three non-lung isolates' genomes lacked a single full-length fadL annotation (1.9%, one from a patient with otitis media, one from a patient with bacteremia, and one other). In contrast, out of 28 isolates from the lower airways (isolates from COPD patient sputum or other pulmonary infections), 9 had truncated fadL genes (33.3%). This demonstrates that an intact fadL gene is the norm among NTHi strains as a whole, but loss-of-function mutations are significantly enriched in lung isolates (Fisher's exact test P value ϭ 1.3 ϫ 10 Ϫ6 ).
As we were completing this study, 269 new NTHi genomes collected from American COPD patients were published by Pettigrew et al. (36), providing us with a second independent test of the hypothesis that fadL null mutations are selected for during long-term COPD-associated infections. Indeed, 43 of all 269 genomes (16.0%) had truncated/absent fadL genes, as determined by Roary gene clustering analysis. Aggregating across assemblies gave consistent results: only 3 of 159 nonlung isolates lacked a predicted full-length fadL gene, whereas 69 of 399 lung/sputum isolates lacked it (1.9% versus 17.3%, Fisher's exact test P value ϭ 4.2 ϫ 10 Ϫ8 ).
A comparison of the frequencies of truncated fadL between the persistent and cleared isolates in the Pettigrew et al. data set found a weak (but not significant) trend toward higher rates of truncated fadL alleles in persistent isolates than in cleared isolates (37 of 202 persistent strains versus 6 of 67 cleared strains; odds ratio, 2.2; Fisher's exact test P value ϭ 0.084). Finally, because persistent isolates had been carefully chosen to pair the first and last isolates with matched MLSTs from the same patient, we also asked whether paired isolates tended to show loss or gain of fadL function over time. Of 101 MLST-matched pairs, 11 were polymorphic for fadL between the first and last isolation dates, and all 11 had a full-length allele at the first visit and a truncated allele at the last visit (binomial test P value ϭ 0.0001), providing strong evidence that fadL truncation mutations arise mostly in situ in the COPD lung.
Taken together, our results strongly suggest that intact fadL function is normally beneficial for H. influenzae but becomes deleterious in the context of lower respiratory tract infections, likely due to the fadL-dependent bactericidal effect of LCFAs being upregulated in COPD patient lungs.

DISCUSSION
Our comparative genomic analyses of longitudinally collected NTHi isolates from COPD patients provides novel insights into the evolutionary dynamics of a bacterium as it adapts to becoming a chronically infecting lung pathogen. This work and that of Pettigrew et al. (36) present data from two independent cohorts (Spanish and U.S. patients, respectively), and both identify an overlapping set of NTHi genes that change during long-term COPD infection. COPD strains show no phylogenetic clustering based on clinical source, geography, duration of persistence, or year of isolation, and phasevariable genes were a frequent source of genetic variation among clonally related isolates. Beyond the major findings in the work of Pettigrew et al. (36), we identified probable recombination and interpatient transmission events, as well as provided experimental evidence of the molecular mechanism behind selection for COPD-specific recurrent loss-of-function mutations in the NTHi gene ompP1, or fadL.
This study was motivated by an earlier one in which we conducted a limited genomic analysis of three serially collected NTHi isolates belonging to a single CT recovered from a single bronchiectasis patient, which showed increasing resistance to antimicrobial peptides and also carried amino acid substitutions in the SapABCDFZ system, which confers antimicrobial peptide resistance (68). Here, we expanded to a larger set of longitudinally collected COPD isolates. Notably, we identified intraclonal polymorphisms from 20 multi-isolate CTs collected over time, rather than in only the first and last in each patient series (36). We found evidence for parallel evolution in several NTHi genes, in which multiple CTs carried coding polymorphisms in the same genes, suggestive of in situ natural selection. In particular, the ompP1or fadL locus stood out because of the high frequency of independent loss-of-function alleles seen in different CTs and in different patients. Our observations are consistent with those of other studies that have identified amino acid variation in the FadL protein among H. influenzae and Neisseria meningitidis strains (51,69,70) and are congruent with reports on the key role of outer membrane proteins for within-host adaptive strategies (71,72).
The recurrent null mutations seen in fadL within NTHi CTs isolated from COPD patient sputum samples suggest a strong selective pressure to lose this gene function during adaptation to the COPD lung, consistent with NTHi's reductive evolution toward becoming a long-term resident of the lower airways (73). In support of this, our survey of Ͼ450 NTHi genomes found that null mutations in fadL were extremely rare except in isolates from lung infections. Given that the FadL NTHi -hCEACAM1 interaction was known (51) and strongly supported by our computational 3D models, we were initially puzzled how losing this interaction could be beneficial, but finding that loss of FadL also confers resistance to microbicidal fatty acids suggests antagonistic pleiotropy, wherein a single gene affects multiple phenotypes, each of which may be beneficial or deleterious depending on conditions. Antagonistic pleiotropy can be modulated by gene inactivation, such that loss-of-function mutants can outcompete the WT parental strains once a trait becomes deleterious (74). This may be the case for FadL: loss of the functional gene avoids triggering bacterial cell death by the combined fatty acid detergent effect and NTHi's inability to metabolize these types of molecules. The molecular basis for the fatty acid detergent effect on NTHi viability is unknown. Of note, Helicobacter pylori has features common to NTHi's susceptibility to free fatty acids, the lack of a ␤-oxidation pathway, and the presence of an active natural competence pathway (75)(76)(77)(78). Despite occupying separate niches of the human body, both pathogens might have common adaptive strategies contributing to colonization linked to chronic disease conditions. Other explanations for recurrent fadL loss-of-function mutations are plausible given that, as a surface protein, FadL may be subject to immune pressure and so bacteria may benefit from its loss.
Several lipids are molecular mediators of respiratory diseases (79). In particular, AA metabolites are key players in COPD airway inflammation (64). Phospholipase A 2 (PLA 2 )-catalyzed hydrolysis of membrane phospholipids results in production of free fatty acids, most importantly AA, which serves as a precursor for inflammatory mediators, such as platelet-activating factor (PAF) (80). In turn, lipoprotein-associated PLA 2 mediates PAF hydrolysis and dampens PAF-mediated inflammation (81). H. influenzae exploits molecular mimicry to evade host inflammation through GlpQ, an enzyme that hydrolyzes phosphorylcholine (ChoP) moieties shared by both PAF and the bacterial LOS molecule, therefore removing PAF from the airways and suppressing inflammation. Such bacterially driven regulation of PAF signaling might be overcome by a therapeutic increase of PAF levels at the site of infection through PLA 2 targeting (82), which in turn, might also dampen free AA release in COPD patient lungs, therefore limiting selection against FadL function during NTHi pathoadaptation (for a model, see Fig. 9).
Finally, we acknowledge that our strain collection presented some limitations. First, we are unlikely to have sampled the initial colonizing bacterial isolates; therefore, we did not order the observed genetic variations or determine the most recent common ancestor (13). Second, we analyzed single isolates from each sampling point, yet we know that NTHi colonization and infection are often polyclonal (24,32,(37)(38)(39)(40). Despite these considerations, our results show the utility of comparing the genomes of bacterial isolates collected over time from long-term chronic infections, showing that genetic polymorphisms within clonally related NTHi isolates can help tease apart transmission of strains between subjects and the genetic changes that result in changes in pathogenic traits. Here, only a small number of genes showed recurrent mutations in different bacterial CTs, and our follow-up studies support a bifunctional role for FadL in both interactions with host epithelia and interactions with free fatty acids.
We conclude that the environment of the COPD patient lung exerts a specific selective pressure to lose functional fadL, and this is likely by providing resistance to bactericidal fatty acids, in spite of gene truncation also compromising the ability of NTHi to invade epithelial surfaces. We speculate that FadL-hCEACAM1 binding may be essential for infection establishment but that, once established, persistence may favor the loss of FadL to overcome free fatty acids pressure. This hypothesis aligns with the observed low hCEACAM1 expression in lung tissue (83). Lastly, we highlight PLA 2 as a potential target for host-directed therapeutics against NTHi chronicity in the setting of COPD.
Bacterial growth conditions and genetic manipulations. H. influenzae COPD isolates (Data Set S1a) and reference strains RdKW20 and NTHi375 were grown at 37°C with 5% CO 2 on chocolate agar (bioMérieux) or supplemented brain heart infusion (sBHI), made from BHI supplemented with 10 g/ml hemin and 10 g/ml NAD, with or without agar for solid and liquid culture, respectively. Spectinomycin 30 g/ml (Spec 30 ) was used when needed. E. coli was grown on Luria Bertani (LB) broth, with or without agar, at 30°C or 37°C and supplemented with ampicillin at 100 g/ml (Amp 100 ) or Spec 50 , as needed. The fadL gene was knocked out by replacement of the coding sequence with a spectinomycin-resistant (Spec r ) cassette as described previously (86) (see Text S1 in the supplemental material and Table 4).
Genome sequencing, assembly, annotation, and data deposition. DNA was extracted from each isolate and sequenced using standard procedures (Text S1). All 92 genomic DNAs were subjected to Illumina paired-end short-read sequencing, and a subset of 17 DNAs were also subjected to Pacific Biosciences RSII long-read sequencing to generate covalently closed circular genomes. Genome assembly was platform specific (Text S1). Assemblies were taxonomically classified with Taxator-tk v1.2 (87). For MLST strain identification, the adk, atpG, frdB, fucK, mdh, pgi, and recA genes were extracted from assemblies and assigned STs (https://github.com/tseemann/mlst). New alleles and STs were submitted as needed (Data Set S1b; https://pubmlst.org/hinfluenzae). To compare PacBio and Illumina assemblies of . This study provides evidence for recurrent fadL mutations being a likely case of antagonistic pleiotropy associated with bacterial pathoadaptation in COPD patient lungs. In particular, loss of FadL function reduces NTHi's ability to infect epithelia but also increases its survival in the presence of bactericidal LCFAs, which are inflammatory markers enriched within the COPD lung. FadL-hCEACAM1 binding may be essential for virulence, but its persistence in the lower airways may indeed be favored by the inability of FadL to, thus, overcome free-LCFA selective pressure. OM, outer membrane; IM, inner membrane; CoASH, coenzyme A. the same isolates' genomes, progressiveMauve (snapshot_2015-02-25) (88) was used to reorder Illumina contigs against the complete PacBio assembly, and perfect agreement was found in all cases. To annotate strains with consistent gene names, a database of annotated Pasteurellaceae genomes from the NCBI database was created (https://github.com/rehrlich/prokka_database_maker). Assemblies were annotated using this database and Prokka (v1.11) (89). All raw reads and assemblies were deposited in NCBI under BioProject number PRJ282520; BioSample accession numbers for each strain are included in Data Set S1b and c. Publicly available genome assemblies at the NCBI and the Sanger Centre were downloaded (n ϭ 572) and reannotated with Prokka to maintain consistency. Assemblies were used only if they had Ͻ50 contigs, except for the Pettigrew et al. genomes (BioProject accession number PRJNA358390), where assemblies with Ͻ100 contigs were retained (except that all 269 assemblies were used in the fadL truncation analysis). This combined data set of 391 publicly available assemblies and 92 assemblies sequenced for this study were used for gene clustering, phylogenetic analysis, and clonal typing.
Gene clustering, phylogenetic analysis, and clonal typing. Homologous protein-coding genes among all genomes were clustered together using Roary at a BLASTp threshold of 75% (90). Phylogenetic tree reconstruction used the hybrid MPI/Pthreads version of RAxML (8.2.4, 10 processes, 16 central processing units [cpu] per process, -f a -x 1234 -m GTRGAMMA -p 7 -N autoMRE -j set thread affinity) on a concatenation of all codon-aware alignments of protein-coding genes with at most one copy per strain (91). PFGE and MLST are low-resolution strain typing methods, so we applied goeBurst via PHYLOViZ 2.0 (92) to cluster isolates into CTs, based on their allelic identity across 309 core protein-coding genes shared by all H. influenzae isolates studied here, as well as several H. haemolyticus and H. parainfluenzae strains. CTs were defined as sets of strains differing by no more than 15 allelic differences in the core gene set.
Detecting intraclonal genetic variation in protein-coding genes. Because methods to call genetic variants by short-read alignment to reference assemblies are relatively mature compared to methods that use whole-genome alignment, we chose a reference sequence from each CT (based on assembly quality, a Pacific Biosciences assembly when available) and then aligned short-read pairs from all isolates belonging to that CT to this reference (Burrows-Wheeler Aligner [bwa] mem version 0.7.5a-r405). Single-nucleotide and short-insertion/deletion variants (SNVs and indels, respectively) were extracted using freebayes (v1.1.0) (filtering out variants with a quality score of Ͻ1), and the impact of coding variants was evaluated using SnpEff (v4.3k) (93). To identify homologs across CTs, the LiftOver tool was used on the Roary gene possession table described above to determine gene identities between the CT reference genomes. Artifacts of gene annotation failure (for example, due to frameshift and nonsense mutations), particularly for the CT references, were reconciled by considering intraclonal variation in gene possession as a high-impact variant.
FadL allelic identities for whole collection. To extract the fadL gene sequence from each of the new NTHi genomes independently of the Prokka annotation, a local database containing the 92 whole-genome sequences was created and the local BLAST tool was employed (94). Best hits to fadL were extracted from each genome using bedtools (95), followed by translation and alignment using Clustal Omega.
Protein and gene expression analysis. Whole-protein extracts from bacterial cultures were analyzed on 12% SDS-PAGE 20-by 20-cm gels designed for long-range protein separation with maximum resolution under reducing conditions (Bio-Rad, PROTEAN II xi Cell), followed by Coomassie brilliant blue staining. Reverse transcriptase real-time quantitative PCR used total RNA extractions from exponentially grown bacterial cultures (optical density at 600 nm [OD 600 ] ϭ 0.6), with primers for fadL and primers for gyrA as endogenous controls (Table 4).
Free fatty acid susceptibility testing. Bacterial suspensions were incubated with various concentrations of free fatty acids or vehicle solution (ethanol), serially diluted, and plated on sBHI agar. Results are expressed as percentages of bacterial survival ([CFU ml Ϫ1 of LCFA/CFU ml Ϫ1 of the vehicle solution] ϫ 100).
NTHi mouse lung infection. A CD1 mouse model of NTHi pulmonary infection was used (97). At least five mice per treatment (genotype by time point) were intranasally inoculated with ϳ2 ϫ 10 8  (98), and subsequent MD simulations were performed with AMBER14 (98). Structure refinement of hCEACAM1 (PDB accession number 4QXW) was performed by means of MD simulations. Protein-protein docking of the three FadL NTHi structures with hCEACAM1 was performed with ZDOCK (56). LCFA docking calculations were performed with AutoDock 4.2.2 (99). Detailed protocols for the calculations can be found in Text S1.
Ethics approval and consent to participate. NTHi isolates were collected at the Hospital Universitari Bellvitge (HUB), a tertiary-care hospital in Barcelona, Spain. NTHi strains were isolated from good-quality sputum samples of 13 COPD patients over a 14-year period (2000 to 2014) at regular check-ups or during visits requiring hospitalization due to an exacerbation of the disease. Informed consent was not required, as sputum sampling and the microbial isolation process are part of the standard microbiological routine. Patient confidentiality was protected in all cases. For animal experiments, CD1 female mice aged 4 to 5 weeks were purchased from Charles River Laboratories (France) and housed under pathogen-free conditions at the Institute of Agrobiotechnology facilities (registration number ES/31-2016-000002-CR-SU-US). Animal handling and procedures were in accordance with the current European (Directive 86/609/EEC) and National (Real Decreto 53/2013) legislations, according to the FELASA and ARRIVE guidelines, and with the approval of the Universidad Pública de Navarra (UPNa) Animal Experimentation Committee (Comité de Ética, Experimentación Animal y Bioseguridad) and the local government authorization (approval with reference number PI-022/15).
Availability of data. All raw reads and assemblies were deposited in the NCBI database under BioProject number PRJ282520; BioSample accession numbers for each strain are included in Data Set S1B and C.
TEXT S1, DOCX file, 0.1 MB. We declare that we have no competing interests.