Cross-Immunity and Community Structure of a Multiple-Strain Pathogen in the Tick Vector

Ixodes ricinus ticks were sampled monthly from the Bois de l'Hôpital site in a deciduous forest near Neuchâtel, Switzerland (47°00′55.6″N, 6°94′16.7″E; surface area of 1 ha), over a period of 11 years (2000 to 2010). We collected questing ticks by dragging a white flag (surface area of 1 m²) over the vegetation along four transects of 100 m. After collection, the species and stage of each tick were identified. In the laboratory, each tick was cut into two halves using sterile scissors. One tick half was used to determine spirochete infection via immunofluorescence microscopy (56). The other tick half was used to grow Borrelia isolates as follows. Half-ticks were placed in culture tubes containing BSK II medium (57), incubated at 34°C, and examined by dark-field microscopy every 10 days for 2 months. DNA was extracted from all spirochete-positive BSK II medium cultures. Briefly, 1 ml of culture medium was centrifuged and washed two times with phosphate-buffered saline before the pellet was suspended in 50 ml of ultrafiltered water. Pellets were incubated at 100°C for 10 min, and the thermolysates were stored at −20°C. Borrelia species were identified using a PCR-reverse line blot (RLB) assay following the method of Schouls et al. (58). The PCR protocol amplifies the 23S-5S spacer gene of B. burgdorferi sensu lato and was previously described by Alekseev et al. (59), and the RLB assay probes were described by Burri et al. (60).

As of 2012, the isolate archive at the University of Neuchâtel contained ∼3,000 tick-derived spirochete isolates for which the Borrelia species had been identified by RLB assay. In this region of Switzerland, I. ricinus ticks are commonly infected with the following Borrelia species: B. afzelii, B. bavariensis, B. burgdorferi sensu stricto, B. garinii, and B. valaisiana (14). In this study, we focused on B. afzelii and B. garinii because they are the two most common species in the isolate archive at the University of Neuchâtel and because they are the most common etiological agents of LB in Europe. We included only Borrelia isolates derived from questing nymphs because this stage takes only one blood meal during the previous larval stage. We included only those isolates that, according to the RLB assay, were infected with a single Borrelia species (either B. afzelii or B. garinii) so that we could assign major ospC groups with confidence to the correct Borrelia species. We established a database of all the nymph-derived spirochete isolates of either B. afzelii or B. garinii obtained from 2000 to 2010. For each of the two Borrelia species, we randomly sampled a maximum of 20 nymph-derived isolates for each of the 11 years. Using this sampling strategy, we obtained 193 isolates of B. afzelii and 190 isolates of B. garinii (Table 1).

The ospC gene was amplified using a nested PCR method consisting of two consecutive amplification steps that was described by Bunikis et al. (49). We modified the oligonucleotides for the second PCR to include 454 Life Science's A or B sequencing adapter, the key, and the template-specific primers (see the supplemental material). The forward primer was tagged for forward sequencing. The reaction was performed in an automatic thermocycler (Eppendorf, USA). We checked the success of the PCRs by visualizing the amplicons on a 1% agarose gel stained with Midori green (Labgene, Switzerland). The PCR products were purified using a DNA purification kit (Wizard SV Gel and PCR cleanup system; Promega, Switzerland). We created eight different pools of Borrelia isolates for the 454 sequencing run (40 to 50 isolates per pool). Isolates were randomly assigned to pools subject to the constraint that each pool contained a balanced number of isolates from each year. Sequencing was done on a 454 Roche GS FLX sequencing apparatus and was outsourced to Microsynth AG (Balgach, Switzerland).

The 454 sequencing run produced 352,882 sequences of the ospC gene. After cleaning the data set (see the supplemental material), we retained 240,410 ospC gene sequences that were each 521 bp long.

We defined an a priori set of seed sequences prior to conducting the cluster analysis. Seed sequences were defined as those sequences that occurred at least 50 times in our cleaned data set. The seed sequences were selected using the cd-hit-dup function of CD-HIT software (61, 62) with an error rate of 1%, which produced a set of 331 seed sequences. These sequences were used as seeds to perform a cluster analysis of the data from the cleaned data set with similarity thresholds ranging from 98% to 91% using the software CD-HIT. In addition, we created a phylogenetic tree (1,000 bootstraps) of the 331 seed sequences using the PhyML program (63). This tree confirmed the reliability of the clustering of the seeds. To have a measure of genetic similarity of the major ospC groups, we calculated the mean pairwise genetic distances using the distmat program (64) with a Kimura model. Each mean pairwise genetic distance was based on 30 sequences of each major ospC group in the pair. We also calculated the focal pairwise genetic distance for each major ospC group. This coefficient measures the average genetic distance between a focal major ospC group and all other major ospC groups.

