We used P. aeruginosa whole genome sequencing of isolates longitudinally sampled from 33 young Danish CF patients, as described previously (Marvig et al., 2013, 2015b). Patients were either infected with one clone type that was sampled repeatedly over many years, harbored co-occurring clone types, or were first colonized by one clone type that was replaced by another with no detected co-occurrence (Marvig et al., 2015a,b). The collection contained sequence data from 441 isolates of 50 clone types. The genomes of P. aeruginosa from the young patients were de novo assembled, as described previously (Marvig et al., 2015a) and are available at the European Nucleotide Archive, accession number PRJEB5438 (www.ebi.ac.uk/ena). The genomes were not assembled to one chromosome but were comprised of multiple contigs (for clinical isolates on average 136 ± 274 SD, Table S1), which may result in not all TA systems being identified. For each contig open reading frames (ORFs) were predicted with the program Prodigal (Prodigal.ornl.gov; Hyatt et al., 2010). In addition, we included two transmissible clone types, DK1 (41 isolates) and DK2 (44 isolates), which were transmitted between patients from 1973 up to the 1990's as previously described (Jelsbak et al., 2007; Yang et al., 2011; Marvig et al., 2013; Markussen et al., 2014). These were primarily from older patients with one or a few isolates from each. In younger patients, following measures to separate patients in the clinic, the majority of Danish CF infections are caused by environmental strains. To test whether TA systems of the P. aeruginosa that establish CF infection differ from those that do not, we included 12 environmental isolates collected from non-infectious sources described in Hilker et al. (2015; sequences available at ENA accession no. PRJEB4961 and de novo assemblies kindly provided by J. Klockgether) and Ghoul et al. (2015; sequences and de novo assemblies available at http://pubmlst.org/paeruginosa/; accession no. WTCHG-44355_209 to 213). The clone types of environmental isolates were determined following Marvig et al. (2015b); three were found to overlap with clone types found in clinical isolates while nine were distinct new types.

