Elucidating the factors that shape species distributions has long been a fundamental goal in ecology and evolutionary biology. In spite of significant theoretical advancements, empirical studies of range limits have lagged behind. Specifically, little is known about how the attributes that allow species to expand their ranges and become widespread vary across phylogenies. Here, we studied the ascidian Botryllus schlosseri, a worldwide invasive species that is also characterized by marked genetic subdivision. Our study includes phylogenetic and population genetic data based on mitochondrial and nuclear genes, as well as polymorphic microsatellites for B. schlosseri colonies sampled from the southern and northern coasts of Europe and the eastern and western coasts of North America. We demonstrate that this well-known model organism comprises three highly divergent and probably reproductively isolated cryptic species (A, D and E), with two more (B and C) being suggested by data retrieved from GenBank. Among these, species A, recovered in all of the surveyed regions, is by far the most common and widespread. By contrast, species B–E, occurring mostly in sites from northern Europe, are considerably more geographically restricted. These findings, along with inferences made on transport opportunity, suggest that divergent evolutionary histories promoted differences in invasive potential between B. schlosseri sibling species, indicating that attributes that facilitate dramatic shifts in range limits can evolve more easily and frequently than previously thought. We propose environmental disturbance as a selective force that could have shaped the evolution of invasiveness in the B. schlosseri complex.
A fundamental goal in ecology and evolutionary biology is to understand the factors that set the limits of species ranges . Current perspectives on range limits are integrated in a wide variety of theoretical models that distinguish between the presence or absence of environmental heterogeneity and/or evolution . Models without environmental gradients or evolution consider Allee effects, whereby ranges are restricted by the inability of marginal populations to grow when small . Models with environmental heterogeneity but with no evolution incorporate dispersal barriers or fundamental niche limitation along ecological gradients . Finally, models with environmental gradients and evolution consider the effect of gene flow and selection on adaptation to discrete habitats at a range's edge [5,6]. While considerable attention has been directed towards the theory of range limits, empirical advancements have lagged behind. A key aspect that remains to be extensively addressed is the degree to which range limits are phylogenetically conserved (i.e. do closely related lineages have similar distributional limits?).
An ideal opportunity to test the theory and patterns of range limits is offered by biological invasions. These can be viewed as ‘natural experiments’ during which species are often placed in settings outside their recent evolutionary context, where some spread extensively while others remain spatially restricted or fail to become established altogether. Early surveys across a large spectrum of taxa revealed that invasive species are often clustered within certain families like ducks (Chordata, Anatidae) or grasses (Magnoliophyta, Poaceae ), indicating that characteristics responsible for the invasive status of some species are highly conserved, transcending deep phylogenetic divisions. Recently, however, notable differences in the degree of spatial spread and invasion success have been reported at lower taxonomic levels between sister species [8–10], suggesting that the ability of taxa to spread and become invasive once introduced to a new environment might be a more variable attribute than previously thought.
Here, we explore these issues in the colonial ascidian Botryllus schlosseri Pallas 1766 (Chordata: Ascidiacea). Also known as the golden star tunicate, B. schlosseri is a powerful model system in the fields of immunobiology, developmental biology and evolutionary biology [11,12]. It is generally considered of European origin [13,14] (but see ), and was probably introduced via European shipping to the east coast of North America in the early 1800s . Since then, it has spread to coastal waters of all non-polar continents, with new records of established populations being added frequently . Importantly, although B. schlosseri is one of the most widespread marine invaders, it also presents the characteristics of a species likely to experience significant genetic breaks across small spatial scales and strong local adaptation. These include a sessile adult stage, a brief free-swimming larval stage  and gregarious settlement of kin larvae . Consistent with this expectation, previous genetic surveys have documented high levels of genetic differentiation between B. schlosseri populations [16,19,20]. Despite these observations, the geographical and/or biological barriers to gene flow within this taxon remain largely unexplored. Moreover, little is known about how B. schlosseri genetic groups differ in the extent of their geographical distribution and invasion success.
We characterize the phylogenetic and population genetic structure of European and North American populations of B. schlosseri, using the mitochondrial cytochrome c oxidase subunit I (COI), the nuclear 18S rRNA and 10 polymorphic microsatellites. We aim to (i) determine the extent of genetic subdivision within this species, (ii) evaluate the degree of independence of genetic groups (i.e. do identified genetic breaks delineate unrecognized cryptic species?), and (iii) examine the spatial distribution of genetic assemblages/sibling species and determine their invasive potential. Finally, we interpret these findings in the light of evolutionary mechanisms known to enhance the ability of species to spread and become invasive.
2. Material and methods
(a) Sampling and molecular methods
Botryllus schlosseri colonies were sampled from the southern and northern coasts of Europe, and the eastern and western coasts of North America (figure 1 and table 1). DNA extractions were performed as described by Lejeusne et al. . A fragment of COI was amplified using the LCO1490/HCO2198 primers . Failed amplifications were repeated using the species-specific primer BsCOIR (5′-GTATTTTATTTTTAGAATTTGGTCAAG-3′) and HCO2198. For phylogenetic purposes, we used 24 additional COI sequences: 19 from GenBank [19,22] and 5 from unpublished data (A. Lacoursière-Roussel 2011, unpublished data; electronic supplementary material, appendix 2). A subset of 42 individuals from the major COI clades was analysed for the nuclear 18S rRNA gene, using the 18S1/18S4 primers . Reaction chemistry and cycling parameters followed Bock et al.  for COI and Xu et al.  for 18S rRNA. Sequences were aligned in CodonCode Aligner v. 2.0.6 (CodonCode Corporation, Dedham, MA).
Nuclear genetic variation was further investigated using 10 polymorphic microsatellites: BS321 , PB29, PB49, PB41 and PBC1 , and Bsm1, Bsm2, Bsm4, Bsm6 and Bsm9 . Amplification conditions, fragment visualization and scoring followed Bock et al. . For samples that failed to amplify, we used a modified protocol (5 cycles at 40°C annealing, and 35 cycles at 43°C annealing).
(b) Data analysis
Phylogenetic relationships were examined using neighbour-joining (NJ) analysis in MEGA v. 4 , and Bayesian inference (BI) in MrBayes v. 3.1 . For detailed methods see electronic supplementary material, appendix 1. Relationships among COI haplotypes were also investigated using a haplotype network in TCS v. 1.21 . For each population, COI haplotype and nucleotide diversities were calculated in Dnasp v. 4.50 . Additionally, the COI data were used to calculate population pairwise ΦST values, with 104 permutations in Arlequin v. 3.0 , using the Tamura–Nei substitution model . Significance levels were adjusted using sequential Bonferroni corrections .
Microsatellite data were checked for deviations from Hardy–Weinberg equilibrium (HWE) and linkage disequilibrium (LD) with Genepop v. 4.0 . Genetic diversity indices of allelic richness and gene diversity were calculated in FSTAT v. 18.104.22.168 . Genetic relationships between multilocus genotypes were assessed using factorial correspondence analysis (FCA) in GENETIX v. 4.05 . Genetic structure was investigated in Arlequin, by calculating pairwise FST values with 104 permutations and significance levels of p-values adjusted as above. To infer population structure without predefined population subdivision, we used the Bayesian clustering approach implemented in Structure v. 2.3.2 . For detailed methods see electronic supplementary material, appendix 1.
(a) Sequence variation and phylogenetic reconstructions
The COI alignment was 524 bp long, and contained 160 polymorphic sites, of which 142 were parsimony informative. Within 586 sequences, we identified 60 haplotypes: 19 new (GenBank accessions: JN248359–JN248377) and 41 reported previously ([19,20,22]; A. Lacoursière-Roussel 2011, unpublished data). The 18S rRNA alignment was 817 bp long and contained 82 polymorphic sites, of which 81 were parsimony informative. Five unique 18S rRNA sequences were identified among 42 individuals (GenBank accessions: JN248378–JN248382).
The COI phylogenetic reconstructions (identical for NJ and BI algorithms) recovered five strongly supported monophyletic clades (clades A–E; figure 2), with high interclade sequence divergences (10.8–16.5%) and comparatively lower intraclade distances (0.8–3.8%). While members of clades A, D and E were sampled in this study, clades B and C consisted of haplotypes retrieved from GenBank (figure 2; electronic supplementary material, appendix 2). The TCS network also recovered five major groups of COI haplotypes. Within each group, considerable genetic structure was detected, with haplotypes separated by up to 47 mutation steps (electronic supplementary material, appendix 3). The 18S rRNA phylogeny was congruent with that based on COI in recovering with high support the monophyly of the three lineages sampled in this study (figure 2), with high divergences (2.3–10.1%). No clear differences in colour morphs were observed between members of different phylogenetic clades.
(b) Geographical patterns of cytochrome c oxidase subunit I diversity
Clade A was geographically widespread. It characterized most (98.5%) individuals from southern Europe and all of the ones from North America. It was also present in five sites from northern Europe (figure 1 and table 1). Clades B and C, not recovered here, were previously sampled in two European locations on the coast of Spain . Clade D, also geographically restricted, was recovered in one site in the English Channel. Similarly, Clade E was mainly recovered in six populations concentrated in the English Channel. Only three individuals of clade E were sampled outside this region, in two southern European sites (figure 1 and table 1). The geographical distribution of clades was largely allopatric, with individual populations containing only one major lineage. However, in six of the sites surveyed, clades A and E were found in sympatry, with one of the two lineages generally dominating each site (table 1). These sympatric locations provided us with a unique opportunity to investigate the level of contemporary gene flow (and indirectly, reproductive isolation) between members of the two divergent clades using microsatellite data in Structure (see below).
Genetic diversity indices of haplotype diversity (h) and nucleotide diversity (π) computed for clade A populations were higher for southern European sites (h = 0.511; π = 0.011) than for those in northern Europe (h = 0.165; π = 0.004), or eastern (h = 0.254; π = 0.008) or western North America (h = 0.323; π = 0.008). However, this general trend masked substantial variation, as all surveyed regions contained monomorphic as well as highly diverse populations. Pairwise ΦST values computed for clade A populations also revealed heterogeneity in genetic structure, with instances of high and significant differentiation between sites within the surveyed regions, as well as low and non-significant differentiation between sites at macrogeographical scales (electronic supplementary material, appendix 4). For clade D, all individuals shared the same haplotype, while for clade E all populations showed high levels of genetic variation (h = 0.605; π = 0.005; table 1). With two exceptions, all pairwise ΦST values computed for clade E populations were high and significant (electronic supplementary material, appendix 5).
(c) Geographical patterns of microsatellite diversity
For the widely distributed clade A, we discerned 135 alleles among the 10 microsatellite markers assayed (mean: 13.5 alleles/locus). While most loci conformed to HWE expectations, 36 out of 156 cases showed significant deviations from HWE. However, no systematic deviations were observed across loci or populations (electronic supplementary material, appendix 6). Also, there was little evidence for LD, with only 6 of 720 locus pairs remaining significant after Bonferroni corrections in one of the sites surveyed. Genetic diversity indices of allelic richness (A) and gene diversity (HE) computed for clade A populations were generally higher for southern European sites (A = 5.2; HE = 0.636) than for those in northern Europe (A = 3.0; HE = 0.520), or eastern (A = 3.2; HE = 0.503) or western North America (A = 3.3; HE = 0.560; for detailed clade A microsatellite results see table 1; electronic supplementary material, appendix 6). The FCA computed using microsatellite data from clade A revealed weak clustering of genotypes based on geographical origin. Higher genetic similarities were sporadically revealed for individuals sampled between rather than within biogeographic regions (figure 3). This result was also evidenced by the low and non-significant pairwise FST values occasionally recorded between sites at macrogeographical scales (electronic supplementary material, appendix 7). For the divergent clades D and E, consistent PCR amplifications were obtained for 7 of the 10 microsatellite loci. However, the genotypic profiles of these individuals were characterized by a low number of alleles (mean: 6.3 alleles/locus) and a high incidence of homozygotes. Since these features are generally indicative of datasets with high rates of null alleles, microsatellite data for clades D and E were used only for the Structure analysis, which minimizes equilibrium requirements.
For the Structure analysis, when data from the seven microsatellite loci amplified in all samples from clades A, D and E were considered, the most likely model was the one with three clusters (electronic supplementary material, appendices 6 and 7). At K = 3, each cluster corresponded perfectly to the three major clades identified by COI and 18S rRNA (i.e. clades A, D and E; figure 4). Admixture coefficients ranged from 0.94 to 1.00, with all individuals showing 90 per cent probability intervals encompassing coefficients of 1, suggesting no contemporary gene flow between the three clusters (figure 4).
Two important observations emerge from these results. First, our findings show that B. schlosseri, one of the most widely researched and well-known ascidians, comprises at least three previously unrecognized and probably reproductively isolated cryptic species. Second, results presented here indicate that only one B. schlosseri species is widespread globally, while its sibling species are highly geographically restricted.
(a) Cryptic speciation in B. schlosseri, a model ascidian
Botryllus schlosseri was initially described over two and a half centuries ago, representing the first ascidian species ever recognized . Since then, it has experienced repeated taxonomic revisions owing to extensive colour polymorphism . As a result, multiple synonymies are available (Ascidiacea World Database, http://www.marinespecies.org/ascidiacea/), all of which are currently lumped under a single species designation. Our results oppose these assumptions, and demonstrate that B. schlosseri is a complex of evolutionarily distinct species, which are at least superficially morphologically indistinguishable. Both mitochondrial (COI) and nuclear (18S rRNA and microsatellite) markers consistently recovered three well-supported clades (i.e. clades A, D and E; figures 2 and 4), with two more being confirmed by COI data retrieved from GenBank (i.e. clades B and C; figure 2). The genetic divergences recorded between these lineages (10.8–16.5% COI divergence), corresponding to about 4.3–11.0 Myr of divergent evolutionary history (1.5–2.5% per Myr ), are significant. These levels far exceed COI divergences reported between cryptic species of other colonial tunicates, like Pseudodistoma crucigaster (2.12% ), Clavelina lepadiformis (5% ) or Pycnoclavella communis (8.55% ), and are comparable with those identified between widely recognized cryptic species of the solitary tunicate Ciona intestinalis (11.1–18.4% [10,46,47]). Furthermore, the complete lack of contemporary gene flow between clades, even in locations where their distribution intersects (figure 4), suggests that reproductive isolation is probably complete. Altogether, these findings indicate the need for a taxonomic re-evaluation of B. schlosseri.
Botryllus schlosseri sensu lato has been used extensively in pioneering research on the heritability of senescence and mortality, the evolutionary importance of self/non-self recognition or the genetic control of organogenesis (reviewed by Manni et al.  and Rosengarten & Nicotra ). Since accurate classification of model species is a major component for corroborating data drawn from different laboratories, the identification of cryptic speciation in B. schlosseri can be expected to impact a number of research fields. Beyond these considerations, our results highlight the importance of confirming species status via molecular data and/or crossing experiments. When such information is not available, species identifications should be interpreted with caution, even in well-known model organisms such as B. schlosseri.
(b) Contrasting distribution ranges in closely related species
Results presented here further indicate that B. schlosseri sibling species differ markedly in the extent of their geographical distribution and, probably, invasive potential. Species A is by far the most common and widespread, with a range spanning all biogeographic regions surveyed here (figure 1). Comparisons with previous molecular studies further indicate that this species is also present at additional sites in the eastern Mediterranean Sea , the east and west coasts of North America [49,50], South America , Japan, Australia  and New Zealand . Apart from the remarkable cosmopolitan distribution, additional lines of evidence suggest that species A experiences recurrent long-distance dispersal events (most probably human-mediated), and is highly invasive. These include frequent occurrence of distantly related haplotypes within populations (table 1; electronic supplementary material, appendix 3), extensive haplotype sharing among sites within biogeographic regions (table 1) and patchy distribution of genetic variation between sites at macrogeographical scales (figure 3 and table 1; electronic supplementary material, appendices 4 and 5). By contrast, species B and C, not recovered here, have only been found in two European locations on the coast of Spain . Similarly, with the exception of three individuals sampled in southern Europe, species D and E were only found in the English Channel (figure 1 and table 1). Although an important next step to confirm our findings would be to perform global sampling for B. schlosseri, particularly covering the centre of botryllid diversity in the western and southern Pacific Ocean , results presented here suggest that species A is much more geographically widespread than species B–E.
In colonial ascidians, because of the short-lived larval stage and sessile adult stage, the natural potential for range expansion is limited and is represented mainly by the rafting of individuals on floating debris . Although this mechanism of dispersal could have contributed to intracoastal spread within the surveyed biogeographic regions, it is unlikely that it also facilitated instances of dispersal at macrogeographical scales observed here, particularly for species A. Human-mediated dispersal is a more likely candidate, as colonial ascidians are known to be problematic foulers of ship hulls and aquaculture equipment . In this context, the sharp biogeographic asymmetries observed here between B. schlosseri sibling species can be explained by two scenarios that are not mutually exclusive: (i) stochastic events such as biases in human-mediated transport routes and opportunity facilitated by chance widespread introductions for B. schlosseri species A, but not for species B–E; or (ii) divergent evolutionary histories resulted in marked differences in the ability of B. schlosseri species to survive transport, establish, spread and become invasive. These scenarios are difficult to separate, but a number of clues suggest that the first one is less likely. First, biases in transport opportunity are expected to be more important in recent invaders, reflecting the stochastic nature of the invasion process. In species with long histories of anthropogenic transport like B. schlosseri, access to dispersal vectors should have eventually counterbalanced differences in spatial spread between sibling species of similar invasive potential. Second, the invasive species A co-occurs with species D and E in the English Channel, a region with extensive opportunity for transport via ship traffic. Third, species E individuals that are probably the results of recent introductions were occasionally found in southern European sites (figure 1 and table 1) . Collectively, these findings oppose the possibility that lack of transport prevented the widespread occurrence of B. schlosseri species other than species A.
The alternative scenario that divergent evolutionary histories shaped differential invasive potential in B. schlosseri sibling species is intriguing. The classical perception is that attributes determining range limits are highly evolutionarily conserved . This view is supported by the few available surveys reporting that closely related species have similar range sizes , and by the observation that invasive species are often clustered within certain families . Results presented here refute this earlier assumption, and add to a growing body of literature reporting marked differences in range limits and invasion success in closely related species [8–10]. Collectively, these studies suggest that species' range limits and invasion abilities are much more evolutionarily volatile than previously thought.
An ensuing question is what prompted the evolution of differential invasion abilities in B. schlosseri species. Much of the current literature on the evolution of invasiveness stems from studies of plant taxa [53,54]. Four main types of evolutionary change are thought to contribute to plant invasiveness: hybridization, polyploidy, bottlenecks and stress-induced genome modification (reviewed by Prentis et al. ). In the case of B. schlosseri, environmental disturbance appears to be the most likely scenario. Recently, it has been proposed that contrasting invasion abilities can evolve between populations or sibling species subjected to different levels of disturbance in their native range . Theoretical models also suggest that environmental fluctuations spanning different timescales can select for organismal flexibility and evolvability , both of which are favourable for invading new habitats . Indeed, owing to their sessile nature, colonial ascidians such as B. schlosseri are prone to be exposed to a multitude of environmental stressors. Although the ancestral distribution range and, implicitly, historical levels of disturbance are difficult to infer for B. schlosseri because of its long and complex history of anthropogenic range expansions, the disturbance scenario provides a promising working hypothesis that remains to be tested. Ecological genomics experiments can be used in future studies to reveal the extent of ecological tolerance for different cryptic species, quantify additive genetic variance for critical invasive traits, identify genes of major phenotypic effect and ultimately clarify the evolution of invasiveness in this system.
We would like to kindly acknowledge C. W. McKindsey, A. Lacoursière-Roussel as well as J. T. Carlton, R. P. Walter, A. Zhan and students in the Cristescu and MacIsaac laboratories for stimulating discussions. We are grateful to our colleagues J. Bishop, M. Carman, P. Chevaldonné, S. Dashfield, E. Grey, A. Izquierdo-Muñoz, D. Jollivet, C. Lejeusne, L. Lévêque, M. Litvaitis, L. Manni, C. Menniti, M. Molis, J. M. Nicolas, M. Rius, P. Sordino, C. Sorte, T. W. Therriault, X. Turon, B. Vercaemer, S. Williams as well as the Roscoff Biological Station Specimen Supply Service (CNRS/FR2424) for generously providing samples. This work was supported by the NSERC Canadian Aquatic Invasive Species Network (CAISN), Fisheries and Oceans Canada, NSERC Discovery Grants to H.J.M. and M.E.C., and by an Early Researcher Award from the Ontario Ministry of Research and Innovation to M.E.C.
- Received December 14, 2011.
- Accepted January 19, 2012.
- This journal is © 2012 The Royal Society