The nature of the genetic variation observed within each major ospC group was investigated. Only the first 350 nucleotides of the ospC gene sequences were analyzed because the end of the forward sequences contained too much sequencing error. For each major ospC group, 30 sequences were selected from the seed sequences that had been used for the clustering (a total of 690 sequences). For major ospC groups with less than 30 seed sequences, nonseed sequences (sequences that occurred less than 50 times in our data set) were added. Sequences were aligned manually, and the numbers of synonymous and nonsynonymous substitutions, insertions, and deletions were recorded.

Several papers have dealt with the nomenclature of the major ospC groups in B. burgdorferi sensu lato (46, 47, 49, 54, 65, 66). In the United States, workers have developed a well-established system of names for the major ospC groups of B. burgdorferi sensu stricto (45). In Europe, in contrast, there is no universally adopted naming system for the major ospC groups of B. afzelii and B. garinii (47, 49). In the current study, we present a summary of all the different names that have been used depending on the author (see Tables S1 and S2 in the supplemental material). In addition, we suggest that future work on the major ospC groups of these two Borrelia species follow the nomenclature developed by Bunikis et al. (49), Hellgren et al. (50), and Strandh and Råberg (54) for B. afzelii and that developed by Lagal et al. (47) for B. garinii. Here, we use “ospC strain” to refer to a genetically distinct entity characterized by a particular major ospC group.

To confirm the reliability of the RLB assay, we performed additional sequencing of a subset of 120 isolates. According to the RLB assay, 60 of these isolates were singly infected with B. afzelii and 60 were singly infected with B. garinii. For each of these 120 samples, we used PCR to amplify the recA gene and the hbb gene according to the protocol described by Richter et al. (67). The amplicons were purified using a Qiagen purification kit. Purified amplicons were sent to Microsynth AG (Balgach, Switzerland) for Sanger sequencing.

All statistical analyses were performed using R (68).

The absolute prevalence of a pathogen strain is its frequency in the entire population (uninfected and infected individuals). The relative prevalence of a pathogen strain is its frequency in the subpopulation of infected individuals. In the present study, the absolute prevalence of infection was determined for each Borrelia species by the RLB assay. We then used 454 sequencing to estimate the relative prevalence of each major ospC group in a randomly selected subset of Borrelia isolates. Some analyses of multiple-strain pathogen communities require estimates of the absolute prevalence. We therefore converted our relative prevalences to absolute prevalences by including the fraction of uninfected ticks that corresponded to the fraction of infected ticks that had been processed by 454 sequencing for each of the 22 combinations of Borrelia species and year (Table 1). For example, if we subsampled half of the infected ticks for a particular combination of Borrelia species and year, then we included half of the uninfected ticks for that combination of species and year. In this way, our data set was a random and representative sample of the populations of infected and uninfected ticks (Table 1). All analyses described below were conducted separately for each Borrelia species.

Strain richness refers to the number of different ospC strains that were found within a single tick (including uninfected ticks). To test whether strains were aggregated inside ticks, we compared the observed distribution of strain richness to the Poisson distribution, where the variance is equal to the mean (69). Aggregation occurs when the number of ticks with multiple strains is greater than the Poisson expectation (variance/mean > 1). We used a bootstrap analysis to calculate the 95% confidence limits of the variance-to-mean ratio.

The strength and the direction of each pairwise association were estimated using the association coefficient described by Pielou (70). We used the absolute frequencies to calculate this pairwise association coefficient, which has a range from −1 (a complete negative association) to +1 (a complete positive association) and is undefined if all hosts are uninfected or infected with both strains. For each pair of strains, the statistical significance of the pairwise association was assessed using a two-tailed Fisher exact test. The significance level was adjusted for multiple comparisons using the Bonferroni correction (for B. afzelii, P = 0.05/171 pairs = 0.000292; for B. garinii, P = 0.05/153 pairs = 0.000327).

We also calculated the focal pairwise association coefficient for each major ospC group. This coefficient is a measure of whether a focal major ospC group is more or less likely to co-occur with other ospC strains. Strains with high focal pairwise association coefficients (+1) are more likely to be found in multiple infections than strains with low focal pairwise association coefficients (−1).