De novo genome assemblies of the first isolate of each clone type from each patient were analyzed to find the number and family of putative TA loci, for a total of 69 first isolates of 52 clone types from 35 patients, and 12 environmental isolates of 12 different clone types. Putative type II toxin-antitoxin systems were identified with the program TAfinder (http://202.120.12.133/TAfinder/index.php; Shao et al., 2011). The cut-off was set at a score of an average shared identity of the toxin and the antitoxin at 80% with previously computationally and experimentally identified TA systems. Systems classified by TAfinder were subjected to BLASTP analysis to determine conserved domains, and a BLASTP analysis excluding Pseudomonas hits to see if they were found, as a pair, in other species. TAfinder failed to annotate as open reading frames one toxin (HigB) and one antitoxin (ParD) that had previously been identified. The antitoxin HigA (PA4674) and the toxin ParE (PA0729) were collected from the ORF files as described below. The coding sequences of the unannotated toxin and antitoxin were identified by a BLASTN search against a local database of all the clinical isolates, and these were translated to protein sequences with EMBOSS Transeq (http://www.ebi.ac.uk/Tools/st/emboss_transeq/) for the first isolate from each clone type. A third previously identified pair, HicBA, was identified by BLAST search of the toxin sequence against the de novo assemblies. The diversity within a given toxin and antitoxin was calculated in Clustal Omega (http://www.ebi.ac.uk/Tools/msa/clustalo/), and unique protein sequences were identified with GenomeTools (http://genometools.org/) and compared to known sequences with BLASTP. The unique sequences were aligned with the reference sequence from TAfinder with Clustal Omega. Protein sequences without a 100% identity match were submitted to GenBank. We refer to all of these putative TA systems as TA systems, but note that for the majority their function still needs to be confirmed experimentally.

The TA systems not found in all isolates were considered to be part of the accessory genome, on GIs and potentially mobile elements. The location of TA systems in the P. aeruginosa accessory genome was determined by BLASTP analysis restricted to the PAO1 reference genome. We attempted to identify where GIs were inserted into the core-genome: If a TA pair was found on a 100 kb contig, the predicted ORFs of the entire contig were subjected to BLASTP analysis against the reference PAO1 genome. A potential GI insertion site was determined when genes with high shared identity to genes in PAO1 (>90%) were found next to genes with lower shared identity (>40%). The insertion site was confirmed by identifying an upstream gene with high similarity to a gene in PAO1, neighboring the suspected insertion site. The inferred insertion sites of GIs were compared to those identified by Mathee et al. (2008).

To identify mutations in TA loci acquired during infection we compared isolates of the same clone type sampled longitudinally from patients. We aligned de novo assembled contigs of the region of interest using BLAST align or mapped raw reads to representative sequences of the clone type. The reads were mapped to reference sequences with the Burrows-Wheeler alignment tool (http://bio-bwa.sourceforge.net, version 07.4) and alignments generated with the paired-end reads or single-end reads setting, respectively. The alignments were filtered to remove unmapped reads, sorted and indexed, and assigned to a read group using Picard Tools (http://picard.sourceforge.net, version 1.91). Differences between isolates of the same clone type were identified using Samtools (http://samtools.sourceforge.net, version 01.19) and mutations manually checked using IGV (http://www.broadinstitute.org/igv, version 2.3.61).

Deletions and insertions of >50 bp in GIs were identified in the sequence alignments described above. To identify gain of TA systems during infection we used the first 5–8 amino acids of a toxin or antitoxin gene, provided by TAfinder, as a search motive to collect the genes from all isolates from the young patients from the files with predicted ORFs. The use of these short motives allowed for discovery of TAs already identified in other isolates that potentially had been vertically transmitted between clone types within patients, and the identity of positive hits was verified by BLASTP search.

Nine clone types found in multiple patients has been shown not to be transmitted among patients, but independently acquired by different patients from the environment (Marvig et al., 2015b; Table S1). To estimate the timescale of gain and loss of GIs with TAs in the environment we used the variation in presence of GIs of these isolates. The number of SNPs between isolates (Marvig et al., 2015b) was used as a measure of time since the last common ancestor, given an estimated average rate of 2.5 SNPs/year (Yang et al., 2010). The three environmental isolates of the same clone types as clinical isolates were not included in this analysis because the sequence coverage was lower, and the SNPs identified therefore not reliable (but see Table S2 for comparisons).

All clone types from the young patients were classified based on their success in establishing an infection; where “low” was used for clone types sampled from a patient once, or multiple times over less than a year (n = 29), “medium” for clone types sampled multiple times infecting for longer than a year but either co-occurring with other clone types or being replaced (n = 16) and “high” for clone types consistently sampled with a frequency above 80% of sampled isolates from a patient infecting for more than a year (n = 24). These were compared with the 12 environmental samples. We used 1-Way ANOVA to test if the number of TA systems were significantly different between categories and a χ²-test to analyse if the presence of GIs varied significantly between categories.

All isolates harbored TA systems, and we found extensive diversity across the studied isolates. TAfinder identified 32 distinct TA systems, by screening all first sampled isolates of clone types from all patients (69 isolates of 52 clone types from 35 patients) and 12 environmental isolates. Of these, 8 antitoxins, 18 toxins, and 7 pairs were verified by BLASTP analysis, as they produced a specific hit to conserved domains of toxin or antitoxins. Sixteen T and As had an unspecific hit to conserved domains of toxin or antitoxins [E < 0.01, Table S3 and Supporting Online Material (Data Sheet 1)]. For the putative TAs without a clear match we did a BLASTP analyses with an Entrez Query specified as toxin or antitoxin, and all except two shared >27% identity with putative TAs from various species [Table S3 and Supporting Online Material (Data Sheet 1)].

We classify the TA systems into four groups based on whether the T and A in a pair are both verified, only one partner is verified or both unverified: in category 1 both the toxin and antitoxin were verified by a specific hit by BLASTP analysis (seven pairs), or experimentally confirmed (two pairs); in category 2 either the toxin or antitoxin has a specific hit but the partner does not (10 pairs); in category 3 one partner has an unspecific hit and the partner does not (seven pairs); and in category 4 neither were verified (six pairs). These latter six pairs, however, all had conserved domains that Makarova et al. (2009) identified as strong candidates for TA systems. Makarova et al. (2009) suggest that neighboring genes to known T or As present strong candidates for novel components, and in the subsequent analyses, we include all pairs from TAfinder where one or both components were verified by unrestricted BLASTP analysis. We leave the category 4 pairs out of the subsequent analyses, although we note that these may be functional TAs given their sequence similarity to putative systems from other species. For that reason, they are listed in Tables S1, S3 and Figure 1, and when mentioned in the text referred to as unverified systems. As further support for the identified pairs being TA systems, all were also found, as pairs, in a diverse range of other species (Table S3). Each isolate harbored between 4 and 11 systems of categories 1–3 (mean 5.89 systems ± 1.64 SD, Table S1). Four TA systems were part of the core genome: PA1029/PA1030 (Xre/COG5654, category 3), PA1878/PA1879 (Xre/PIN, category 3), PA0124/PA0125 (RHH/RelE, category 2), and PA4674/ unannotated (HigAB, category 1; Tables 1, 2; found in all environmental and first sampled clinical isolates, except one that was missing the HigAB pair). For the Xre/COG5654 and the Xre/PIN systems the toxins were not classified as such by BLASTP analysis, but the PIN toxin shared 68% identity with previously identified toxins in the restricted search. The RHH/RelE and HigAB pairs had been identified previously (Pandey and Gerdes, 2005), and the latter confirmed experimentally (Wood and Wood, 2016). All of these core genome pairs were also found in the P. aeruginosa reference genomes PA1, PA7, UCBPP-PA14, and LESB58. The other 22 TA systems varied in prevalence from 1% (found in one clone type in one isolate) to 65% (1–2 copies in 48 first sampled isolates of 34 CF clone types and 4 environmental isolates). Our analyses suggest that these are located on GIs with specific insertion sites (Table S3). Therefore, variation in TA systems across isolates was caused by those that reside in the accessory genome, on potentially mobile elements.

The TA loci were categorized by conserved domains, and for some as belonging to a (super-) family. Based on TAfinders classification there were six different toxin families/domains, and eight different antitoxin families/domains represented (and additionally three unverified toxins and four antitoxins). We named each TA pair according to the family or domain (Makarova et al., 2009) as identified by TAfinder, but note that these are not always identical to the conserved domains identified by BLASTP (Table S3), and protein length in amino acids. For each toxin and antitoxin family, 1–15 proteins were identified (Figure 1; Table 1). All but 12 protein sequences had a 100% identity match in GenBank (Table S3). The 12 unique sequences were given the accession numbers KY855415-KY855426. Alignments of the unique sequences to the reference sequences from TAfinder are available in the Supplementary Online Material (Data Sheet 1). Three TA systems previously identified in PAO1 were not annotated by TAfinder. When the protein sequences of these were analyzed in TAfinder, they shared 67% (HigAB, classified as Xre/RelE by TAfinder, category 1), 72% (ParDE, classified as PHD/RelE, category 1) and 53% (HicBA, category 1) identity with pairs recognized by TAfinder. This suggests that the chosen cut-off at 80% identity is somewhat conservative and leaves some TA systems unidentified, however, lowering this obviously increases the risk of including false positives.

Comparison of the genes flanking TA pairs showed that they were often highly conserved among clone types, and the GIs with a given TA pair were in most cases found in one or two locations in the genome (Tables S1, S3). Many of these insertion sites have already been identified as sites with “pathogenicity islands” (Liang et al., 2001; He et al., 2004; Qiu et al., 2006; Battle et al., 2008) and “regions of genomic plasticity” (Mathee et al., 2008; Table S3). BLASTP analyses of the GIs showed that they frequently harbored “mobility genes” such as transposases and integrases, suggesting that they are, or have been, mobile elements. The frequently observed PHD/RelE pair (category 1) is known to reside on the prophage Pf4 (Webb et al., 2004), which is important for biofilm development (Rice et al., 2009).

We commonly observed two TA pairs on the same GI, as marked in pink in Table S1. For example, the Xre/RelE 96/99 pair (category 2) was only found together with RHH/RelE 90/115-140 (category 1), while the latter was alone on a GI in 75% of the clone types. The AbrB/RelE system (category 1) was sometimes found in the place of the unverified system Xre/COG5654 208/210 (category 4), with the same flanking genes. This could suggest that the latter may also be a functional TA pair.

Not all TA pairs were found on contigs long enough to determine their insertion site in the genome and provide the full-length sequence of the GI. This was often the case when two GIs of the same type were present, suggesting that the presence of highly similar regions hinders de novo genome assembly. The insertion sites of an “orphan” GI (TA systems found on a contig with no link to the core genome) could often be inferred when specific insertion sites for that GI had been identified for other clone types: If the known insertion sites were found on separate contigs flanked by GI-specific genes this was assumed to be the insertion site. In a few cases no insertion site was identified. It cannot be ruled out that these TA pairs were found on plasmids, however the sequencing coverage of these regions were comparable to the average of the isolates, suggesting that if so they would be low-copy plasmids.

In the longitudinal samples, we only observed two mutational events in core genome TA loci during infection; the complete loss of a pair as part of a large >40 kb deletion in DK13, and a non-synonymous SNP in a toxin in one isolate of DK1 (Table 2). All samples of DK21 from one patient had ca. 40 kb insertion disrupting the Xre/COG5654 pair (category 4). Further, three clone types, that were sampled only once from a patient and from the environment, harbored stop codons or an insertion in toxin or antitoxin genes in the core genome (DK23, DK49, and MG7), and one had an apparent deletion and insertion covering the HigAB pair (DK28, category 1; Table 2).

The mobility of some GIs was supported by the finding of a circular GI in one patient: the GI comprised an individual contig, where the end-genes were found in the middle, and the core genome flanking genes were found next to each other on a separate contig (RHH/RelE 90/115-140, DK11, category 1). Further, two GIs with unverified TA systems were found to double their copy number during infection (Table 2). However, we found no evidence for horizontal transmission of GIs with TA systems between clone types within patients. When clone types co-occurring in a patient shared a GI type and TA system, the sequences had >50 SNPs, which would not be expected if they had been transmitted recently, given the average rate of mutation of 2.5 SNPs/year (Marvig et al., 2013).

Transmission between patients is evident for three clone types, DK6, DK15 and DK26 (Marvig et al., 2015b; Table S1), where the GIs with TA systems were identical for DK6 and DK15. For DK26, isolates from one patient (P99F4) showed an apparent acquisition of a >212 kb GI not found in other isolates, with a TA pair (COG5606/RelE 111/125, category 3) found in two other clone types. In another patient with DK26 (P92F3), there was a duplication of an already present GI with an unverified TA pair, as mentioned above (Table 2; Table S1).

Complete or partial loss of GIs carrying TA systems during infection was observed 15 times in the sample population, in 7 out of 52 clone types. One isolate of IT11 lost a GI, that is present in two copies in the other isolates from the patient. In four isolates of the transmissible DK2 clone type, we observed the complete excision of the most common GI (with RHH/RelE 90/115-140, category 1) in four independent events, after 29 (three isolates) and 33 (one isolate) years of infection, respectively. This GI was also lost in one isolate of DK44. Partial loss of a GI including two unverified TA systems occurred once in DK26. However, there were eight events in four clonetypes, where TA systems remained despite partial loss of a GI, in an isolate of DK3, an isolate of DK24, two isolates of DK36 and four isolates of DK2 (unverified category 4 pairs; Table 2).

To estimate the rate of loss and gain of GIs with TA systems on an ecological time scale, we used the nine clone types that had been independently acquired by multiple patients from the environment. For these we had information on the number of SNPs in the core genome differentiating them (Marvig et al., 2015b). As an example, for clone type DK6, the isolates most closely related to each other, sampled from two patients, are separated by 421 SNPs (~168 years). Three GIs are only present in one of the DK6 lines giving an estimated rate of 0.033 losses or gains of GIs per year, or one event per 30.4 years. This gives us 11 independent estimates of gain and loss of GIs, as two clone types contributed two data points each (Table S2). We recorded between 0 and 6 events that had occurred over an estimated 21.6–201.6 years. This lead to on average one loss or gain every 52.69 ± 42.21 years std. dev. and a range between an event every 9.9–151.6 years (excluding the data point with 0 events).

There was no significant difference in the number of TA systems of a clone type and its ability to establish a chronic infection, as defined in the Materials and Methods (One-Way ANOVA, df = 3; 33.67; F = 0.65, p = 0.59). The most commonly occurring RHH/RelE 90/115-140 TA system (category 1) residing on GIs was, however, found more often in clinical isolates compared to environmental isolates (70% of first clinical isolates compared to 33% of environmental, χ² = 4.4, df = 1, p < 0.05, with the caveat of a small sample size of environmental isolates). Prevalence of the RHH/RelE 90/115-140 TA system across different categories of clinical isolates was similar, ranging from 64 to 75%.

We observed few mutations in the core genome TA loci during infection. Notably, the few isolates with stop-codons and an insertion in core genome antitoxins were not successful at establishing infection, as they were either found as clone types only sampled once or from an environmental isolate. Two of these, interestingly, had the same mutation introducing a stop codon (Table S3). This corresponds with other studies (Mine et al., 2009) where TA loci are not acquiring loss of function mutations and, therefore, are not selected against but maintained as selfish elements or for other physiological cellular processes.

Complete or partial loss of GIs during infection was rare, suggesting that they are beneficial pathogenicity islands facilitating host colonization (He et al., 2004) and/or that TA systems are successful at stabilizing them. In six out of 15 cases an entire GI was lost, and because these were found in the same locations as documented ICEs this may represent active excision of mobile elements. Excision may be a way for the ICE to “jump ship” in search of a new host in the event of cell damage. This was seen in Pseudomonas knackmussi where the fraction of the cell population in an ICE “transfer-competent” state increased in the presence of reactive oxygen species and cell damage (Reinhard and van der Meer, 2014). In eight of the nine partial losses of GIs, the TA systems were intact, which could indicate that TAs can stabilize a part if not all of a GI. Even though pathogenicity islands are well-documented to be beneficial during infection (He et al., 2004; Battle et al., 2008), horizontal transfer of these or other GIs was not observed between clone types within patients. This may reflect that infecting isolates already have the GIs they need or that spatial segregation of co-occurring clone types within the lung and sinus environment limits the possibility for transfer.

To evaluate the stability of GIs in patients, we used the variation in TA-bearing GIs within clone types that had been independently acquired from the environment to estimate the rate of loss and gain of these on an ecological timescale. This gave an estimate of one loss or gain of a GI every 53 years. There are some obvious caveats with this estimate: the number of SNPs in the core genome only provides a rough estimate of the time since the last common ancestor, and GIs may have been lost and acquired repeatedly in the time between sampling of the compared isolates. That said, the estimate of one event per 53 years is in line with what we observe in the clinical isolates where loss occurs rarely: in the longest sampled clone type, DK2, we recorded losses only four times after 29 and 33 years of infection (Table 2; three of these deletions occurred in isolates with an increased mutation rate, so-called hypermutators with a mutation rate >3.6 times higher than the other isolates (Marvig et al., 2013). Comparable stability of GIs in the environment and in patients suggests that they are not only maintained for pathogenesis, and further supports a stabilizing effect of TAs.

Our findings show that within P. aeruginosa there is no correlation between an incremental number of TA systems and the ability to establish a chronic infection. Previous work in S. typhimurium and E. coli have shown an association between certain TA systems as well as a correlation between the cumulative number of TA systems and level of persistence (Maisonneuve et al., 2011; Helaine and Kugelberg, 2014; Helaine et al., 2014) and there is a strong correlation between the two across eight Pseudomonas species (Vogwill et al., 2016). They have also been implicated in inducing biofilm formation in E. coli and enhancing host virulence in H. influenza, Mycobacterium tuberculosis, and S. typhimuriun (Lobato-Marquez et al., 2015), and inducing persistence in the development of recalcitrant infections in P. aeruginosa (Mulcahy et al., 2010).

We find that GIs with the RHH/RelE 90/115-140 TA pair were found significantly more frequently in clinical than environmental isolates, supporting a correlation between infectivity and these TA-bearing GIs. They were found at two insertion hotspots where multiple pathogenicity islands have been described (in the intergenic regions PA0976//PA0988 and PA4541//PA4542; Table S3), from which some genes have been shown to affect virulence in infection models (He et al., 2004; Battle et al., 2008). The GIs sharing these insertion sites are likely mosaic in nature, with recombination occurring within and between isolates, and even between species (Larbig et al., 2002; He et al., 2004). A similar pattern of some specific TAs being found more often in pathogenic isolates has been observed in the related Pseudomonas putida (Molina et al., 2016).

Our results suggest that some TA systems in P. aeruginosa may serve in the maintenance of GIs and do not appear to play a role in persistence or may not even be functional. P. aeruginosa is highly versatile and ubiquitous in the environment, as well as an opportunistic pathogen (Stover et al., 2000). Therefore, it may rely on a baseline number of TA systems, to mediate stress responses while colonizing/infecting a diverse range of environments and host types with variable abiotic conditions. If so, harboring additional TA systems may not provide greater benefits to strains, specifically in infection compared to other environments, and the core genome TA systems could be all that is needed.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.