|Home | About | Journals | Submit | Contact Us | Français|
Renibacterium salmoninarum is the causative agent of bacterial kidney disease, an important disease of farmed and wild salmonid fish worldwide. Despite the wide spatiotemporal distribution of this disease and habitat pressures ranging from the natural environment to aquaculture and rivers to marine environments, little variation has been observed in the R. salmoninarum genome. Here we use the coverage depth from genomic sequencing corroborated by real-time quantitative PCR to detect copy number variation (CNV) among the genes of R. salmoninarum. CNV was primarily limited to the known dominant virulence factors msa and p22. Among 68 isolates representing the UK, Norway and North America, the msa gene ranged from two to five identical copies and the p22 gene ranged from one to five copies. CNV for these two genes co-occurred, suggesting they may be functionally linked. Isolates carrying CNV were phylogenetically restricted and originated predominantly from sites in North America, rather than the UK or Norway. Although both phylogenetic relationship and geographical origin were found to correlate with CNV status, geographical origin was a much stronger predictor than phylogeny, suggesting a role for local selection pressures in the repeated emergence and maintenance of this trait.
1. The sequence data for all isolates used in this study is available for download from http://trace.ncbi.nlm.nih.gov/Traces/sra/?study=ERP003780
This article identifies expansive duplication of the genes encoding the dominant virulence factors msa and p22 in the fish pathogen Renibacterium salmoninarum, the organism responsible for bacterial kidney disease of salmonid fish. R. salmoninarum is a highly clonal bacterium with a very limited accessory genome, and although duplication of msa is already known as a concept, this study extends the finding to p22, the other major surface protein. The number of identical gene copies may in some cases be as high as five. The data suggest multiple independent duplication events that appear to be much more common in strains circulating in the Pacific Northwest region of North America, pointing to local selection pressures as important for the repeated emergence.
Gene copy number variation in bacteria is probably severely underreported, and there are very few reports on the regional distribution of the phenomenon. It is hoped that the findings and methodology presented in this article may serve to fuel the interest in performing gene copy number variation studies, as the mechanism is increasingly being seen as more frequent and phenotypically important than previously believed.
Renibacterium salmoninarum is the causative agent of bacterial kidney disease (BKD) in cultured and wild salmonid fish. BKD can result in acute morbidity or mortality, or it can be a slowly progressive disease causing an often dramatic decline in growth. BKD is economically important in aquaculture, where it can spread horizontally throughout sea pens of juvenile and subadult Atlantic salmon (Salmo salar) (Murray et al., 2012) or vertically through transferred broodstock or eggs (Evelyn et al., 1986). It is also a concern for conservation and restoration efforts for endangered fish stocks because infections are prevalent among more susceptible free-ranging Pacific salmon in river and marine systems (Pascho et al., 1993; Rhodes et al., 2011; Sandell et al., 2015).
Although the pathogenicity of R. salmoninarum is incompletely understood, several antigenic determinants have been described, including the dominant immunogenic protein major soluble antigen (MSA) (Turaga et al., 1987; Wiens and Kaattari, 1991), an abundant heat-stable 57 kDa extracellular protein that makes up 60–70% of all surface proteins in R. salmoninarum (Fredriksen et al., 1997; Wood and Kaattari, 1996), and is involved in immunosuppression (Brown et al., 1996; Fredriksen et al., 1997; Turaga et al., 1987), agglutination (Senson & Stevenson, 1999; Wiens et al., 1999; Wiens and Kaattari, 1991) and virulence (Coady et al., 2006; O’Farrell & Strom, 1999; Senson & Stevenson, 1999). Other antigenic determinants include capsular synthesis, heme acquisition operons, haemolysins and an immunosuppressive 22kDa surface protein provisionally named p22 (Fredriksen et al., 1997). The p22 gene encodes a poorly described loosely associated surface protein (Fredriksen & Bakken, 1994) that has been implicated in suppression of antibody production and a stronger agglutination of leucocytes than that which is seen for the MSA protein (Fredriksen et al., 1997).
The genome of the type strain of R. salmoninarum, ATCC 33209T, contains two identical transcriptionally active copies of the MSA-encoding gene: msa1 and msa2 (O’Farrell & Strom, 1999; Rhodes et al., 2002). Both genes are essential for the development of clinical disease and mortality (Coady et al., 2006). Whilst it seems certain that a single copy was originally acquired through horizontal gene transfer and subsequently duplicated within the bacterial genome (Wiens et al., 2008), the origin of this gene is unclear, as no homologue to the msa gene has ever been found in any other sequenced genome. Both msa loci are flanked by insertion sequences and transposases, and msa2 is additionally flanked by several degraded genes related to conjugation (including traA relaxase, type IV secretion protein and site-specific recombinase resolvase). Because multiple copies of identical genes are unusual in bacterial genomes, O’Farrell & Strom (1999) suggested that multiple msa copies might confer a selective advantage. Subsequently, Rhodes et al. (2004) demonstrated the presence of a third copy in some isolates, and provided clear evidence for a positive correlation between msa copy number and mortality at lower infection doses.
The gene content variation of this species appears to be exceptionally low, with core- and pan-genomes reported to be very similar even for strains sampled over 50 years from a wide range of habitats (Brynildsrud et al., 2014). However, this does not include paralogues, and the findings of Rhodes et al. (2004) suggest that copy number variation (CNV) in the msa genes of R. salmoninarum has phenotypic relevance. Gene duplication has been shown to be adaptive in bacteria (Riehle et al., 2001), and CNV is known to be an important mechanism for dose variation of specific proteins under appropriate environmental conditions (Stranger et al., 2007). As an example, a recent study demonstrated that some strains of Mycobacterium tuberculosis harboured a large, tandem gene duplication and noted greater expression of an anaerobic survival regulon that is contained within the duplication (Domenech et al., 2010).
The aim of the present study was to screen a diverse collection of R. salmoninarum isolates for evidence of CNV in any of the core genes and, if found, to investigate phylogenetic and spatial patterns of the distribution of genetic variants. This work can provide a better understanding of Renibacterium microevolution and may shed light on the mechanisms of differential disease manifestation in different populations.
Sixty-eight isolates whose spatial and temporal origins varied widely were sequenced on an Illumina GAII platform at The Genome Analysis Centre (TGAC), Norwich, UK, as part of a previous effort by the authors, and are available at the Sequence Read Archive of the National Center for Biotechnology Information (NCBI) under the accession numbers listed in Table 1. Non-pairing reads, reads containing ambiguous characters and reads with an average PHRED score of <20 were discarded before alignment to reference genome ATCC 33209 (available from NCBI GenBank under accession number NC010168) with Geneious v7.1 (Biomatters), using the option to randomly map reads with multiple best hits.
CNVs were discovered using the R (R Development Core Team, 2012) package CNOGpro (Brynildsrud et al., 2015) with the following parameters: coverage counted in sliding windows of length 50 bp, prior probability of changing states (for each read count observation) was set to p=1.0×10–10 and the error-rate parameter was set to 0.01. The runHMM method was used to call CNV regions and copy numbers were considered correct if they agreed with credible intervals (percentiles 1–99) from the runBootstrap method. When evaluating results we discarded IS994 tallies, as 69 copies (69 orfA and 67 orfB) of this element are known to exist in the reference genome (Wiens et al., 2008), making it impossible to evaluate copy number variation with our method. We also considered standalone CNV calls in segments shorter than 300 bp as unreliable, as such calls could happen from chance alone (Brynildsrud et al., 2015). When quantifying total msa enrichment, the signal from msa1 and msa2 were added together, and the relative frequencies were inferred by inspecting the signal from the hypothetical protein-encoding gene p12 (Fig. S1, available in the online Supplementary Material). These results were corroborated using real-time quantitative PCR (qPCR) on selected isolates with different copy number multiplicities (detailed in the Supplementary Material).
The presence or absence of msa CNV in an isolate was considered a binary trait, and associations between this trait and year of isolation, host species and saltwater/freshwater habitat were investigated by logistic regression, both bivariable and multivariable with interaction terms using R.
Phylogenetic trees were created from single nucleotide polymorphism alignments with the program MrBayes (Ronquist & Huelsenbeck, 2003) (see Supplementary Material). Pairwise patristic distances between isolates were calculated as the sum of branch lengths between leaf pairs of the consensus tree. Pairwise geodesic distances between isolates' geographical origins were calculated by solving for central angle in the spherical law of cosines and multiplying by the radius of the Earth. The latitude–longitude coordinates were rounded to the nearest degree. In some cases the exact sample origin was not known, so the coordinate pair was set to represent geographical midpoints for the sub-national region. To test for phylogenetic and spatial clustering of CNV presence/absence, we created a binary matrix where equal CNV statuses of isolate pairs were coded as 1 and unequal as 0. In this analysis we regarded isolates with asterisks listed in Table 1 as positive for the duplication in question and the remaining isolates as negative. We then adopted a Mantel test-like approach by performing the Mann–Whitney U test of equal distributions between groups defined by CNV status on patristic/geodesic distance data. This test estimator was subsequently compared with those obtained from 10000 random permutations of the CNV status matrix. The trait was considered to be phylogenetically or spatially clustered if the test estimator fell below the lower 1-percentile limit in the distribution of permuted data set estimators.
Overall, very little CNV was seen in our isolates. In fact, the coverage data of most isolates (57/68) indicated no variation at all. This finding is consistent with previous reports of a high degree of sequence conservation in the R. salmoninarum genome. Nevertheless, CNV was found in 11 isolates, shown in Fig. 1.
A complete list of all CNVs discovered in this study can be found in Table 2. In total, there were nine distinct CNV regions. Four of these were unique to the Carson5b isolate and two to isolate 5006. The remaining CNVs were non-unique and occurred jointly (i.e. the presence of one CNV type also implied the presence of the others) in all 11 CNV isolates. Among these were duplications of the genes encoding the primary surface surface proteins of R. salmoninarum: the msa gene and a 22kDa hypothetical protein (hereafter referred to as p22).
The total number of msa copies in CNV-positive isolates ranged from two to five. This confirms the supposition that the minimum copy number of msa genes is two, as no isolate presented a read coverage that was suggestive of only a single copy. There were two different msa duplication types, for which we provisionally introduce the nomenclature ‘type I’ and ‘type II’. Type II was a subunit of type I, but the two can be differentiated by type II’s lack of a marker gene, p12 (a predicted gene annotated Rsal33209_1032). As this gene is only part of type I duplications, the relative frequency of the two types can be found by inspecting the coverage of the p12 gene (Fig. S1).
Type I msa duplication included the msa gene, the p12 marker gene, the transposase-encoding Rsal33209_0133 and the inactivated insertion sequence (IS) sequence ISRs3, including all intergenic segments and flanking inverted IS994 sequences. Type I msa duplication thus very closely resembles the genomic region roughly between coordinates 110000 and 115000 in ATCC 33209, and is surely a duplication of the msa1 gene.
Type II msa duplication included the msa gene with the intergenic sequence from the terminus of the gene and roughly 800 bp downstream, which resembles two different regions of ATCC 33209: coordinates 110400 to 112901 or 945077 to 947575 in ATCC 33209. We could therefore not determine whether type II duplications represent duplications of msa1 or msa2, and unfortunately read mapping proved unhelpful to investigate this. Although the msa1 and msa2 genes differ very slightly at upstream and downstream sites, the ORFs themselves are identical, and there are several large (130–180 bp) inverted and direct repeats plus one 91 bp perfect palindrome associated with the gene, confounding read mapping (Fig. 2). However, previous experiments have only found duplications of msa1 (Rhodes et al., 2004). An msa1 origin must also be suspected for our data due to the fact that the traA.2 gene neighbouring msa2 was not duplicated in any isolates.
It remains unknown whether msa loci are differentially regulated. Using the terminator prediction tool ARNold (Naville et al., 2011) and the RibEx riboswitch explorer (Abreu-Goodger & Merino, 2005), we discovered that the palindrome at the 3' of the msa ORF contained a predicted rho-independent terminator/riboswitch-like element at both the msa loci, although with ‘G’ as the central loop nucleotide for msa1 and ‘C’ for msa2, opening the possibility for riboswitch-mediated regulation (Fig. 2).
The third non-unique CNV region matched the region between coordinates 2965759 and 2967751 in ATCC 33209. This region is flanked by inverted IS elements and contains a single ORF, encoding the p22 protein (a 22kDa hypothetical protein labelled RSal33209_3334). Also part of the duplication unit was the intergenic segments on both sides of this ORF. The total number of p22 copies in msa-duplicated isolates was estimated as ranging from two to five.
The msa–p22 duplication trait did not correlate with year of isolation, host species or saltwater/freshwater habitat. However, a strong geographical pattern was seen in the presence of gene duplication. CNV was absent in the exclusively European lineage 2 (lineage notation from Brynildsrud et al., 2014), and limited to defined clusters within the widely distributed lineage 1A and the Pacific Northwest-associated lineage 1B. Among the 10 isolates containing additional copies of msa and p22 genes, six are from the Pacific Northwest, three are from Eastern North America (New Brunswick, Canada) and oneis from Norway, corresponding to 55, 43 and 8% of the total investigated isolates from each respective region. Notably, of the 36 UK isolates, not one displayed CNV (Fig. 3).
To test whether CNV was clustered within different phylogenetically and spatially defined groups, we used Mantel correlation analyses (Fig. 1). For geodesic data, we found the Mann–Whitney U estimator to be 255344, compared with the full range 432002–515531 from the permuted dataset, which translates to a p-value of <1.0×10–4 when calculated conservatively as in Diniz-Filho et al. (2013). However, because the distribution of U values follows a near-perfect normal distribution (as calculated by the Anderson–Darling test of normality), a parametric p-value estimation of p<1.0×10–50 can be used (Fig. 4). In other words, CNV was strongly clustered into geographically defined groups. This can also happen because phylogenetically related isolates tend be spatially clustered, so we also investigated whether the pairwise patristic distances between isolates impacted the CNV. For these data, Mann–Whitney's U was computed as 419090, which is also lower than the full range of all permuted-matrix values (424 489–525 215) (non-parametric p<1.0×10–4; Gaussian parameterization p=7.4×10–5). Although this implies association between patristic distance and CNV as well, the pairwise geodesic distance is a much stronger predictor of CNV status, implying that these gene duplications are primarily the result of local selection pressures.
Although duplicates of genes encoding ribosomal and transfer RNA subunits are well known in bacteria, other gene duplication–amplification events are only now gaining attention and have probably been underreported in the literature (Andersson & Hughes, 2009; Elliott et al., 2013). CNV in msa has been documented previously (Rhodes et al., 2004), and the current study extends CNV to p12 and p22.
No isolates were found to have CNV in one gene but not the others, suggesting that these genes have a functional interaction relationship in which increased copies of either are not valuable without concomitant copy number increases of the other, or that the genes are somehow duplicated together due to linkage. The latter possibility is perhaps somewhat marginalized by the known genomic distance between msa and p22, which in ATCC 33209 is around 300 000 bp between msa1 and p22, going through the origin. However, it is possible that these genes are more closely located in strains other than ATCC 33209.
Although we have detected several large gene duplications, we have not been able to predict their relative orientation and distance to each other or to the rest of the chromosome. Wiens and Dale suggest a plasmid context of msa3, based on variable hybridization intensity in Southern blots, and another possible scenario could be the association of msa3 with a phage, as an unconfirmed observation of an R. salmoninarum phage was previously reported (Fryer & Lannan, 1993). Both msa1 and msa2 are flanked by inverted IS sequences, notably IS994, and IS3-like insertion sequences as well as other ORFs with high homology to transposable elements and transposases, suggesting that they could be transferred and integrated through recombination or transposition mechanisms, although duplications have only been documented in msa1 (Rhodes et al., 2004)
Ten of the 68 isolates screened in this study (~15%) displayed an increased copy number msa, p12 and p22 genotype, in stark contrast to the 19 of 26 isolates (~73%) that Rhodes et al. (2004) found to be msa3-positive. In their paper, every isolate except two (MT239 and GL64, which are msa3-negative) were from the Pacific Northwest region of the USA, suggesting that the strains circulating in that particular region have a higher frequency of multiple-copy msa genotypes. A predominantly North American CNV distribution is also consistent with the findings of Wiens & Dale (2008), who observed the msa3 gene in North American but not in European isolates. In this study, the geographical origin was a much stronger predictor of CNV status than the inferred patristic distance from the phylogenetic tree. Isolates 05372K and Cow-Chs-94 for example are of lineage 1B origin, and thus thought to have diverged from lineage 1A isolates such as Carson5b between 100 and 700 years ago (Brynildsrud et al., 2014). In spite of this these isolates all have duplications of the msa, p12 and p22 genes, a trait not shared by phylogenetic neighbours of these isolates. Note, however, that they are all sampled from fish originating from the Columbia River main basin, where multiple fish stocks co-occur. The fact that we observe this pattern of low intra-cluster but high inter-cluster patristic distances and that isolates originate from multiple geographical locations across North America (and, in a single case, Norway), sampled over a 19 year period from five different species of salmon from both freshwater and saltwater habitats, strongly suggests multiple independent introductions of the trait rather than simple inheritance.
Importantly, gene copy numbers varied widely across isolates displaying CNV. This has a number of apparent implications. Firstly, it seems that two is the basic copy number of the msa gene, as this genotype was by far the most common across lineages and ecosystems, was the genotype of the oldest isolates, and no isolates contained fewer than two msa copies. The diverse duplication pattern thus points to a base number of two msa genes with subsequent copy number expansions as a more parsimonious explanation than higher-value msa copy number and subsequent gene loss. Secondly, this heterogeneous duplication pattern indicates locally restricted gene duplication–amplification events rather than prevailing ecotypes as an explanation for the geographical clustering of CNV.
It is not clear to what extent the duplications we have found in the present work impact overall pathogen fitness. One possibility is that the observed duplications in fact represent selfish mobile genetic elements. However, this possibility contradicts the current understanding of the msa gene, as two copies have been proposed to confer selective advantage (O’Farrel & Strom, 1999). The immediate benefit of duplications could be through modulation of protein dosage under variable environmental conditions, while the long-term advantage is that the extra copies can, over time, accumulate mutations and evolve new functions (Conant & Wolfe, 2008; Kondrashov, 2012; Kondrashov et al., 2002). In favour of a selectionist explanation is the observation that these duplications are seemingly not immediately removed from the population, but rather shared by related isolates and thus perhaps maintained in local populations. (Isolates 05372K and Cow-chs-94, for example, are closely related despite being from separate river systems and isolated 11 years apart, and they both have multiple duplications of the msa, p12 and p22 genes, although the exact numbers of each gene appear to vary.)
Rhodes et al. (2004) found that the presence of a third msa copy was clearly associated with increased mortality at lower, environmentally relevant doses. It is therefore tempting to speculate that the additional copies that we have found are increasingly beneficial to the bacterium. Such duplication–amplification events of immunomodulatory genes are now thought to be common under adaptation to new, extreme and variable environments (Elliott et al., 2013), and these results point to a higher extent of such selection pressures in the Pacific Northwest than elsewhere. Our findings suggest that extra msa copies interact with the relatively unknown p22 protein, as the two were always duplicated together. The nature of this interaction remains unknown and more research is needed to conclusively determine the relative fitness- and virulence relationships between different duplication-value R. salmoninarum isolates.
O.B.B. and S.F.N. are funded through NMBU. S.G. is funded through the Norwegian Research Council and VaxxiNova Norway. E.J.F. is funded by BBSRC/NERC grant WGS-AQUA (BB/M026388/1). L.D.R. is supported by NOAA Fisheries.