Finally, we tested whether strains with genetically dissimilar major ospC groups were more likely to occur together inside the same tick. Specifically, we tested for a correlation between the pairwise associations and the pairwise genetic distances using the Mantel test in the vegan package in R (71).

The ospC gene sequence data have been deposited in the Sequence Read Archive under BioProject PRJNA293785 with the accession number SRP063760. A type sequence of each major ospC group of B. afzelii and B. garinii is available in GenBank (see the supplemental material for the GenBank accession numbers). The GenBank accession numbers for major ospC groups V1 of B. valaisiana and Q of B. burgdorferi sensu stricto are AF093483 and JQ253799, respectively.

For each of the two Borrelia species, we randomly sampled a maximum of 20 nymph-derived isolates for each of the 11 years. Using this sampling strategy, we obtained 193 isolates of B. afzelii and 190 isolates of B. garinii (Table 1). The 454 sequencing run produced 352,882 sequences of the ospC gene. After cleaning the data set (see the supplemental material), we retained 240,410 ospC gene sequences that were each 521 bp long. There were 23 different major ospC groups in our local population of I. ricinus ticks, and this number was stable across a range of similarity thresholds (93% to 98%; see Fig. S1 in the supplemental material). According to BLAST analysis of the 23 major ospC groups, 10 belonged to B. afzelii (ospC groups A1, A2, A3, A5, A7, A9, A10, A11, A12, and A14), 11 belonged to B. garinii (ospC groups G2, G4, G6, G7, G8, G9, G10, G11, G13, G14, and G15), 1 belonged to B. valaisiana (ospC group V1), and 1 belonged to B. burgdorferi sensu stricto (ospC group Q). All of the 23 major ospC groups had been previously described in the literature (1, 46, 47, 72).

There were a number of major ospC groups that were found in both Borrelia species. Five of the 10 B. afzelii major ospC groups (A1, A9, A10, A11, A14) were found in isolates that were singly infected with B. garinii according to the RLB assay. Conversely, 7 of the 11 B. garinii major ospC groups (G2, G7, G8, G9, G11, G13, G14) were found in isolates that were singly infected with B. afzelii according to the RLB assay. In what follows, we use the terms “native” and “exotic” to distinguish between these two types of major ospC groups. Thus, B. afzelii had 10 native and 9 exotic ospC groups (7 B. garinii-derived ospC groups and the V1 and Q ospC groups) for a total of 19 major ospC groups (Table 2), and B. garinii had 11 native and 7 exotic ospC groups (5 B. afzelii-derived ospC groups and the V1 and Q ospC groups) for a total of 18 major ospC groups (Table 3).

The major ospC groups were more evenly distributed in B. garinii than B. afzelii. In B. afzelii, the strain of ospC group A10 (referred to here as strain A10) was found in 54.40% of the infected nymphs (Table 2). Three strains (A1, A9, A14) were found in 24.35% to 31.09% of the infected nymphs, and the remaining 15 strains were found in less than 11% of the infected nymphs (Table 2). In B. garinii, strain G8 was found in 42.11% of the infected nymphs, and six other strains were found in 21.05% to 31.58% of the infected nymphs (Table 3). Six strains were found in 10.53% to 14.74% of the infected nymphs, whereas the five remaining strains were found in less than 6% of the infected nymphs. In B. garinii, the exotic major ospC group A10 was very common and occurred in 31.58% of the infected nymphs.

For B. afzelii, the mean pairwise genetic distance within each of the 10 native major ospC groups was 1.06% (n = 44,250 pairwise comparisons, range = 0.42% to 2.16%). For B. garinii, the mean pairwise genetic distance within each of the 11 native major ospC groups was 1.14% (n = 54,285 pairwise comparisons, range = 0.66% to 1.63%).

Across the 690 sequences belonging to the 23 major ospC groups, there were 140 nucleotide substitutions, of which 58.9% were synonymous (83/141) and 41.1% were nonsynonymous (58/141). Single nucleotide insertions and deletions and codon deletions were common, whereas double nucleotide insertions or deletions were rare (see Table S3 in the supplemental material).

