- Research article
- Open Access
"Nested" cryptic diversity in a widespread marine ecosystem engineer: a challenge for detecting biological invasions
BMC Evolutionary Biologyvolume 11, Article number: 176 (2011)
Ecosystem engineers facilitate habitat formation and enhance biodiversity, but when they become invasive, they present a critical threat to native communities because they can drastically alter the receiving habitat. Management of such species thus needs to be a priority, but the poorly resolved taxonomy of many ecosystem engineers represents a major obstacle to correctly identifying them as being either native or introduced. We address this dilemma by studying the sea squirt Pyura stolonifera, an important ecosystem engineer that dominates coastal communities particularly in the southern hemisphere. Using DNA sequence data from four independently evolving loci, we aimed to determine levels of cryptic diversity, the invasive or native status of each regional population, and the most appropriate sampling design for identifying the geographic ranges of each evolutionary unit.
Extensive sampling in Africa, Australasia and South America revealed the existence of "nested" levels of cryptic diversity, in which at least five distinct species can be further subdivided into smaller-scale genetic lineages. The ranges of several evolutionary units are limited by well-documented biogeographic disjunctions. Evidence for both cryptic native diversity and the existence of invasive populations allows us to considerably refine our view of the native versus introduced status of the evolutionary units within Pyura stolonifera in the different coastal communities they dominate.
This study illustrates the degree of taxonomic complexity that can exist within widespread species for which there is little taxonomic expertise, and it highlights the challenges involved in distinguishing between indigenous and introduced populations. The fact that multiple genetic lineages can be native to a single geographic region indicates that it is imperative to obtain samples from as many different habitat types and biotic zones as possible when attempting to identify the source region of a putative invader. "Nested" cryptic diversity, and the difficulties in correctly identifying invasive species that arise from it, represent a major challenge for managing biodiversity.
Biological invasions are a major global threat that can fundamentally and irreversibly modify native communities [1, 2]. Particularly when a biological invasion involves an ecosystem engineer, the consequences for an invaded ecosystem can be catastrophic . Ecosystem engineers monopolise space, accumulate biomass and have strong effects on species interactions by increasing architectural complexity of ecosystems and moderating environmental extremes . Non-indigenous species that function as ecosystem engineers are of major concern because they can replace indigenous habitat-forming species [5, 6] and drastically alter an invaded habitat [7, 8]. To maintain the diversity and integrity of biotic habitats, it is thus of great importance that such species are correctly identified and managed.
Coastal environments are among the most threatened ecosystems, with invasions of coastal assemblages across and between oceans facilitated by the movement of ocean-going ships and aquaculture [9–11]. However, in many of the world's coastal regions, a large proportion of marine species cannot be clearly identified as being either native or introduced due to a lack of systematic, biogeographic and historical evidence . The increasing availability of DNA sequence data has improved this situation to some extent, resulting in an exponential increase in the identification of cryptic biodiversity . Particularly in the case of poorly studied marine invertebrate groups, genetic methods can enable researchers to differentiate between recently introduced exotic species that should be monitored and controlled, and long-established, cryptic species that may have been previously overlooked and that may even require protection.
The ascidians (Chordata: Urochordata) are a group of sessile, filter feeding marine invertebrates that include both important ecosystem engineers and aggressively invasive species [14, 15]. Many ascidians are major occupiers of primary space along temperate coasts, where they provide habitat for numerous other organisms [16, 17] by enhancing habitat complexity when aggregated . Although ascidians have low natural dispersal potential because their lecithotrophic larvae remain in the plankton for very short periods of time (minutes to hours in most species) [19, 20], several species are recognised as pests on a global scale, occurring on multiple continents [21–23]. Dispersal on smaller scales may occur naturally as larvae attach themselves to floating objects that are moved around by currents, but adults attached to vessel hulls are considered to be the most likely vectors facilitating the worldwide spread of these species [24–26]. As for many other marine invertebrate groups, the taxonomy of some ascidians is poorly resolved [22, 27], and recent genetic studies have indicated that several supposedly cosmopolitan species are in fact comprised of two or more genetic lineages that should be treated as distinct species [28–32].
Here we examine the large, solitary ascidian Pyura stolonifera (Heller, 1878), which is an important foundation species particularly in temperate coastal regions of the southern hemisphere [16, 18, 33]. The taxonomic status and origin of the species are unclear, and have generated an extensive debate [34–38]. It remains uncertain whether populations in Africa, Australasia and South America are the fragmented remains of a pan-Gondwanan species [35, 36] or whether the species originated in one region and was recently been introduced to the other two regions [34, 39]. It is also disputed whether P. stolonifera is a single species  or a species complex . We study genetic patterns between regional populations to determine a) levels of cryptic diversity, b) the invasive or native status of each regional population and c) the most appropriate sampling design for identifying the boundaries of each evolutionary unit. Our findings indicate the presence of multiple genetic lineages within regions, which, together with inadequate sampling, can seriously hinder our capacity to detect invasive populations.
Pyura stolonifera is particularly common in southern Africa and Australia , but localised populations have also been reported from South America [16, 40], northwestern Africa [37, 38] and, most recently, New Zealand . The taxon is an important ecosystem engineer that dominates intertidal and subtidal habitats in Africa [17, 37, 41], Australia and New Zealand , and intertidal areas in Chile , where it achieves among the highest biomasses ever reported in such environments . Pyura stolonifera forms extremely large aggregations, resulting in aggressive monopolisation of the available substratum .
Sampling and amplification of genetic markers
A total of 518 ingroup samples were collected in all regions from which there are reliable reports of Pyura stolonifera, except Senegal. Within each region, samples were collected at several sites that span the taxon's entire range, including 16 sites in Africa, 26 sites in Australia, seven sites in New Zealand and one site in Chile (Table 1). A small piece of mantle tissue (< 1 cm3) from each individual sampled was preserved in a solution containing 70% ethanol and 30% TE buffer. This medium was replaced on a daily basis until it no longer changed color and until the tissue had become completely white. Obtaining high quality DNA proved difficult, and even an extraction protocol developed to eliminate contaminants present in ascidian tissues and tested specifically on Pyura stolonifera  did not produce better results than standard extraction protocols. We consequently used a salting-out protocol to extract DNA .
We amplified one mitochondrial DNA (mtDNA) marker, the cytochrome oxidase subunit I gene (COI), and three nuclear DNA (nrDNA) markers: 18S (a component of the 40S cytoplasmic small ribosomal subunit in eukaryotes) and two nuclear genes containing introns, namely ATP synthase subunit α (ATPSα) and Adenine Nucleotide Transporter (ANT, also known as ADP/ATP translocase) (Table 2). We also used unpublished COI sequence data generated previously by some of our collaborators, and incorporated some published sequence data (Table 1). The COI gene was the primary marker used for identifying genetic lineages due to its high variability. ANT was the most variable nuclear marker. It was primarily used to confirm genetic structure identified using COI by amplifying a sub-set of samples, to study genetic diversity in selected populations, and to provide an alternative for COI in the few cases where this marker did not amplify due to a possible mutation in the primer annealing region (two samples from Tasmania and all samples from New Zealand, Table 1). ATPSα and 18S were less informative and for that reason were only used for phylogeny reconstructions.
Most of the primers used are universal. The ANT gene did not amplify readily using published primers , and a forward primer was designed to amplify it in Stolidobranchia ascidians (StolidoANT-F) in conjunction with a universal reverse primer (Table 2). This primer combination proved particularly useful for both phylogenetic and phylogeographic work, as the PCR product amplified reliably and contained a long, variable intron. The primer combination developed here amplified the ANT gene not only in Pyura spp., but also in other genera within the order Stolidobranchiata, including Botrylloides (GenBank accession number JF962229), Botryllus (JF962231) and Styela (JF962232).
PCR reaction conditions comprised 1 μl of template DNA (~150 ng), 3 μl of reaction buffer (Promega), 6 μl of dNTP mixture containing 125 mM of each dNTP, 1.2 μl of each primer (5 mM dilutions), 1 unit of Taq DNA polymerase (Promega, Madison, USA) and ddH20 to a final volume of 30 μl. Concentrations of MgCl2 differed for each marker (Table 2). PCR profiles consisted of an initial denaturing step (94°C for 3 min), 35 cycles of denaturing (94°C for 30 s), annealing (45 s at a primer-specific annealing temperature, Ta; Table 2) and extension (72°C for 45 s), and a final extension step (72°C for 10 min). The problem of PCR reactions being affected by low purity of DNA extractions could be circumvented by diluting DNA templates, which supposedly reduced potential inhibitors to a level at which they no longer interfered with the PCR reaction. Nonetheless, a quality screening procedure was applied in which samples were excluded when the quality of their trace files was too low to identify each nucleotide with absolute certainty in three successive sequencing runs. As ANT tended to amplify more reliably than COI, this explains why at some sites, more ANT sequences than COI sequences were generated, even though only a fraction of samples was sequenced using this marker (Table 1). PCR products were purified using the UltraCleanTM 15 DNA Purification Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA), sequenced in both directions using Big Dye terminator chemistry version 3.1 (Applied Biosystems, USA) and run on a 3130xl Genetic Analyser.
Phase determination and alignment
In heterozygous individuals whose two ANT alleles differed in length, we used CHAMPURU v1.0  to determine each phase. In cases where there were no length differences, each sequence was deduced using default settings for multi-allelic loci without stepwise mutation in PHASE v2.1 . Whenever there were multiple possible phases, we selected the two alleles having the highest probability, which tended to be an order of magnitude greater than the probabilities of all other sequence pairs. Using alternative phases made no obvious difference in terms of estimates of genetic diversity and phylogenetic reconstructions.
COI and 18S sequences were aligned by eye in MEGA4 . The nuclear genes containing introns were aligned using the multiple sequence alignment program T-COFFEE at the BIOHPC website (http://cbsuapps.tc.cornell.edu). Poorly aligned regions were eliminated using GBLOCKS by specifying the least stringent conditions at the GBLOCKS server (http://molevol.cmima.csic.es/castresana/GBLOCKS_server.html).
Phylogenetic relationships among lineages of Pyura stolonifera were reconstructed using an aligned data-set from combined sequence data of four loci, totaling 2611 bp in length. After exploring the phylogenetic signal of data-sets from each locus, phylogenetic trees based on combined sequence data were reconstructed using no more than two representatives per lineage and region (14 ingroup taxa and four outgroup taxa). In several cases, regional lineages identified were excluded from the combined analyses because sequence data were not available for all four loci (Table 1). Three methods of phylogenetic reconstruction were employed and results compared: minimum evolution and parsimony, both employed in MEGA4, and Bayesian inference employed in MRBAYES 3.1 . Support for nodes in the minimum evolution and parsimony trees was assessed by means of 10 000 bootstrap replications. In Bayesian inference, four chains of three million generations each were run simultaneously and trees were sampled every 100 generations. After removing the first 10% of trees as burn-in, posterior probabilities of nodes were assessed by constructing a 50% majority rule consensus tree. To check for consistency of results, the analyses were repeated three times. For Bayesian inference, the data-set was divided into 14 partitions: codon positions 1-3 of COI (534 bp), 18S (1703 bp), codon positions 1-3 for the ANT exon (81 bp) and the ATPSα exon (105 bp), ANT intron (137 bp), ATPSα intron (51 bp), ANT indels (19 characters) and ATPSα indels (8 characters). Rates were allowed to vary among nucleotide partitions, and the GTR+I+Γ model was specified for each. Information from indels was only incorporated when these that had clearly defined alignment boundaries.
In several cases, lower-level phylogenetic relationships were inferred using sequence data from single loci by either constructing neighbour-joining trees  in MEGA4 using maximum composite likelihoods  of Tamura-Nei distances  or by constructing median-joining haplotype networks in NETWORK4516 (2009 version) . Indels were coded as single nucleotide differences irrespective of their length.
In populations that were not clearly differentiated on the basis of being comprised of monophyletic clades or clusters in haplotype networks, we estimated genetic diversity indices and calculated fixation indices to determine whether their allele frequencies differed. We used our two most variable markers, COI and ANT, to calculate the following statistics in ARLEQUIN v3.5 : h (haplotype or gene diversity), π (nucleotide diversity), and pairwise fixation indices as a measure of population differentiation (ΦST for COI and FST for ANT). In addition, we estimated observed and expected heterozygosity of ANT sequences in each population.
We used molecular dating to determine whether divergence of closely related populations that are represented in different regions likely occurred during historical times (i.e. as a result of a human-mediated introduction) or whether they have more ancient origins (i.e. divergence > 2000 years ago). Divergence times among several pairs of populations were estimated under the isolation-with-migration model  using the program IMa . We limited ourselves to populations that were either significantly differentiated on the basis of fixation indices, or that were comprised of reciprocally monophyletic sister clades. As two of the markers (ATPSα and 18S) showed little or no differentiation at the lowest taxonomic level, we used either a combination of COI and ANT, or COI only when ANT showed too little genetic variation to estimate divergence times. The population from New Zealand is represented by a single allele that was not found in its genetically most similar population in Australia (see Results) but is likely to be present there (see Discussion). For that reason, we did not estimate divergence times in this taxon. The HKY model  was selected as the best-fitting model for both markers by FINDMODEL, inheritance scalars of 0.25 for COI and 1 for ANT were specified, and a generation time of one year was assumed . To our knowledge, there are as yet no published evolutionary rates for the COI gene of ascidians, but evolutionary rates between 0.5 and 1.5% per million years (Myr-1) are assumed in most studies on other marine invertebrate taxa [61, 62]. To incorporate uncertainty, we specified a COI rate of 1% Myr-1 and incorporated a range from 0.5-1.5% Myr-1. The program does not require evolutionary rates to be specified for all markers included in an analysis. To check for the consistency of results, IMa runs were repeated three times using the following command line specifications: -l 25 000 (25 000 trees, with trees sampled every 100th generation, i.e. a total of 2.5 × 106 generations) -b 100 000 (deletion of the first 100 000 generations as burn-in), -q1 500 (maximum θ = effective population size parameter scaled to evolutionary rate), -t 2 (maximum divergence time scaled to evolutionary rate), -m1 20 (migration into population 1) -m2 20 (migration into population 2), -f g -n 80 -g1 0.999 -g2 0.3 (geometric heating scheme with 80 heated chains and heating parameters of 0.999 and 0.3).
Identification of evolutionary lineages
A total of 888 DNA sequences were generated, including 403 COI sequences, 411 ANT sequences, 34 ATPSα sequences and 40 18S sequences (Table 1). Most sequences (excluding some ANT sequences that were < 200 bp in length) were submitted to GenBank (accession numbers JF961754 - JF962415, see Additional File 1). Complete data-sets of aligned sequences are available in the following additional files: COI: Additional File 2; ANT: Additional File 3; ATPSα: Additional File 4 and 18S: Additional File 5. Additional File 6 contains combined sequence data used for phylogeny reconstruction. Six monophyletic clades were recovered with high support (bootstrap values and Bayesian posterior probabilities ≥ 99%, Figure 1b). Taxonomic descriptions exist for four of these clades (see Discussion for details). The African species Pyura stolonifera sensu stricto and P. herdmani formed a well-supported clade, and we also found strong support for the monophyly of two clades comprised of samples from Australasia and Chile (P. praeputialis and Pyura sp., Figure 1b). The placement of a third species present in Australasia (P. dalbyi, western and southeastern Australia) remains unresolved (Figure 1b). Support for its monophyly with the other two Australasian species was high in phylogenies of two of the more slowly-evolving partitions (e.g. ANT gene, Minimum Evolution Bootstrap support: 96%, 18S: 96%) and such a taxonomic placement is also supported by morphological data (see Discussion).
Regional sub-structuring was identified within Pyura herdmani and Pyura sp. on the basis of phylogenetic trees or haplotype networks constructed using sequence data from single markers. Pyura herdmani was comprised of four distinct lineages on the basis of mtDNA COI sequences (Figure 1d), although the more slowly-evolving nrDNA ANT sequences showed little differentiation among regions (Figure 1c), possibly due to incomplete lineage sorting. One of the lineages recovered using COI sequences is confined to Morocco, two occur in temperate South Africa, and the fourth is restricted to subtropical and tropical regions of southern Africa (south-eastern and eastern South Africa, and southern Mozambique). We included representatives of only two lineages in the phylogenetic tree based on combined sequence data (Figure 1b) to indicate that genetic differentiation among them is not much lower than among the taxa we considered to be distinct species, but it is important to note that there is presently not enough data to recognise any additional species within P. herdmani. The Australasian species Pyura sp. was comprised of two closely related ANT lineages of which one is found in Australia (Victoria, Tasmania and South Australia) and the other occurs both in Australia (Victoria and Tasmania) and in New Zealand (Figure 1e). Genetic diversity in Australia was high (13 unique alleles in 32 specimens), whereas all 42 individuals from New Zealand had the same ANT allele. We identified four heterozygous individuals in Australia having alleles from both lineages, suggesting that these are not distinct species.
Two pairs of geographically distant populations were not recovered as being distinct on the basis of phylogenetic trees or haplotype networks: Australian vs. Chilean representatives of Pyura praeputialis (Figure 1b) and Western Australian vs. southeastern Australian representatives of P. dalbyi (Figure 1f). Genetic diversity statistics were similar for all four populations of P. praeputialis investigated (three populations from Australia and one from Chile, Table 3). In most cases, genetic diversity estimates for the supposedly recently introduced population from Antofagasta, Chile, were the second highest of all the populations studied, and heterozygosity of this population at the diploid, intron-containing ANT gene was not lower than that of the Australian populations. On the basis of both pairwise ΦST values among the mtDNA COI haplotypes and pairwise FST values among the alleles of the nuclear ANT gene, we found significant structure between a site in southeastern Australia (Kilcunda) and all other sites, including the Chilean site Antofagasta. Significant structure between southeastern Australia and the two other regions, but not between the Australian east coast and Antofagasta, was also found when using a larger COI data-set that included samples from additional sites in eastern (sites 18 - 24 in Table 1) and southeastern (sites 28, 33 and 34) Australia for which no ANT data were generated (N = 232; southeast coast vs. east coast: ΦST = 0.070, P < 0.01, southeast coast vs. Antofagasta: ΦST = 0.070, P < 0.01; east coast vs. Antofagasta: ΦST = -0.003, P = 0.51). This supports the idea that the significant genetic structure found between the two eastern Australian sites and the one southeastern Australian site is not an artifact of small samples sizes.
No structure was identified between representatives of P. dalbyi from Western Australia vs. southeastern Australia (ΦST = 0.06, P = 0.18), and diversity indices were similar (Western Australia: h = 1.000 ± 0.045, π = 0.005 ± 0.003, southeastern Australia: h = 1.000 ± 0.017, π = 0.004 ± 0.003; COI sequence data only).
A divergence time estimate of 1.1 million years ago (95% confidence interval: 0.4 - 2.4 million years ago) was estimated for the Moroccan population of P. herdmani and its temperate southern African sister lineage (Figure 1c). Although the eastern and southeastern Australian populations of P. praeputialis shared haplotypes, they were also estimated to have diverged prior to the historical period (150 thousand years ago with a 95% confidence interval of 79 - 420 thousand years).
In the present study, we show that the widespread ascidian Pyura stolonifera is a species complex that comprises at least five distinct species. Within some of these, we found additional genetic structure at regional scales, and we identified three populations that are likely to be non-indigenous. The species associated with the species complex are ecosystem engineers that create habitat complexity, compete with other sessile species for food and space and tend to be highly abundant once established [15, 63]. As introduced ascidians can fundamentally alter both the structure and composition of benthic communities [15, 64], reduce the diversity of native species  and threaten economically important species [39, 66], it is imperative to define whether they are native or introduced in any particular region. However, phylogenetic and phylogeographic information may often be insufficient to determine this conclusively, even in conspicuous species such as the members of the P. stolonifera species complex.
Species within Pyura stolonifera
In the taxonomic literature, three species associated with the Pyura stolonifera species complex have traditionally been considered, namely Pyura stolonifera sensu stricto, P. herdmani (Drasche, 1884) and P. praeputialis (Heller, 1987), but the validity of the latter two has been challenged [35, 36] and their species names are not consistently applied . For example, Castilla et al.  found genetic differentiation between the African P. herdmani (probably misidentified as P. stolonifera [see 67]), and populations from eastern Australia and Chile, and recommended referring to the latter as P. praeputialis. However, most subsequent studies continued to refer to the eastern Australian population as P. stolonifera, e.g. [68–77]. We found that P. stolonifera sensu stricto is restricted to temperate southern Africa (Figure 1) and that its range overlaps with that of P. herdmani, which occurs in temperate, subtropical and tropical southern Africa, as well as Morocco (Table 4). Confirming the findings of Castilla et al. , we found that P. praeputialis is both morphologically  and genetically (this study) distinct from its African congeners. The extensive sampling in Australasia revealed the existence of two more species within the species complex, namely P. dalbyi Rius & Teske, 2011  and Pyura sp. (Victoria, South Australia, Tasmania and New Zealand), which has yet to be formally described.
Three of the species identified (Pyura herdmani, P. praeputialis and Pyura sp.) can be further subdivided into regional genetic lineages on the basis of reciprocal monophyly or differences in allele frequencies, and the populations at several localities are likely to be the product of long-distance dispersal. These issues are discussed in the following two sections.
Based on our genetic data, four populations may be the product of long-distance colonisation events, as they are genetically very similar to populations that are isolated from them by large geographic distances. These include the populations in Chile, New Zealand, Morocco and Western Australia. However, these populations differ considerably from each other both in terms of genetic diversity and in terms of how genetically distinct they are from their putative source populations.
Castilla et al.  found that the populations in eastern Australia and Antofagasta, Chile, are genetically very similar, but we considered this to be insufficient to conclude that the Chilean population has recently been introduced, as it could also indicate incomplete lineage sorting among significantly differentiated populations. It was thus considered necessary to obtain larger sample sizes from each population and to compare genetic diversity of the different populations. Lack of genetic structure between P. praeputialis from sites in eastern Australia and Chile supports the idea that the Chilean population is the product of a recent introduction, despite its high genetic diversity. Invasive marine invertebrates can show high genetic diversity due to multiple introductions of large numbers of individuals from different sources, e.g. [78, 79] and the Chilean population includes a random sample of haplotypes found in eastern Australia. A very different result was found for Pyura sp., which has only recently been reported from New Zealand . All 42 individuals from New Zealand had the same ANT allele, indicating loss of diversity through genetic drift or a strong bottleneck effect. The allele found in New Zealand is part of a cluster of haplotypes present in southeastern Australia, and the population in New Zealand clearly does not represent a distinct species. The Western Australian population of P. dalbyi was probably recently founded by individuals from southeastern Australia. Distribution records indicate that this species is absent from the Great Australian Bight  and thus has a disjunct distribution typical of an introduced species, with c 2500 km between its two regional populations. Even more compellingly, while P. dalbyi is common in Victoria (Table 1), it seems to be confined to only two Western Australian sites that are more than 1000 km apart . At one of these (Albany), it has been found exclusively inside the harbour, suggesting that it has failed to spread beyond this point of introduction.
Thus we have at least three populations (in Chile, New Zealand and Western Australia, Table 4) that have apparently been recently introduced through human activities and that should be controlled if possible, even though they may provide biogenic habitat for other species.
Lastly, the Moroccan population of P. herdmani was recovered as a distinct lineage with high nodal support using the most quickly evolving marker used in this study, mtDNA COI. Molecular dating indicated that, like several other marine invertebrates with similar antitropical distributions [80, 81], it diverged from its southern African sister lineage prior to the Holocene.
Cryptic divergence within regions
In several cases, we identified genetic sub-structure within individual taxa (Pyura herdmani, P. praeputialis and Pyura sp.) that may point to the existence of additional cryptic species. While there is no evidence that any of these have become invasive elsewhere, their existence highlights the importance of sampling throughout the entire native range of a taxon suspected of having become invasive. Failure to capture all of the genetic diversity present within a particular region will result in a recently introduced species being mistaken for an indigenous species that was previously overlooked.
The most clear-cut example of cryptic divergence was found in Pyura herdmani between the temperate and subtropical/tropical provinces in southern Africa, which are inhabited by distinct lineages whose ranges overlap on the southeast coast. Phylogeographic disjunctions that coincide with water temperature have been documented in this region for various other marine organisms, and claims that the genetic lineages identified constitute cryptic species could in several cases be supported by morphological and physiological data [82–84]. In Pyura sp., two closely related lineages were identified. The fact that both were present at the same sites in Victoria and Tasmania and that several individuals had ANT alleles from both lineages, indicates that if regional genetic structure ever existed, it has almost completely eroded as a result of subsequent high levels of gene flow. Lastly, in P. praeputialis we found a very recently established genetic disjunction (based on allele frequency differences) across the Bass Strait. Phylogeographic breaks in this region have been documented for a large number of marine species [85, 86] and in most cases, allopatric speciation due to the rise of the former Bassian Isthmus that connected Tasmania with the Australian mainland during periods of low sea-level has been invoked .
Distinguishing a putative invader from a previously overlooked cryptic species can be a challenging task. Our results highlight the importance of extensive sampling to differentiate between native and introduced ranges in widespread marine invertebrates and illustrate the difficulty of correctly identifying non-indigenous species in marine invertebrates with poorly resolved taxonomy.
When attempting to match a population that is suspected of having been recently introduced to a source population, samples meant to represent a particular region usually originate from only a small portion of a species' local range, e.g. [28, 34, 87]. As there may be considerable variation in habitat quality along the range of widely distributed coastal species , such a sampling design can result in incorrect conclusions being drawn on whether populations are exotic or native when multiple genetic lineages are present within regions. In our case, some of the species identified have a preference for sheltered conditions (Pyura herdmani, P. dalbyi and Pyura sp.), whereas others can also be found at exposed sites on the open coast (P. stolonifera and P. praeputialis) . Even more importantly, the fact that co-distributed coastal invertebrates in Australia, South Africa and North America tend to have congruent phylogeographic patterns that are often linked to well-documented marine biogeographic disjunctions, e.g. [84–86, 89, 90], indicates that it is crucial to collect samples in all biogeographic provinces in which a widespread species is represented (e.g. P. herdmani). To achieve good sampling cover, it is thus necessary to collect samples at as many sites as possible rather than obtaining large numbers of samples from a small number of sites. The latter approach is commonly used in population genetic studies in order to accurately estimate genetic diversity at each site, but such information is of little value when the aim of a study is to identify the source population of a putative invader.
Failure to identify and control a non-indigenous species could lead to habitat monopolisation at the expense of native species (e.g. in our study the populations in Chile, New Zealand and Western Australia), while the removal of an organism mistakenly identified as being invasive would constitute habitat destruction and may even result in the extinction of a native species (e.g. Pyura herdmani in Morocco). An inadequate sampling design, in which large numbers of sequences are generated but not all of the evolutionary lineages present in a particular region are recovered, can give researchers a false sense of confidence about the alien or indigenous status of poorly known marine organisms. This may obstruct management efforts aimed at controlling an introduced species during the critical early stages of an invasion.
Grosholz E: Ecological and evolutionary consequences of coastal invasions. Trends in Ecology and Evolution. 2002, 17: 22-27. 10.1016/S0169-5347(01)02358-8.
Geller JB, Darling JA, Carlton JT: Genetic Perspectives on Marine Biological Invasions. Annual Review of Marine Science. 2010, 2: 367-393. 10.1146/annurev.marine.010908.163745.
Coleman FC, Williams SL: Overexploiting marine ecosystem engineers: potential consequences for biodiversity. Trends in Ecology & Evolution. 2002, 17: 40-44. 10.1016/S0169-5347(01)02330-8.
Harley CDG, Hughes AR, Hultgren KM, Miner BG, Sorte CJB, Thornber CS, Rodriguez LF, Tomanek L, Williams SL: The impacts of climate change in coastal marine systems. Ecology Letters. 2006, 9: 228-241. 10.1111/j.1461-0248.2005.00871.x.
Branch GM, Steffani CN: Can we predict the effects of alien species? A case-history of the invasion of South Africa by Mytilus galloprovincialis (Lamarck). Journal of Experimental Marine Biology and Ecology. 2004, 300: 189-215. 10.1016/j.jembe.2003.12.007.
Ross CA, Auge H: Invasive Mahonia plants outgrow their native relatives. Plant Ecology. 2008, 199: 21-31.
Dukes JS, Mooney HA: Disruption of ecosystem processes in western North America invasive species. Revista Chilena de Historia Natural. 2004, 77: 411-437.
Wallentinus I, Nyberg CD: Introduced marine organisms as habitat modifiers. Marine Pollution Bulletin. 2007, 55: 323-332. 10.1016/j.marpolbul.2006.11.010.
Carlton JT, Geller JB: Ecological roulette: the global transport of nonindigenous marine organisms. Science. 1993, 261: 78-82. 10.1126/science.261.5117.78.
Grosholz ED: Recent biological invasion may hasten invasional meltdown by accelerating historical introductions. Proceedings of the National Academy of Sciences of the United States of America. 2005, 102: 1088-1091. 10.1073/pnas.0308547102.
McQuaid CD, Arenas F: Biological invasions: insights from marine benthic communities. Marine Hard Bottom Communities. Edited by: Wahl M. 2009, Ecological Studies 206, Springer-Verlag Berlin, Heidelberg, 309-320.
Carlton JT: Deep Invasion Ecology and the Assembly of Communities in Historical Time. Biological Invasions in Marine Ecosystems. Edited by: Rilov G, Crooks JA. 2009, Berlin Heidelberg: Springer-Verlag, 13-56.
Bickford D, Lohman DJ, Sodhi NS, Ng PKL, Meier R, Winker K, Ingram KK, Das I: Cryptic species as a window on diversity and conservation. Trends in Ecology & Evolution. 2006, 22: 148-155.
Lambert G: Invasive sea squirts: A growing global problem. Journal of Experimental Marine Biology and Ecology. 2007, 342: 3-4. 10.1016/j.jembe.2006.10.009.
Castilla JC, Guiñez R, Caro AU, Ortiz V: Invasion of a rocky intertidal shore by the tunicate Pyura praeputialis in the Bay of Antofagasta, Chile. Proceedings of the National Academy of Sciences of the United States of America. 2004, 101: 8517-8524. 10.1073/pnas.0401921101.
Castilla JC, Guiñez R, Alvarado JL, Pacheco C, Varas M: Distribution, population structure, population biomass and morphological characteristics of the tunicate Pyura stolonifera in the Bay of Antofagasta, Chile. Marine Ecology. 2000, 21: 161-174. 10.1046/j.1439-0485.2000.00695.x.
Branch GM, Griffiths CL, Branch ML, Beckley LE: Two Oceans: A guide to the marine life of southern Africa. 2010, Cape Town: Struik Publishers, Revised edn
Rius M, Branch GM, Griffiths CL, Turon X: Larval settlement behaviour in six gregarious ascidians in relation to adult distribution. Marine Ecology Progress Series. 2010, 418: 151-163.
Svane I, Young CM: The ecology and behaviour of ascidian larvae. Oceanography and Marine Biology Annual Review. 1989, 27: 45-90.
Millar RH: The biology of ascidians. Advances in Marine Biology. 1971, 9: 1-100.
Rius M, Pascual M, Turon X: Phylogeography of the widespread marine invader Microcosmus squamiger (Ascidiacea) reveals high genetic diversity of introduced populations and non-independent colonizations. Diversity and Distributions. 2008, 14: 818-828. 10.1111/j.1472-4642.2008.00485.x.
Zhan A, Macisaac HJ, Cristescu ME: Invasion genetics of the Ciona intestinalis species complex: from regional endemism to global homogeneity. Molecular Ecology. 2010, 19: 4678-4694. 10.1111/j.1365-294X.2010.04837.x.
Locke A: A screening procedure for potential tunicate invaders of Atlantic Canada. Aquatic Invasions. 2009, 4: 71-79. 10.3391/ai.2009.4.1.7.
Carlton JT: Patterns of transoceanic marine biological invasions in the Pacific Ocean. Bulletin of Marine Science. 1987, 41: 452-465.
Lambert G: Ecology and natural history of the protochordates. Canadian Journal of Zoology. 2005, 83: 34-50. 10.1139/z04-156.
Coutts ADM, Dodgshun TJ: The nature and extent of organisms in vessel sea-chests: A protected mechanism for marine bioinvasions. Marine Pollution Bulletin. 2007, 54: 875-886. 10.1016/j.marpolbul.2007.03.011.
Pérez-Portela R, Duran S, Palacín C, Turon X: The genus Pycnoclavella (Ascidiacea) in the Atlanto-Mediterranean region: a combined molecular and morphological approach. Invertebrate Systematics. 2007, 21: 187-205. 10.1071/IS06015.
Caputi L, Andreakis N, Mastrototaro F, Cirino P, Vassillo M, Sordino P: Cryptic speciation in a model invertebrate chordate. Proceedings of the National Academy of Sciences. 2007, 104: 9364-9369. 10.1073/pnas.0610158104.
Iannelli F, Pesole G, Gissi C: Mitogenomics reveals two cryptic species in Ciona intestinalis. Trends in Genetics. 2007, 23: 419-422. 10.1016/j.tig.2007.07.001.
Nydam ML, Harrison RG: Genealogical relationships within and among shallow-water Ciona species (Ascidiacea). Marine Biology. 2007, 151: 1839-1847. 10.1007/s00227-007-0617-0.
Brunetti R: Botryllid species (Tunicata, Ascidiacea) from the Mediterranean coast of Israel, with some considerations on the systematics of Botryllinae. Zootaxa. 2009, 2289: 18-32.
Barros RC, Rocha RM, Pie MR: Human-mediated global dispersion of Styela plicata (Tunicata, Ascidiacea). Aquatic Invasions. 2009, 4: 45-57. 10.3391/ai.2009.4.1.4.
Dalby JEJ: Consequences of aggregated living in the ascidian Pyura stolonifera: Evidence for non-contact intraspecific competition. Marine and Freshwater Research. 1995, 46: 1195-1199. 10.1071/MF9951195.
Castilla JC, Collins AG, Meyer CP, Guiñez R, Lindberg DR: Recent introduction of the dominant tunicate, Pyura praeputialis (Urochordata, Pyuridae) to Antofagasta, Chile. Molecular Ecology. 2002, 11: 1579-1584. 10.1046/j.1365-294X.2002.01524.x.
Kott P: The Australian Ascidiacea, Part 1. Phlebobranchia and Stolidobranchia. Memoirs of the Queensland Museum. 1985, 23: 1-438.
Kott P: Observations on non-didemnid ascidians from Australian waters. Journal of Natural History. 2006, 40: 169-234. 10.1080/00222930600621601.
Monniot C, Bitar G: Sur la présence de Pyura stolonifera (Tunicata, Ascidiacea) à Ras Achaccar (côte nord atlantique marocaine). Comparaison anatomique distinctive avec Pyura praeputialis . Bulletin de l'Institut Scientifique, Rabat. 1983, 7: 83-91.
Lafargue F, Wahl M: Contribution to the knowledge of littoral ascidians (Ascidiacea, Tunicata) of the Senegalese coast. Bulletin de l'IFAN. 1986, 46: 385-402.
Hayward BW, Morley MS: Introduction to New Zealand of two sea squirts (Tunicata, Ascidiacea) and their subsequent dispersal. Records of the Auckland Museum. 2009, 46: 5-14.
Van Name WG: The North and South American ascidians. Bulletin of the American Museum of Natural History. 1945, 84: 1-476.
Fielding PJ, Weerts KA, Forbes AT: Macroinvertebrate communities associated with intertidal and subtidal beds of Pyura stolonifera (Heller) (Tunicata: Ascidiacea) on the Natal coast. South African Journal of Zoology. 1994, 29: 46-53.
Kumar S, Degnan BM, Ross IL, Hawkins CJ, Lavin MF: Isolation of DNA and RNA from ascidians. Marine Biology. 1988, 98: 95-100. 10.1007/BF00392663.
Sunnucks P, Hales D: Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus Sitobion (Hemiptera: Aphididae). Molecular Biology and Evolution. 1996, 13: 510-524.
Jarman SN, Ward RD, Elliott NG: Oligonucleotide primers for PCR amplification of coelomate introns. Marine Biotechnology. 2002, 4: 347-355. 10.1007/s10126-002-0029-6.
Flot J-F: Champuru 1.0: a computer software for unraveling mixtures of two DNA sequences of unequal lengths. Molecular Ecology Notes. 2007, 7: 974-977. 10.1111/j.1471-8286.2007.01857.x.
Stephens M, Smith N, Donnelly P: A new statistical method for haplotype reconstruction from population data. American Journal of Human Genetics. 2001, 68: 978-989. 10.1086/319501.
Tamura K, Dudley J, Nei M, Kumar S: MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0. Molecular Biology and Evolution. 2007, 24: 1596-1599. 10.1093/molbev/msm092.
Notredame C, Higgins DG, Heringa J: T-Coffee: a novel method for fast and accurate multiple sequence alignment. Journal of Molecular Biology. 2000, 302: 205-217. 10.1006/jmbi.2000.4042.
Talavera G, Castresana J: Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments. Systematics Biology. 2007, 56: 564-577. 10.1080/10635150701472164.
Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.
Saitou N, Nei M: The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular Biology and Evolution. 1987, 4: 406-425.
Tamura K, Nei M, Kumar S: Prospects for inferring very large phylogenies by using the neighbor-joining method. Proceedings of the National Academy of Science of the United States of America. 2004, 101: 11030-11035. 10.1073/pnas.0404206101.
Tamura K, Nei M: Estimations of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans in chimpanzees. Molecular Biology and Evolution. 1993, 10: 512-526.
Bandelt H-J, Forster P, Roöhl A: Median-joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution. 1999, 16: 37-48.
Excoffier L, Lischer HEL: Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources. 2010, 10: 564-567. 10.1111/j.1755-0998.2010.02847.x.
Nielsen R, Wakeley J: Distinguishing migration from isolation. A Markov chain Monte Carlo approach. Genetics. 2001, 158: 885-896.
Hey J, Nielsen R: Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proceedings of the National Academy of Science of the United States of America. 2007, 104: 2785-2790. 10.1073/pnas.0611164104.
Hasegawa M, Kishino H, Yano T: Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution. 1985, 22: 160-174. 10.1007/BF02101694.
Posada D, Crandall KA: Selecting the best-fit model of nucleotide substitution. Systematic Biology. 2001, 50: 580-601.
Egan EA: The seasonal reproductive cycle of the nemertean Gononemertes australiensis Gibson in relation to that of its ascidian host, Pyura pachydermatina (Herdman). Journal of Experimental Marine Biology and Ecology. 1984, 76: 225-246. 10.1016/0022-0981(84)90190-4.
Knowlton N, Weigt LA: New dates and new rates for divergence across the Isthmus of Panama. Proceedings of the Royal Society B. 1998, 265: 2257-2263. 10.1098/rspb.1998.0568.
Marko PB: Fossil calibration of molecular clocks and the divergence times of geminate species pairs separated by the Isthmus of Panama. Molecular Biology and Evolution. 2002, 19: 2005-2021.
Monteiro SM, Chapman MG, Underwood AJ: Patches of the ascidian Pyura stolonifera (Heller, 1878): structure of habitat and associated intertidal assemblages. Journal of Experimental Marine Biology and Ecology. 2002, 270: 171-189. 10.1016/S0022-0981(02)00020-5.
Rius M, Pineda MC, Turon X: Population dynamics and life cycle of the introduced ascidian Microcosmus squamiger in the Mediterranean Sea. Biological Invasions. 2009, 11: 2181-2194. 10.1007/s10530-008-9375-2.
Blum JC, Chang AL, Liljesthröm M, Schenk ME, Steinberg MK, Ruiz GM: The non-native solitary ascidian Ciona intestinalis (L.) depresses species richness. Journal of Experimental Marine Biology and Ecology. 2007, 342: 5-14. 10.1016/j.jembe.2006.10.010.
Bullard SG, Lambert G, Carman MR, Byrnes J, Whitlatch RB, Ruiz G, Miller RJ, Harris L, Valentine PC, Collie JS, et al: The colonial ascidian Didemnum sp. A: Current distribution, basic biology and potential threat to marine communities of the northeast and west coasts of North America. Journal of Experimental Marine Biology and Ecology. 2007, 342: 99-108. 10.1016/j.jembe.2006.10.020.
Rius M, Teske PR: A revision of the Pyura stolonifera species complex (Tunicata, Ascidiacea), with a description of a new species from Australia. Zootaxa. 2011, 2754: 27-40.
Raftos DA, Robbins J, Newton R, Nair S: A complement component C3a-like peptide stimulates chemotaxis by hemocytes from an invertebrate chordate - the tunicate, Pyura stolonifera. Comp Biochem Physiol Mol Integr Physiol. 2003, 134: 377-386. 10.1016/S1095-6433(02)00287-8.
Knott NA, Davis AR, Buttemer WA: Passive flow through an unstalked intertidal ascidian: Orientation and morphology enhance suspension feeding in Pyura stolonifera. Biol Bull. 2004, 207: 217-224. 10.2307/1543210.
Barnes PB, Davis AR, Roberts DE: Sampling patchily distributed taxa: a case study using cost-benefit analyses for sponges and ascidians in coastal lakes of New South Wales, Australia. Mar Ecol Prog Ser. 2006, 319: 55-64.
Pinto MR, Melillo D, Giacomelli S, Styroera G, Lambris JD: Ancient origin of the complement system: Emerging invertebrate models. Adv Exp Med Biol. 2007, 598: 372-388. 10.1007/978-0-387-71767-8_26.
Ross DJ, Keough MJ, Longmore AR, Knott NA: Impacts of two introduced suspension feeders in Port Phillip Bay, Australia. Mar Ecol Prog Ser. 2007, 340: 41-53.
Dafforn KA, Glasby TM, Johnston EL: Links between estuarine condition and spatial distributions of marine invaders. Diversity and Distributions. 2009, 15: 807-821. 10.1111/j.1472-4642.2009.00587.x.
Harasti D, Glasby TM, Martin-Smith KM: Striking a balance between retaining populations of protected seahorses and maintaining swimming nets. Aquatic Conservation: Marine and Freshwater Ecosystems. 2009, 20: 159-166.
Beck HJ, Styan CA: Colour patterns in the sea urchin, Heliocidaris erythrogramma, suggest limited connectivity across the Southern and Pacific Ocean coastlines of Australia. Mar Freshwater Res. 2010, 61: 143-152. 10.1071/MF08156.
Wilson GS, Raftos DA, Corrigan SL, Nair SV: Diversity and antimicrobial activities of surface-attached marine bacteria from Sydney Harbour, Australia. Microb Res. 2010, 165: 300-311. 10.1016/j.micres.2009.05.007.
Marshall DJ, Evans JP: Context-dependent genetic benefits of polyandry in a marine hermaphrodite. Biol Lett. 2007, 3: 685-688. 10.1098/rsbl.2007.0438.
Roman J: Diluting the founder effect: cryptic invasions expand a marine invader's range. Proceedings of the Royal Society B. 2006, 273: 2453-2459. 10.1098/rspb.2006.3597.
Simon-Bouhet B, Garcia-Meunier P, Viard F: Multiple introductions promote range expansion of the mollusc Cyclope neritea (Nassariidae) in France: evidence from mitochondrial sequence data. Molecular Ecology. 2006, 15: 1699-1711. 10.1111/j.1365-294X.2006.02881.x.
D'Amato ME, Harkins GW, de Oliveira T, Teske PR, Gibbons MJ: Molecular dating and biogeography of the neritic krill Nyctiphanes. Marine Biology. 2008, 155: 243-247. 10.1007/s00227-008-1005-0.
Palero F, Lopes J, Abello P, Macpherson E, Pascual M, Beaumont MA: Rapid radiation in spiny lobsters (Palinurus spp) as revealed by classic and ABC methods using mtDNA and microsatellite data. BMC Evolutionary Biology. 2009, 9: 263-10.1186/1471-2148-9-263.
Ridgway TM, Stewart BA, Branch GM, Hodgson AN: Morphological and genetic differentiation of Patella granularis (Gastropoda: Patellidae): recognition of two sibling species along the coast of southern Africa. Journal of Zoology. 1998, 245: 317-333. 10.1111/j.1469-7998.1998.tb00107.x.
Teske PR, McLay C, Sandoval-Castillo J, Papadopoulos I, Newman BK, Carvalho D, Griffiths CL, McQuaid CD, Barker NP, Borgonie G, et al: Tri-locus sequence data reject a "Gondwanan origin hypothesis" for the African/South Pacific crab genus Hymenosoma. Molecular Phylogenetics and Evolution. 2009, 53: 23-33. 10.1016/j.ympev.2009.05.031.
Teske PR, Papadopoulos I, Newman BK, Dworschak PC, McQuaid CD, Barker NP: Oceanic dispersal barriers, adaptation and larval retention: an interdisciplinary assessment of potential factors maintaining a phylogeographic break between sister lineages of an African prawn. BMC Evolutionary Biology. 2008, 8: 341-10.1186/1471-2148-8-341.
Waters JM: Marine biogeographic disjunction in temperate Australia: historic landbridge, contemporary currents, or both?. Diversity and Distributions. 2008, 14: 692-700. 10.1111/j.1472-4642.2008.00481.x.
Ayre DJ, Minchinton TE, Perrin C: Does life history predict past and current connectivity for rocky intertidal invertebrates across a marine biogeographic barrier?. Molecular Ecology. 2009, 18: 1887-1903. 10.1111/j.1365-294X.2009.04127.x.
Stefaniak L, Lambert G, Gittenberger A, Zhang H, Lin S, Whitlatch RB: Genetic conspecificity of the worldwide populations of Didemnum vexillum Kott, 2002. Aquatic Invasions. 2009, 4: 29-44. 10.3391/ai.2009.4.1.3.
Sotka EE, Wares JP, Barth JA, Grosberg RK, Palumbi SR: Strong genetic clines and geographical variation in gene flow in the rocky intertidal barnacle Balanus glandula. Mol Ecol. 2004, 13: 2143-2156. 10.1111/j.1365-294X.2004.02225.x.
Teske PR, Winker H, McQuaid CD, Barker NP: A tropical/subtropical biogeographic disjunction in southeastern Africa separates two Evolutionarily Significant Units of an estuarine prawn. Marine Biology. 2009, 165: 1265-1275.
Pelc RA, Warner RR, Gaines SD: Geographical patterns of genetic structure in marine species with contrasting life histories. Journal of Biogeography. 2009, 36: 1881-1890. 10.1111/j.1365-2699.2009.02138.x.
Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R: DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine Biology and Biotechnology. 1994, 3: 294-299.
Weekers PHH, Gast RJ, Fuerst PA, Byers TJ: Sequence variations in small-subunit ribosomal RNAs of Hartmanella vermiformis and their phylogenetic implications. Molecular Biology and Evolution. 1994, 11: 684-690.
Samraoui B, Weekers PHH, Dumont HJ: Two taxa within the North African Lestes virens complex (Odonata, Zygoptera). Odonatologica. 2003, 32: 131-142.
We thank Rocío Pérez-Portela, Syd Ramdhani, Peter Tung, Fabien Forget and Isabelle Papadopoulos for providing additional samples. Sampling in NSW was conducted under NSW DPI research permit 05⁄0090. PR Teske was supported by a postdoctoral research fellowship for overseas study by the National Research Foundation, an overseas study grant from the Ernest Oppenheimer Memorial Trust and a Rhodes University postdoctoral research fellowship. M Rius was supported by the 'Agencia Española de Cooperación Internacional para el Desarrollo' from the Spanish 'Ministerio de Asuntos Exteriores y de Cooperación'. The work was partially funded by a grant from the DST-NRF Centre of Excellence for Invasion Biology (to M Rius), by the South African Research Chairs Initiative of the Department of Science and Technology and the National Research Foundation (to CD McQuaid), by grants from Macquarie University (MQ A006162) and the Australian Research Council (DP110101275) to LB Beheregaray, and by Flinders University. This contribution represents manuscript no. 40 of the Molecular Ecology Group for Marine Research (MEGMAR).
PRT designed the study, collected most of the Australian samples, generated most of the sequence data, did the analyses and prepared the manuscript. MR helped with the design of the study, provided most of the South African samples and helped with the laboratory work. CAS, MPP, SB, CFG, KW, MP, CFM, CRMA, GMC and SCB provided additional samples and generated additional sequence data. CMQ, NPB and LBB provided conceptual guidance and logistical support. All authors contributed to the preparation of the manuscript, and read and approved the final version.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.