We calculated the pairwise genetic distances between the 21 native major ospC groups (i.e., the exotic major ospC groups were excluded). The conspecific or allospecific pairwise genetic distance refers to whether the two native major ospC groups belonged to the same Borrelia species or to different Borrelia species, respectively. The mean conspecific pairwise genetic distances of B. afzelii (n = 45 pairwise comparisons, mean = 17.84%, range = 8.86% to 23.59%; Fig. 1) and B. garinii (n = 55 pairwise comparisons, mean = 18.54%, range = 12.77% to 26.36%; Fig. 1) were not significantly different (t = −1.13, degrees of freedom [df] = 98, P = 0.261). The mean allospecific pairwise genetic distance (n = 110 comparisons, mean = 20.94%, range = 10.71% to 25.93%; Fig. 1) was significantly larger than the mean conspecific pairwise genetic distance (t = 6.921, df = 208, P < 0.001).

The strain richness of the infected nymphal ticks was calculated using both the native and exotic major ospC groups for each Borrelia species. Most nymphal ticks were infected with multiple ospC strains. For B. afzelii (19 major ospC groups), 78.76% (152/193) of the infected nymphal ticks carried multiple ospC strains. For B. garinii (18 major ospC groups), 84.74% (161/190) of the infected nymphal ticks carried multiple ospC strains. The mean strain richness of B. garinii (2.96 strains per infected tick, range = 1 to 11 strains per tick) was 25% higher than that of B. afzelii (2.37 strains per infected tick, range = 1 to 7 strains per tick), and this difference was statistically significant (t = −3.764, df = 381, P < 0.001).

The variance-to-mean ratio of strain richness was calculated using the absolute frequencies of both the native and exotic major ospC groups for each Borrelia species. The variance-to-mean ratio of strain richness was significantly greater than 1.0 for both Borrelia species, indicating that the ospC strains were highly aggregated in the nymphal ticks (Fig. 2). The variance-to-mean ratio of strain richness was 47.4% higher in B. garinii (mean = 3.98, 95% confidence interval [CI] = 3.44 to 4.51) than in B. afzelii (mean = 2.70, 95% CI = 2.50 to 2.91).

The pairwise association coefficients within each Borrelia species were calculated using the absolute frequencies of both the native and exotic major ospC groups. For B. afzelii (19 major ospC groups), 66.7% (114/171) of the pairwise association coefficients were positive (mean = 0.08, range = −0.01 to 0.35). All of the 39 statistically significant (P < 0.05/171) pairwise associations were positive (mean = 0.20, range = 0.11 to 0.35). For B. garinii (18 major ospC groups), 77.8% (119/153) of the pairwise association coefficients were positive (mean = 0.15, range = −0.004 to 0.48). Again, all of the 78 statistically significant (P < 0.05/153) pairwise associations were positive (mean = 0.25, range = 0.11 to 0.48). The mean pairwise association coefficient of B. garinii was significantly greater than that of B. afzelii (t = 6.10, df = 322, P < 0.001). The results were the same for the subset of 117 statistically significant pairwise association coefficients (t = 3.54, df = 115, P < 0.001).

The mean focal pairwise association index of B. garinii (n = 18 major ospC groups, mean = 0.35, range = 0.12 to 0.63) was 20.7% higher than that of B. afzelii (n = 19 major ospC groups, mean = 0.29, range = −0.005 to 0.62), but this difference was not significant (t = 1.35, df = 35, P = 0.187). Some ospC strains, such as B. afzelii ospC strain A10 and B. garinii ospC strains G8, G14, and A10, had very high focal pairwise association coefficients (>0.50).

The relative frequency (F2 in Tables 2 and 3) of each major ospC group (transformed by the square root) was a highly significant predictor of the focal pairwise association index in B. afzelii (n = 19 major ospC groups, F_{1, 17} = 565.7, P < 0.001, r² = 0.971; Fig. 3) and in B. garinii (n = 18 major ospC groups, F_{1, 16} = 1,688.0, P < 0.001, r² = 0.991; Fig. 3). An analysis of covariance found no significant interaction between Borrelia species and relative frequency on the pairwise association index (F_{1, 33} = 0.9, P = 0.354). The slope of the relationship between the focal pairwise association index and the square root-transformed relative frequency was therefore the same for each Borrelia species (0.90 ± 0.044). After controlling for the relative frequency, the focal pairwise association index remained significantly higher for B. garinii than B. afzelii (F_{1, 34} = 4.5, P = 0.041).

The correlation between the pairwise genetic distance between the major ospC groups and the pairwise association index of the major ospC groups inside the ticks was not statistically significant for either B. afzelii (n = 171 pairwise elements; Mantel test, r = −0.041, P = 0.630; Fig. 4) or B. garinii (n = 153 pairwise elements; Mantel test, r = −0.053, P = 0.664; Fig. 4). Similarly, the correlation between the focal pairwise genetic distance and the focal pairwise association index was not statistically significant for either B. afzelii (n = 19 major ospC groups; Pearson correlation test, r = 0.032, t = 0.13, df = 17, P = 0.897; Fig. 5) or B. garinii (n = 18 major ospC groups; Pearson correlation test, r = −0.175, t = −0.71, df = 16, P = 0.485; Fig. 5).

For the B. afzelii-infected nymphs (as determined by the RLB assay), 42.0% (81/193) of the infections contained exotic major ospC groups. Conversely, for the B. garinii-infected nymphs (as determined by the RLB assay), 55.3% (105/190) of the infections contained exotic major ospC groups. In B. afzelii, 1.88% of the ospC gene sequences (2,150/114,432) were of exotic origin. In B. garinii, 3.09% of the sequences (3,897/125,978) were of exotic origin.

The five native B. afzelii major ospC groups that had undergone horizontal transfer (to B. garinii) accounted for 84.62% of the native ospC gene sequences of B. afzelii (95,598/112,596). Similarly, the seven native B. garinii major ospC groups that had undergone horizontal transfer (to B. afzelii) accounted for 91.45% of the native ospC gene sequences of B. garinii (111,641/122,081). For the two Borrelia species combined, we compared the mean relative frequency between native major ospC groups that had experienced horizontal transfer or not (Tables 2 and 3). The 12 native major ospC groups that had experienced horizontal transfer had a significantly higher relative frequency (t = 4.89, df = 19, P < 0.001; Fig. 3) than the 9 native major ospC groups that had not been transferred. For the 12 major ospC groups that had experienced horizontal transfer, the relative frequency in the native species was always higher than the frequency in the exotic species (paired t test, t = 8.13, df = 11, P < 0.001). There was a significant correlation between the relative frequencies of these 12 major ospC groups in their native species and in their exotic species (Pearson correlation test, r = 0.783, t = 3.99, df = 10, P = 0.002).

For a subsample of 120 Borrelia isolates, we obtained high-quality sequences for both the recA gene and the hbb gene for 110 isolates. This genetic analysis confirmed the Borrelia species classification of the RLB assay in 95.5% (105/110) of our samples. In general, this additional analysis showed that the RLB assay provided a reliable classification of the Borrelia species.

Cross-reactive acquired immunity should reduce the co-occurrence of antigenically similar strains in the same reservoir host (5, 73, 74). Assuming that the genetic similarity between major ospC groups predicts the probability of cross-reactive acquired immunity between their respective OspC antigens, we expected a positive relationship between the pairwise genetic distance and the pairwise association index. A study of multiple B. afzelii infections in a population of wild voles found that dissimilar major ospC groups were more likely to be found together in the same host than similar major ospC groups (1). In contrast, we found no relationship between the pairwise genetic distance and the pairwise association index inside the ticks. A simple explanation for the difference between these two studies is that one looked at the reservoir host, whereas the other looked at the tick vector. The vertebrate immune system structures the community of ospC strains inside the reservoir hosts, and the tick is expected to acquire a subsample of this strain community (13, 75). However, the community structure of the colonizing strains might be modified by subsequent events inside the tick. The spirochete load can change dramatically following the blood meal and during the larva-to-nymph molt, and these changes in abundance may alter the community of strains inside the tick (14, 76). There is currently very little information as to how Borrelia strains interact inside the tick over time (75).

Despite a large sequencing effort (240,410 sequences from 389 ticks over a period of 11 years), our study did not find any ospC alleles that were intermediately divergent (2% to 8%) from the community of major ospC groups. The average pairwise genetic distance between conspecific major ospC groups was always greater than 8%. Thus, our study confirmed the highly discontinuous pattern of genetic variation at the ospC locus, as was reported in previous studies of B. burgdorferi sensu lato (45–47, 49). The original studies demonstrating the 8% divergence cutoff between the major ospC groups in B. afzelii and B. garinii were based on a total of 140 sequences (46, 47, 49). Thus, it was reasonable to expect that a larger sequencing effort would discover the missing intermediate major ospC groups. The present study suggests that strains carrying ospC alleles that are 2% to 8% divergent rarely coexist in the same Borrelia species. Interestingly, our results also suggest that such intermediately divergent ospC alleles rarely coexist in Borrelia species occupying different niches in the vertebrate host community. The allospecific pairwise genetic distances between B. afzelii and B. garinii had the same discontinuous pattern, suggesting that the evolution of the ospC allele is not independent in each Borrelia species. One explanation is that frequent horizontal transfer of the ospC gene (see below) essentially synchronizes the community of major ospC groups between these two Borrelia species. Thus, a new ospC mutant can rise to a detectable frequency only if it is at least 8% divergent from all the major ospC groups of all the Borrelia pathogens in the local community.

Purifying selection by cross-reactive antibodies is the most plausible explanation for the missing intermediate major ospC groups. Mutant major ospC groups that are genetically similar and thus antigenically similar to their conspecific counterparts are vulnerable to cross-reactive antibodies, which prevent these mutants from invading the local population (11, 74, 77). Indirect evidence for this hypothesis comes from studies of B. burgdorferi sensu stricto that have shown that antibodies against a given OspC antigen provide protection against strains belonging to the same major ospC group but not against strains belonging to a different major ospC group (42, 44). Thus, in B. burgdorferi sensu stricto, we see the end result of divergent selection on the ospC antigen to escape from the cross-reactive antibodies directed against its neighbors. However, immune escape between strains carrying different major ospC group alleles is not always perfect. A recent study of B. afzelii found that OspC antigens belonging to different major ospC groups (that were 20.68% divergent in DNA sequence) could still interfere with one another via the cross-reactive antibodies of the rodent host (43). In summary, cross-reactive antibodies are the most plausible source of purifying selection that removes intermediately divergent major ospC groups from the local population of Borrelia pathogens.

Each major ospC group contains a cluster of ospC alleles that are less than 2% divergent from each other with respect to their DNA sequence. An earlier study of the ospC polymorphism found that there was very little sequence variation within the major ospC groups (9 of the 11 major ospC groups were isosequential) and concluded that the effective population size of each group was very small (45). The present study found that strain richness within the major ospC groups was high (331 unique seed sequences clustered into 23 major ospC groups, giving a mean richness of 14.4 seed sequences per group). Examination of the sequence variation within the major ospC groups found that the majority (∼60%) of the nucleotide substitutions were synonymous. Sequencing error was the likely cause of all insertions and deletions that resulted in a frameshift, as the downstream parts of the OspC protein were highly conserved. The present study suggests that most of the sequence variation observed within a major ospC group does not affect the amino acid sequence and/or the function of the OspC protein. Under this scenario, natural selection cannot discriminate among the different allelic variants, and the frequencies of these silent alleles depend on genetic drift alone (78). Thus, the observation that the major ospC groups contain allelic variants does not contradict the hypothesis that cross-immunity selects against insufficiently divergent ospC alleles.

We found that the B. afzelii and B. garinii isolates in our local area contained 10 and 11 native major ospC groups, respectively (strain richness increased to 19 and 18, respectively, after the exotic major ospC groups were included). These findings were similar to those of previous studies of the ospC polymorphism in local tick populations (45, 53, 55, 79). In the United States, local populations of B. burgdorferi sensu stricto have between 13 and 16 different major ospC groups (13, 45, 53, 79). In Europe, local populations of B. afzelii have 10 different major ospC groups (55). Why do populations of B. burgdorferi sensu lato contain 12 and not 120 different major ospC groups? The relatively small length of the ospC gene (∼600 bp) in combination with the constraints of remaining functional and of avoiding cross-reactive acquired immunity with all of its neighbors must set an upper limit to the number of major ospC groups that can exist in a single, local population.

Ixodes ricinus nymphs infected with B. afzelii or B. garinii contained a mean richness of 2.37 and 2.97 different ospC strains per nymph, respectively. Our estimates of the mean ospC strain richness in questing I. ricinus nymphs were considerably higher than those from previous studies (1.07 to 1.44 strains per nymph) of B. burgdorferi sensu lato infections of questing nymphs in I. scapularis, I. pacificus, and I. ricinus (53, 55, 80). However, our estimates of the mean ospC strain richness are similar to those from previous studies (1.0 to 3.7 strains per tick) that sampled larval ticks feeding on infected hosts (13, 55, 81, 82) or questing adults (45, 53). We found a very high prevalence of multiple ospC strains in both B. afzelii (78.76%) and B. garinii (84.74%). Our estimates of the prevalence of multiple ospC strain infections were much higher than those in previous studies (0.0% to 50.0%) (1, 45, 55, 81–85). These earlier studies estimated the mean ospC strain richness and the prevalence of multiple strain infections by using molecular methods, such as the RLB assay, cold single-strand conformation polymorphism (SSCP) analysis, and Sanger sequencing. A previous study showed that the RLB assay was more sensitive at detecting infections with multiple strains than SSCP analysis (45, 53). Thus, differences in strain richness among studies are partially confounded by differences in the sensitivity of the molecular methods. In the present study, we used next-generation sequencing (NGS) to analyze an average of 618 ospC gene sequences per tick (240,410 sequences, 389 ticks). This high level of coverage for each tick allowed us to detect major ospC groups that represented a tiny fraction of that particular amplicon (<0.2%). Thus, NGS methods estimate higher levels of strain diversity than the older molecular methods because these approaches are more sensitive at detecting rare strains.

The ospC strains were highly aggregated in I. ricinus ticks. There are a number of explanations for this pattern of aggregation. First, vertebrate hosts vary in their competence in acquiring infections and transmitting spirochetes to vector ticks (18, 30, 86, 87). If a substantial proportion of ticks feed on incompetent hosts and the remaining ticks feed on infected reservoir hosts, the Borrelia strains will appear to be aggregated inside the ticks when these two subpopulations of ticks are analyzed together. Second, the aggregation of ticks on reservoir hosts can also cause the aggregation of Borrelia strains in the tick vector. Numerous studies have observed that a minority of wild rodent hosts feeds a majority of the I. ricinus ticks (88, 89). There are a variety of ecological factors that cause aggregation of ticks in wild rodents (90). Third, the immune system of the vertebrate host might have more difficulty controlling and/or clearing infections with multiple strains than infections with a single strain (74, 91). Under this hypothesis, infections with multiple strains have a higher probability of establishing infection in the reservoir host and/or have a higher host-to-tick transmission success than infections with a single strain.

The mean ospC strain richness of B. garinii (2.96 strains per nymph) was 25% higher than that of B. afzelii (2.37 strains per nymph). Borrelia garinii ospC strains were more often found in multiple infections and had higher pairwise association indices than their B. afzelii counterparts. The higher mean ospC strain richness for B. garinii was not caused by differences in total strain richness, which was the same between B. afzelii (n = 19 ospC strains) and B. garinii (n = 18 ospC strains). Borrelia garinii is mostly found in avian hosts, whereas B. afzelii is mostly found in rodent hosts (16, 20–22, 24, 25). Thus, one explanation is that the richness of B. garinii ospC strains in birds is higher than the richness of B. afzelii ospC strains in rodents. To date, most of the studies investigating ospC strain richness in vertebrate reservoir hosts have focused on mammals (1, 13, 50, 52, 55, 66, 82, 83). The two studies that investigated birds (81, 85) found a mean ospC strain richness (1.00 to 1.76 strains per tick) that was low compared to that found in the studies with mammals, but this could be an artifact of the relatively inefficient molecular methods that were used in the bird studies (RLB assay and Sanger sequencing). A recent study characterized the spirochete loads of different Borrelia species inside questing I. ricinus nymphs (14). This study found that the median spirochete load of B. garinii (5,080 spirochetes per nymph) was 1.6 times higher than that of B. afzelii (3,140 spirochetes per nymph) (14). Thus, another explanation for the discrepancy in strain richness between the two species is that the tick vector has a higher carrying capacity for B. garinii than for B. afzelii.

In the present study, we used only nymph-derived isolates that, according to the RLB assay (which targeted the 23S-5S spacer gene), were singly infected with either B. afzelii or B. garinii, and we therefore excluded those isolates that, according to the RLB assay, were coinfected with multiple Borrelia species. A surprising result was therefore that 42% of the B. afzelii-infected ticks and 55% of the B. garinii-infected ticks contained exotic major ospC groups from other Borrelia species. While the exotic major ospC groups were widespread, they were not abundant. For B. afzelii and B. garinii, the exotic major ospC groups accounted for 1.88% and 3.09% of all the ospC sequences, respectively. Two alternative explanations for these results are the coinfection explanation and the horizontal gene transfer explanation.

In the coinfection explanation, the presence of native and exotic ospC groups in the same isolate represents a true coinfection with B. afzelii and B. garinii that was not detected by the RLB assay. In this scenario, the RLB assay detected the highly abundant Borrelia species (corresponding to the abundant native ospC alleles from the 454 sequencing) but not the Borrelia species of low abundance (corresponding to the low-abundance exotic ospC alleles from the 454 sequencing). Although we have recently shown that the sensitivity of our RLB assay is similar to that of a widely used quantitative PCR assay for detecting Borrelia infections in ticks (14, 92), we cannot exclude this possibility. Previous studies have shown that nymphs and adult ticks carry coinfections with B. afzelii and B. garinii (14, 32–34). In our local I. ricinus population, the RLB assay found that 0.75% (13/1,731) of the nymphs are coinfected with B. afzelii and B. garinii (14). If the coinfection explanation is true, the implication is that the true frequency of coinfections with B. afzelii and B. garinii has been underestimated by a factor of 56 to 73 (42% to 55% in the present study versus 0.75% in the study by Herrmann et al. [14]). Interestingly, we recently showed that the spirochete load of B. afzelii and B. garinii coinfections was much lower than the additive expectation (14). We suggested that the complement system of the vertebrate host reduced the spirochete load of the maladapted Borrelia species inside the nymphal tick (14). Thus, the observation that the exotic ospC groups had low abundances (1.88% and 3.09%) is consistent with our understanding of how the vertebrate complement system would interact with coinfections containing bird- and rodent-adapted Borrelia species inside the same tick.

If the coinfection hypothesis is true, the observation that >50% of our nymphs had coinfections with B. afzelii and B. garinii requires an additional explanation, given that these two Borrelia species are adapted to different vertebrate hosts (24, 32). We recently presented four mechanisms by which B. afzelii and B. garinii could coinfect the same tick (14). The specificity of B. afzelii for rodents and of B. garinii for birds is not 100% complete. Studies in France (93) and England (94) found that introduced species of rodents were infected with the bird-adapted B. garinii. Similarly, other studies found that birds can transmit the rodent-adapted B. afzelii (95, 96). Rodents or birds that are coinfected with B. afzelii and B. garinii could transmit both species to feeding ticks. Larval ticks could also acquire coinfections of B. afzelii and B. garinii by taking multiple blood meals from different hosts (43, 44), by a combination of cofeeding and systemic transmission, or by a combination of vertical transmission and systemic transmission (14). For example, a larval tick cofeeding next to a B. afzelii-infected nymph on a B. garinii-infected bird would acquire both Borrelia species. Cofeeding transmission has been observed in both B. afzelii and B. garinii (90). In contrast, vertical transmission of Borrelia pathogens is believed to be rare (36).

In the horizontal gene transfer explanation, the RLB assay result is correct (all nymphs were infected with a single Borrelia species) and the exotic major ospC group alleles were horizontally transferred into the recipient Borrelia species. Previous studies have found that major ospC group alleles can be transferred between Borrelia species (46, 97, 98). In the present study, we found that 12 of the 21 native ospC major groups occurred as exotic major ospC groups in the other Borrelia species. If the horizontal gene transfer explanation is true, we make the following three observations. First, the present study suggests that horizontal transfer of the ospC gene is much more common than previous reports in the literature would suggest (45, 72, 97–100). Second, horizontal gene transfer was more likely for major ospC groups that were common in the donor Borrelia species. Third, some of the exotic ospC alleles increased to an appreciable frequency in the recipient Borrelia species. For example, the major ospC group A10, which is native to B. afzelii, had a relative frequency of 31.58% (60/190 infected nymphs) in B. garinii (Table 3). If the horizontal gene transfer explanation is true, one puzzling observation is why the exotic major ospC groups are so much less abundant (1.88% to 3.09%) than the native ospC groups (98.12% and 96.91%) in the nymph-derived isolates. This observation suggests that the exotic strains are underperforming relative to the native strains. Future studies should perform deep sequencing of other Borrelia genes that differ between B. afzelii and B. garinii (e.g., recA or hbb) in order to differentiate between the coinfection hypothesis and the horizontal gene transfer hypothesis.

In conclusion, despite a large sequencing effort (240,410 sequences), our study was unable to find any ospC alleles that were intermediately divergent (>2% to <8%) from the community of major ospC groups. Our study suggests that cross-reactive antibodies prevent the invasion of intermediately divergent ospC alleles. In contrast, we found no evidence that cross-immunity in the vertebrate host was structuring the community of ospC strains in the tick vector. Deep sequencing found that more than 50% of our Borrelia isolates contained exotic major ospC group alleles that had a low abundance (1.88% to 3.09% of all ospC gene sequences) and that had been derived from other Borrelia species. Two alternative explanations for these exotic major ospC group alleles are coinfections, where only the numerically dominant Borrelia species was detected by the RLB assay, or horizontal transfer of the ospC gene between Borrelia species.