Cryptic speciation or global spread? The case of a cosmopolitan marine invertebrate with limited dispersal capabilities.

The existence of globally-distributed species with low dispersal capabilities is a paradox that has been explained as a result of human-mediated transport and by hidden diversity in the form of unrecognized cryptic species. Both factors are not mutually exclusive, but relatively few studies have demonstrated the presence of both. Here we analyse the genetic patterns of the colonial ascidian Diplosoma listerianum, a species nowadays distributed globally. The study of a fragment of a mitochondrial gene in localities worldwide revealed the existence of multiple cryptic species. In addition, we found a complex geographic structure and multiple clades occurred in sympatry. One of the species showed strong population structure irrespective of geographical distances, which is coherent with stochastic dispersal linked to human transport. The present study shows the complexity of discerning the role of cryptic diversity from human-driven range shifts worldwide, as well as disentangling the effects of natural and artificial dispersal.

The application of molecular methods to species delimitation has uncovered an overwhelming amount of unrecognized cryptic diversity. Instances of purportedly widespread species have been shown to correspond to distinct species, which failed to be discriminated by conventional taxonomy based on morphology12. This casts doubt on the status of many cosmopolitan species. On the other hand, genetic studies have shown that true cosmopolitanism is often associated with anthropogenic transport (e.g.3). Man-related vectors (e.g. shipping, agriculture, aquaculture) can boost the dispersal potential of the species, allowing transoceanic and transcontinental transport and range expansions well outside their natural capabilities. This provides a fertile ground for studies using genetic markers456 addressed at unravelling the history of primary and secondary introductions, genetic admixture, and other processes inferred from the present-day distribution of genetic diversity in introduced species (e.g.78). This is particularly interesting in groups with short dispersal abilities (e.g., species with lecithotrophic larvae), as genetic tools have the potential of discerning between natural and artificial dispersal.

There is increasing evidence that complex patterns underlie the genetics of widespread species, with structure found at diverse, nested levels (e.g.9). Highly divergent genetic lineages have been detected in organisms whose distributions are cosmopolitan or nearly so, belonging to diverse groups (e.g. Bugula neritina10, Mytilus galloprovincialis11, Amphipholis squamata12). Ascidians in particular are among the most important introduced species in the sea13, and while populations of some introduced ascidians clearly belong to a single species (e.g.81415), in other cases cryptic speciation has been found. A paradigmatic instance is the case of Ciona intestinalis, a model species studied in hundreds of laboratories, which in fact consists of four cryptic species of which two are invasive1617. Another prominent instance is the colonial ascidian Botryllus schlosseri, which comprises at least three different species, with contrasting invasive capabilities18. How evolutionary histories can result in differences in invasive potential remains a challenging field of study18.

For such widespread species, two possibilities arise: they can be a species complex containing cryptic diversity or a single species that has been widely transported by man. Both possibilities are non-mutually exclusive, and can in fact co-occur (e.g.1718). The use of molecular markers can help inferring the processes that led to the present distribution of the species or group of species considered. However, there is a dearth of studies that describe the combined effects of cryptic diversity and human-driven range shifts in cosmopolitan species. This is surprising, given the ecological and economic impact caused by introduced species, for which an accurate assessment of their hidden diversity and cryptic speciation is necessary for adequate management and mitigation policies5.

In this study we analyse the genetic variability of the widespread colonial ascidian Diplosoma listerianum, a species described from NE Atlantic1920 and nowadays distributed globally21. D. listerianum is commonly found in intertidal and subtidal rocks and brown algae in natural habitats around the British Islands (e.g.2223) but also abundantly found in fouling communities within marinas and harbours all around the world (e. g.2425). D. listerianum ranks among the most widely distributed marine invertebrates, but paradoxically the very short free-swimming period of its larva, that lasts only a few hours before settlement and rapid metamorphosis2026, justifies the prediction of genetic structuring even at local scales2728. In addition, D. listerianum presents some biological particularities such as the ability to retain and select exogenous sperm in the oviduct2930 and colonial fusion23, features that may have relevance during the colonization processes.

The worldwide expansion of this species therefore provides a good natural experiment to test the interplay between genetic structure and the homogenizing effect of stochastic artificial transport. Additionally, as in other cosmopolitan species or species complexes, the list of synonymies of Diplosoma listerianum is a long one, comprising over 30 names given at different geographic regions and times31. Hence the aims of our study were a) to investigate whether D. listerianum is a single species or a species complex, and b) to unravel the distribution and genetic structure of genetically homogeneous clades.

Results

We sequenced a fragment of the mitochondrial gene Cytochrome c Oxidase subunit I (COI) with a total length of 531 bp from 234 colonies collected from 14 different localities (Fig. 1). A total of 216 variable sites (40.7%), and 43 haplotypes were found in all the sequences analysed (Table 1). A total of 34 (79%) were private haplotypes.

Figure 1

Diplosoma listerianum distribution.

Left: Sampling locations of D. listerianum. The grey shadow indicates known distribution of the species, pie charts on the map represent clade frequencies for each locality, and pie size is proportional to sample size. AN: Antofagasta, Chile; AR: Arenys, Spain; BA: Bastimentos Island, Panama; BO: Bocas del Toro, Panama; CA: Bodega Bay, California; CO: Coquimbo, Chile; CT: Cape Town, South Africa; HB: Hout Bay, South Africa; JA: Misaki, Japan; PA: Port Alfred, South Africa; PL: Plymouth, UK; SA: Santander Bay, Spain; MB: Melbourne Bay, Australia; WAS: Snog Harbour, Washington. Right: Haplotype frequencies of clade A for the 11 populations analysed. White haplotypes represent private haplotypes and black ones represent haplotypes shared with populations which have not been considered for population genetics analyses. Dashed arrows suggest the most likely way of spreading. This map has been created by R.P.-P. in Adobe Illustrator CS3 Software.

Table 1

Phylogenetic clade, sampling locality, locality code, coordinates, number of individuals analysed (N), number of haplotypes (Nh) including number of private haplotypes in brackets, haplotype codes, haplotype diversity (Hd) and nucleotide diversity (π) with standard deviation. For populations of clade A with N ≥ 10, haplotype richness after rarefaction (Hr[10]), Tajima's D and Fu's FS statistics are given. These are also the populations considered for further population genetics analyses

Clade	Locality	Code	Coordinates	N	Nh	Haplotypes	Hd	π	Hr[10]	Tajima's D	Fu's FS
A	Bodega Bay, California. USA	CA	38°19′44.78″N, 123°03′23.99″W	23	4 (2)	H_1, H_2, H_3, H_4	0.486 ± 0.105	0.0313 ± 0.0313	1.831	1.196	16. 906
	Snog Harbour, Washington. USA	WS	48°34′14.91″N, 123°10′02.07″W	19	2	H_5, H_6	0.351 ± 0.111	0.0013 ± 0.0004	0.967	0.541	2.033
	Queen Anne's Battery, Plymouth. UK	PL	50°22′01.00″N, 04°07′52.00″W	10	5 (2)	H_2, H_4, H_7, H_8, H_9	0.667 ± 0.163	0.0323 ± 0.0139	4.000	−0.046	5. 703
	Cape Town, South Africa	CT	33°55′09.77″S, 18°26′35.06″E	15	3	H_2, H_4, H_9	0.562 ± 0.095	0.0406 ± 0.0077	1.666	2. 249	17. 579
	Hout Bay, South Africa	HB	34°02′59.59″S, 18°20′56.04″E	10	4	H_2, H_4, H_9	0.622 ± 0.138	0.0458 ± 0.0081	2.000	2. 427	13. 331
	Port Alfred, South Africa	PA	33°35′36.08″S, 26°53′33.00″E	8	5 (4)	H_14, H_15, H_16, H_17, H_18	0.857 ± 0.108	0.0601 ± 0.0151	-	-	-
	Bocas marina, Panama	BO	9°20′08.16″N, 82°13′08.76″W	18	9 (7)	H_2, H_18, H_19, H_20, H_21, H_22, H_23, H_24, H_25	0.824 ± 0.075	0.0218 ± 0.0099	4.866	−1.798*	2.404
	Bastimentos Island, Panama	BA	9°29′13.28″N, 82°10′35.84″W	3	2	H_2, H_18	0.667 ± 0.314	0.0703 ± 0.0331	-	-	-
	Santander Bay, Spain	SA	43°27′43.68″N, 03°47′40.59″W	17	4 (2)	H_2, H_4, H_28, H_29	0.625 ± 0.108	0.0306 ± 0.0097	2.485	0.509	13.178
	Arenys, Spain	AR	41°34′38.83″N, 02°33′24.18″E	23	4 (2)	H_2, H_4, H_42, H_43	0.320 ± 0.121	0.0208 ± 0.0094	1.561	−1.515	12. 460
	Coquimbo, Chile	CO	29°57′17.86″S, 71°20′05.49″W	21	11 (7)	H_2, H_4, H_5, H_6, H_34, H_35. H_37, H_38, H_39, H_40, H_41	0.824 ± 0.084	0.0397 ± 0.0085	5.500	0.669	2.954
	Antofagasta, Chile	AN	23°38′34.61″S, 70°23′58.94″W	11	3	H_4, H_5, H_34	0.564 ± 0.134	0.0162 ± 0.0121	1.000	−2.125	7. 887
	Melbourne Bay, Australia	MB	38°08′46.60″S, 144°56′56.58″E	12	2	H_2, H_4	0.530 ± 0.076	0.0459 ± 0.0066	1.000	2. 753	19. 994
	Kanagawa, Japan	JA	35°08′33.67″N, 139°37′22.26″E	2	1	H_33	0.000	0.000	-	-	-
	TOTAL Clade A			192	34 (26)		0.826 ± 0.019	0.0582 ± 0.0022		1.097	17.710
B	Kanagawa, Japan	JA	35°08′33.67″N, 139°37′22.26″E	10	3 (3)	H_30, H_31, H_32
	TOTAL Clade B			10	3 (3)		0.378 ± 0.181	0.0169 ± 0.0101
C	Bastimentos Island, Panama	BA	9°29′13.28″N, 82°10′35.84″W	2	1	H_11	0.000	0.000
	Hout Bay, South Africa	HB	34°02′59.59″S, 18°20′56.04″E	11	1 (1)	H_10	0.000	0.000
	Port Alfred, South Africa	PA	33°35′36.08″S, 26°53′33.00″E	13	3 (2)	H_11, H_12,H_13	0.295 ± 0.156	0.0023 ± 0.0017
	TOTAL Clade C			26	4 (3)		0.591 ± 0.053	0.0342 ± 0.0026
D	Bastimentos Island, Panama	BA	9°29′13.28″N, 82°10′35.84″W	6	2 (2)	H_26, H_27	0.533 ± 0.172	0.0020 ± 0.0007
TOTAL				234	43 (34)		0.877 ± 0.014	0.0984 ± 0.0049

*Significant at P < 0.05.

Phylogenetic trees reconstructed based on Bayesian Inference (BI), Maximum Likelihood (ML), and Maximum Parsimony (MP) criteria grouped the COI haplotypes of Diplosoma listerianum into four well supported monophyletic clades (posterior probabilities of 1.00, bootstrap values ≥90%) henceforth named clades A, B, C and D (Fig. 2). Interclade genetic divergence ranged from a 17% between clades A and D to a 20% between clades C and B, and C and A (see Table 2), and intraclade variability was always comparatively lower (0.2%–7.3%). However, the phylogenetic relationships among clades could not be completely resolved with this gene fragment (see Fig. 2), as reflected by overall low support and some unstable between-clade relationships depending on the reconstruction method. Likewise, relationships within the major clade (A) were not clearly defined and varied according to the method used, although a subclade A.1, grouping 8 haplotypes, appeared in all cases with strong support (posterior probability of 1, bootstrap values of 100%, Fig. 2).

Figure 2

Phylogenetic tree.

BI consensus tree of haplotypes of D. listerianum. Four main clades, (A) (and subclade A.1), (B), (C) and (D) are highlighted. For inter-calde relationships and for the internal arrangement within clade A, as there were differences among the three methods, the BI topology (and associated posterior probabilities) is shown. Branches retrieved by the three methods are indicated by three support values on the nodes. Values represent posterior probabilities for BI when >0.5, and bootstrap supports when >50% for ML and MP analyses, in that order. A sequence of Diplosoma spongiforme (Acc. number AY600972.1) was included as an outgroup.

Table 2

Percentage of genetic divergence (based on p-distances) between cryptic clades of Diplosoma listerianum and D. spongiforme for the COI gene. Intraclade variability is also shown

	Clade A	Clade B	Clade C	Clade D
Clade B	17.3%
Clade C	20.0%	20.0%
Clade D	17.0%	19.3%	17.6%
D. spongiforme	21.5%	20.0%	21.7%	20.0%
Intraclade	7.3%	1.7%	0.2%	0.2%

An unrooted network constructed using the complete dataset of sequences also supported the existence of four markedly divergent clades lacking intermediate haplotypes (Fig. 3). Clade A, which appeared in all localities analysed, was the most diverse and widely distributed, including 34 haplotypes and a number of “missing” haplotypes that had to be inferred to fully connect the network. In addition, clade A was the most frequent in most localities, with the exception of Bastimentos Island (Panama), Hout Bay and Port Alfred (both from South Africa), and Kanagawa, the only locality from Japan. Within this clade, two distant haplotypes, haplotypes H_2 and H_4, were the most frequent and were present in most localities. Clade C grouped four haplotypes from only three sites of Panama and South Africa (Bastimentos Island, Hout Bay and Port Alfred) (see Table 1 and Fig. 3), and was the most frequent clade encountered in the South African localities with a 54% of the specimens belonging to this clade. This is despite the fact that H_10 found in Hout Bay was separated by more than 30 mutation steps from other haplotypes of the same clade. Clades B and D were the most geographically restricted clades. Clade B, which grouped three genetically distant haplotypes (H_30, H_31 and H_32), appeared only in Japan but it was the most important genetic clade in that particular locality (83.3% of individuals, see Table 1 and Fig. 1). Clade D, which grouped two closely related haplotypes, only appeared in Bastimentos (Panama) but it was the most frequent clade at this particular locality (54.5 % of individuals, see Table 1 and Fig. 1). Overall, the network did not reveal clear geographic structuring. Only a small group of haplotypes (H_18, H_19, H_20, H_21, H_22, H_23, H_24 and H_25), mostly from Panama, formed a cluster (subclade A.1) within clade A, separated by more than 40 mutation steps from the nearest haplotype. This same cluster was also observed in the phylogenetic trees (Fig. 2).

Figure 3

Haplotype network for Diplosoma listerianum from COI data.

Areas of the circles are proportional to the number of sampled individuals. Partitions inside the circles represent the proportion of each population within each haplotype. Small white crossed dots without name represent missing, probably unsampled haplotypes or extinct sequences. Lines between circles represent one mutational step, and roman numerals are the number of mutations between haplotypes when more than one.

Population genetics

Population genetics analyses of clade A were conducted in 11 populations where we obtained ≥10 individuals of this clade. Details about populations and number of samples considered for these analyses are presented in Table 1. These localities presented different population sizes, and haplotype richness after rarefaction showed that Coquimbo (Chile), Bocas de Toro (Panama), and Plymouth (UK) were in that order the most diverse populations, with a high percentage of private haplotypes (64%, 78%, and 40% respectively) as well. In addition, for the Bocas population, eight of the nine haplotypes (89%) belonged to subclade A.1. On the other hand, Washington (EEUU), Antofagasta (Chile) and Melbourne (Australia) had the lowest values of genetic diversity without presence of private haplotypes in any of them. Some populations, such as Hout Bay, Cape Town and also Melbourne presented low values of haplotype diversity and richness but relatively high values of nucleotide diversity indicating the presence of highly divergent haplotypes within these populations. The neutrality tests computed for the 11 populations did not detect clear evidences of recent demographic events (Table 1).

AMOVA results based on grouping populations of clade A into four different marine basins and two hemispheres did not reveal significant differences among groups (P = 0.654 and P = 0.331, respectively, Table 3). However, the AMOVA showed significant differences among populations either when pooled into groups (P < 0.001) or without grouping (P < 0.001). In all cases most of the genetic variability (69–68%) was observed within populations (Table 3). Further analyses based on pairwise comparisons between populations using the FST and D estimators revealed significant differences in genetic structure between most of the populations (Table 4), with some exceptions: the populations of Hout Bay and Melbourne did not display significant differences with most of the other populations, and Antofagasta (Chile) and Plymouth (UK) did not show significant differences in genetic structure with four populations each (Table 4).

Table 3

AMOVA grouping populations of Diplosoma listerianum (Clade A) according to basins (Atlantic, Mediterranean, Indian and Pacific), hemispheres, and among populations without grouping

SOURCE OF VARIATION	DF	SUM OF SQUARES	VARIANCE COMPONENTS	PERCENTAGE OF VARIATION
4 Basins
Among basins	3	6.674	−0.012	−2.96 (FCT = −0.030, p = 0.654)
Among populations within basins	7	17.335	0.140	34.00 (FSC = 0.330, p < 0.001)
Within populations	168	47.734	0.284	68.96 (FST = 0.310, p < 0.001)
2 Hemispheres
Among hemispheres	1	2.979	0.007	1.60 (FCT = 0.016, p = 0.331)
Among populations within hemispheres	9	21.030	0.128	30.51 (FSC = 0.310, p < 0.001)
Within populations	168	47.734	0.284	67.89 (FST = 0.321, p < 0.001)
Total
Among populations without grouping	10	24.009	0.131	31.59 (FST = 0.316, p < 0.001)
Total	178	71.743	0.419

Table 4

Genetic differentiation (FST and D) between populations of Diplosoma listerianum (Clade A) for the COI gene. FST values are represented below the diagonal, and differentiation D values and CIs (bounded between 0 and 1) are shown above the diagonal. P-values for significance of the FST and for the CI intervals of D were set at 0.0109 following FDR correction. Significant values of D (CI not enclosing 0) and FST are indicated by an asterisk

	CA	WS	PL	CT	HB	BO	SA	AN	CO	MB	AR
CA	-	1.000*	−0.037	0.750*	0.395	0.888*	−0.007	0.742*	0.634*	0.153	0.598*
		(1.000–1.000)	(0.000–0.175)	(0.487–1.000)	(0.000–0.813)	(0.678–1.000)	(0.000–0.147)	(0.494–0.990)	(0.338–0.930)	(0.000–0.504)	(0.311–0.884)
WS	0.577*	-	1.000*	1.000*	1.000*	1.000*	1.000*	1.000*	0.885*	1.000*	1.000*
			(1.000–1.000)	(1.000–1.000)	(1.000–1.000)	(1.000–1.000)	(1.000–1.000)	(1.000–1.000)	(0.708–1.000)	(1.000–1.000)	(1.000–1.000)
PL	−0.019	0.522*	-	0.654*	0.437	0.869*	−0.046	0.855*	0.720*	0.216	0.744*
				(0.274–1.000)	(0.000–0.941)	(0.644–1.000)	(0.000–0.183)	(0.574–1.000)	(0.344–1.000)	(0.000–0.646)	(0.398–1.000)
CT	0.410*	0.552*	0.297*	-	0.183	0.988*	0.759*	0.572*	0.480*	0.511*	0.497*
					(0.000–0.629)	(0.949–1.000)	(0.483–1.000)	(0.223–0.921)	(0.154–0.895)	(0.140–0.882)	(0.127–0.877)
HB	0.255*	0.541*	0.194*	0.114	-	0.960*	0.405	0.174	0.037	−0.023	0.030
						(0.856–1.000)	(0.000–0.835)	(0.000–0.522)	(0.000–0.313)	(0.000–0.267)	(0.000–0.293)
BO	0.328*	0.416*	0.221*	0.300*	0.257	-	0.881*	1.000*	0.985*	0.928*	0.989*
							(0.668–1.000)	(1.000–1.000)	(0.934–1.000)	(0.777–1.00)	(0.956–1.000)
SA	0.004*	0.517*	−0.025	0.341*	0.197	0.250*	-	0.756*	0.624*	0.176	0.627*
								(0.484–1.000)	(0.287–0.961)	(0.000–0.549)	(0.305–0.949)
AN	0.399*	0.527*	0.344*	0.286*	0.041	0.293*	0.327*	-	0.103	0.276	0.162
									(0.000–0.381)	(0.000–0.634)	(0.000–0.476)
CO	0.253*	0.378*	0.188*	0.181*	0.006	0.174*	0.188*	−0.004	-	0.164	0.163
										(0.000–0.453)	(0.000–0.428)
MB	0.132	0.573*	0.138	0.297*	−0.016	0.295*	0.111	0.130	0.062	-	0.099
											(0.000–0.384)
AR	0.469*	0.666*	0.476*	0.404*	0.057	0.441*	0.426*	0.088	0.112*	0.135	-

The MDS based on the FST values did not reveal any clear grouping among populations. Only the populations of Plymouth and Santander, which displayed the lowest value for the pairwise FST estimator, overlapped on the graphical representation. Washington, which clearly differed from all the other populations in haplotype composition with only two haplotypes, both shared with Coquimbo and one with Antofagasta, appeared separated from other populations (Fig. 4). The pattern of genetic differentiation observed between populations was unrelated to geographic origin and no signal of isolation by distance was detected (Mantel test: r = 0.015, P = 0.405).

Figure 4

MDS graph.

Multidimensional scaling plot based on FST values between populations of Diplosoma listerianum.

Discussion

The present study shows a complex scenario that comprises, on the one hand, cryptic diversity and, on the other hand, population genetic signatures that suggest human-mediated dispersal of the species Diplosoma listerianum. Our results show the utility of genetic tools to discern cryptic speciation from human-driven range shifts worldwide, as well as to disentangle the effects of natural and artificial dispersal.

The analysis of the sequences of COI recovered four monophyletic clades in D. listerianum. The large genetic divergences between them (>17%) indicate that D. listerianum is a species complex of, at least, four evolutionarily distinct lineages. We were unable to substantiate any clear diagnostic morphological character due to the scarcity of candidate characters (f.i., lack of calcareous spicules in the tunic of this genus), and the high variability of shape and colour of the colonies. In species of Diplosoma the morphology of the larvae and incubation place can assist species identification32; however, no differences in these characters were found in the brooding colonies examined.

Although COI is the gene of choice in studies of barcoding of animal species, a genetic divergence threshold to separate species has not been defined for ascidians, and the cut-off value varies from group to group. For this reason, we compared the levels of genetic divergence obtained between genetic clades of D. listerianum with those obtained in recent studies of the genus Diplosoma3233. In both studies the COI gene showed a large range of interspecific divergence, from 8.7% to 22%. Values between 17% and 20% found in D. listerianum are then in the upper range of interspecific divergence assessed for this genus. Unfortunately, the fragment sequenced by those authors did not overlap with the gene fragment sequenced in this study, so the extrapolation of their results to our data is not straightforward. Comparing levels of divergence obtained in this study with other ascidians for the same gene fragment demonstrates that divergence levels between clades of D. listerianum are comparable with interspecific ranges of other colonial tunicates, such as the genera Clavelina (15–20%)34, Pycnoclavella (10–21%)34, Botryllus (10–16.5%)18, or the solitary Ciona intestinalis (11.1–18.4%)1735. In addition, our divergence levels exceed by far those reported between sibling ascidian species of Pseudodistoma crucigaster (2.12%)36, Clavelina lepadiformis (5%)37, Pycnoclavella communis (8.55%)28, and Pyura praeputialis (10%)38. Overall, our results indicate that the clades detected can be identified as sibling species that have gone unnoticed. The lack of diagnostic characters renders them as morphologically cryptic species.

There are no calibrated mutation rate estimates for COI in ascidians (due to lack of adequate fossil records), but data from other marine invertebrates is in the range of 1.6 to 2.6% per million years39. This places the split between clades in D. listerianum at ca. 6–12 MYA. A rapid speciation event could explain the lack of resolution of inter-clade relationships. Initial differentiation could have been boosted by geographic discontinuities, and natural range shifts could have occurred afterwards. However, the wide present-day distribution of some of the clades is, most likely, the result of recent translocation by man as the species is mostly found in harbours and artificial structures. The finding of private haplotypes is attributable to sampling effects, not to in situ evolution of these haplotypes. Likewise, the genetic diversity observed in some localities has likely originated in the native area and its current distribution is the result of multiple introductions with associated bottlenecks and, possibly, selection.

The results obtained in this study show important differences in geographical distribution of the different clades within D. listerianum (sensu lato). Differences as such have been related with different levels of invasiveness among clades1718. The clade A of D. listerianum is the most abundant, the one that shows a wider distribution (both sides of the Atlantic, Mediterranean, both sides of the Pacific, and Indian Ocean), and is found in all populations surveyed. This clade is the only one present in European waters and holds the highest genetic diversity. Historical reports of D. listerianum situate its centre of distribution around the British Islands20 where it was first described from the English Channel19. Considering this and our genetic findings, we assign clade A as the original nominal species D. listerianum.

In some cases the geographical distribution of the cosmopolitan species may assist understanding their original native range, but in many others the information is blurred by a long history of introductions that prevents to assign unambiguously a native or introduced status in a given area144041. This may be the case of clade C, the most important species in the easterner populations of South Africa (HB at the Atlantic side and PA at the Indian Ocean), only found outside this area by two specimens in Bastimentos (Panama). The current distribution and levels of genetic diversity observed may be coherent with clade C being native from South Africa and introduced in Panama, although the opposite may be true, or both areas can have been populated from other unsampled sources. Unfortunately, with the data here presented we cannot reach further conclusions about the geographical origin of this clade. The first South African record of D. listerianum is by Millar in 194942, and it has been detected in subsequent surveys2543, but the lack of taxonomic distinction among species in the complex prevents a meaningful historical interpretation. In contrast, the other two cryptic species, clade B and clade D, were restricted to their presumably native ranges. Specimens of clade B were only found in Kanagawa at the Sagami Bay in Japan, where they appeared in sympatry with one specimen of clade A. The less frequent species in our samples, clade C, was the most common one in Bastimentos, but was not found in the nearby locality of Bocas del Toro. It coexists with clade A and clade D in Bastimentos. The large differences in species composition between Bastimentos and Bocas de Toro and the mosaic of species found in the former locality may result from introduction events from multiple sources combined with the presence of the native species.

One relevant question is whether it is possible or advisable to rescue some of the synonymies31 to refer to the sibling species detected. For instance, the Japanese species Diplosoma mitsukurii Oka, 1892, which is undistinguishable from D. listerianum44, could actually be a valid species corresponding to our clade B. Likewise, Diplosoma macdonaldi Herdman, 1886, described from Brazil and common in tropical and subtropical W Atlantic waters45 is another synonym46 that could conceivably correspond to our clade D. In our view, however, it is not advisable to resuscitate these old names because several clades occur in these areas, and because existing morphological descriptions are useless to establish correspondences with our genetic lineages. Given that there is no way to disentangle literature references to D. listerianum and synonyms with respect to the species here reported, it seems adequate for now to note their existence, and to await further morphological studies informed by genetic data before making any formal description of new species.

The analyses of the 11 populations of D. listerianum (clade A) revealed the existence of a complex pattern of population structure influenced by both its limited dispersal potential and human transport. In general, populations of D. listerianum displayed significant differences in allele frequencies that were not related to geographical distances. Although the scenario of an increased human-mediated transport due to an ever increasing shipping activity could counteract genetic differentiation, other factors related to the introduction processes (e.g., bottlenecks), genetic drift and strong selection may be involved in rapid divergence of the populations even if they are derived from the same sources515. Furthermore, we have also observed several evidences of long-distance transport, such as the sharing of haplotypes between distant populations (e.g. Coquimbo and Washington, Melbourne, South Africa and Europe), and genetic homogeneity between very distant geographical areas such as Plymouth, California and Melbourne, which cannot be explained by natural dispersal. The presence of highly differentiated haplotypes within some populations also suggests that multiple introductions from several sources have occurred relatively frequently. Therefore, human-mediated gene flow is likely playing an important role in populations' connectivity and spreading. Since the likelihood of larvae surviving in ballast water is very low due to their short free-swimming period, the most plausible vectors in this species are the transfer of adult colonies in ship hulls and aquaculture41340. As a telling instance, a ferry service connects twice per week the ports of Plymouth and Santander (http://www.brittany-ferries.co.uk/), two localities showing no genetic differentiation (Fig. 4).

The worldwide introduction of D. listerianum as defined here and its success in the new environments likely reflects a high adaptability of the species to new environments47. Biological traits such as the ability of the species to retain and select exogenous sperm in the oviduct, probably from genetically distant partners, during weeks2930 might confer advantages during the introduction and colonization process. Cross-fertilized zygotes can be produced up to a month after the sperm reaches the oviduct, and the brooded larvae are released only when they are competent to settle48, which reduces larval mortality at the expense of dispersal potential. If colonies introduced in a new habitat already contain sperm from distant genotypes and/or embryos they may be able to release cross-fertilized larvae in very short periods of time decreasing the risk of inbreeding depression and allowing a rapid colonization of the new habitat. Additionally, colonies of D. listerianum frequently fuse resulting in chimeras with up to six different genotypes23, therefore a single colony of D. listerianum may contain an important fraction of the genetic diversity of the source populations, which may have important consequences on the long-term survival of the newly established populations15.

We acknowledge that using a single mitochondrial marker could provide an incomplete picture, and studies using nuclear markers are sought in order to better clarify the evolutionary relationships found in this species complex. Nonetheless, the COI gene has repeatedly been shown to be extremely useful and accurate in detecting cryptic speciation in ascidians (e.g.92834363738). Phylogenetic analyses and network structure support that D. listerianum is actually a species complex consisting of at least four highly divergent species with different invasive potential. Two of these species, the putative original D. listerianum (clade A) and clade C, seem to have been introduced in different biogeographical areas, whereas clade B and clade D remain restricted to their putative native ranges. Population genetic analyses of D. listerianum revealed the existence of a complex pattern of population genetics probably influenced by a limited dispersal potential of the larva, large-distance dispersal by human transport and population divergence caused by demographic events during the colonization episodes. The invasive success of D. listerianum may reside on its ability to form colonial chimeras, exogenous storage of sperm, and brooding. Further studies focusing on relating these characteristics with the genetic groups detected could provide insights to understand the differential invasive potential in this group of species. The results presented in this study also point out the need of a complete re-description and evaluation of the D. listerianum species complex. D. listerianum is a remarkable example of how the introduction of a species worldwide and the absence of proper diagnostic characters to distinguish closely related species have hindered the recognition of the presence of local species. Our study highlights how the application of molecular tools to purportedly cosmopolitan taxa is necessary to disentangle the effects of human introductions from those of cryptic speciation.

Methods

Sampling collection

Colonies morphologically attributable to Diplosoma listerianum were collected on artificial substrates (aquaculture settings, harbours and marinas) of 14 different localities covering most of the global species distribution (see details in Table 1 and Fig. 1). Colonies of Diplosoma were collected from ropes and biota attached to the floats, at least 1 m apart to each other to avoid sampling clonal fragments of the same colony. Colonies were directly preserved in absolute ethanol, and stored at −20°C. Once in the laboratory, preserved colonies were dissected under the binocular microscope and several zooids were separated with forceps for DNA extraction.

DNA extraction and sequencing

Total DNA was extracted using a REDExtract-N-Amp kit (Sigma-Aldrich) from one zooid per colony to circumvent the potential presence of chimeras. The tunicate-specific primers Tun_forward, 5′ TCGACTAATCATAAAGATATTAG 3′, and Tun_reverse2, 5′ AACTTGTATTTAAATTACGATC 3′49 were used for the amplification of a fragment of the COI mitochondrial gene. PCR amplification reactions were performed in a total volume of 25 μl with 0.5 μl of each primer (10 μM), 0.5 μl dNTPs (10 μM), 1.25 μl MgCl2 (50 mM), 0.25 μl HotMaster taq polymerase (Invitrogen, www.invitrogen.com), 2.5 μl of 10× buffer, 18.75 μl ultrapure water, and 2 μl of template DNA. A single denaturation step at 94°C for 2 min was followed by 48 cycles (denaturation at 94°C for 1 min, annealing at 38–40°C for 1.5 min and extension at 72°C for 1 min) and a final extension at 72°C for 7 min in a PCR System 9700 thermal cycler (Applied Biosystems). The same primers (Tun_forward and Tun_reverse2) were used for the sequencing reaction. The forward strand was initially sequenced, and whenever any ambiguous base call or poor quality segments were found, the reverse strand was sequenced to correct any problem. The PCR products were sequenced with an ABI Big-Dye Ready-Reaction Perkin Elmer kit on an ABI Prism 377XL automated sequencer (Applied Biosystems) in the Scientific and Technical Services of the University of Barcelona (Barcelona, Spain).

All sequences were edited and aligned using BioEdit Sequence Alignment Editor v. 7.0.8.050 and the results from the alignment verified by eye. Sequences of the haplotypes found in this study have been deposited in GenBank (accession numbers KF791867-KF791909).

Phylogenetic and genetic structure analyses

The complete dataset of sequences was collapsed in haplotypes for phylogenetic analyses. We also included one sequence of Diplosoma spongiforme as an outgroup (Genbank Acc. Number AY600972.1). For phylogenetic reconstruction we used three different criteria: Bayesian Inference (BI), Maximum Likelihood (ML) and Maximum Parsimony (MP). The best-fit model of nucleotide substitution for the data was selected by statistical comparison of 56 different models of evolution with Modeltest v. 3.051 using the Akaike Information Criterion (AIC), and values of the selected evolutionary model were implemented for BI and ML. The BI tree was estimated after 2 million generations with a sample frequency of 100 (20,000 final trees) with the software MrBayes v. 3.1.252. After verifying that stationarity had been reached, the first 2,000 trees were discarded and an independent majority-rule consensus tree was generated from the remaining trees (18,000 trees). The ML analysis was run in PhyML53, and the MP analysis was done as implemented in Seaview. For the last two analyses 1,000 bootstrap replicates were run to assess the robustness of the nodes. Genetic divergences between phylogenetic clades obtained from the reconstructed tree were calculated in Mega v. 5.154 based on the nucleotide p-distances to allow for comparisons with published data on intra- and interspecies distances in related taxa.

The complete data sets of sequences obtained for the COI was used to construct an unrooted network, under the null hypothesis of no genetic differentiation among localities. We used the Network vs. 4.6.1.1 program, which assumes the median-joining network method in the absence of recombination55. This method begins by combining the minimum spanning trees within a single network. With a parsimony criterion, median vectors (which represent missing intermediate haplotypes) are added to the network. Three loops observed in the networks were solved using criteria derived from coalescent theory.

Number of haplotypes (Nh), number of private haplotypes, haplotype diversity (Hd) and nucleotide diversity (π) values were computed with DnaSP v. 5.1056. Further population genetics analyses were restricted to clade A, as we could not mix different genetic pools in these analyses, and clade A was the only one for which enough number of individuals and populations were found. We further selected for analysis populations of clade A for which we had ≥10 colonies (11 populations, Table 1). For these, a corrected haplotype richness (Hr) was calculated after rarefaction (adjusted to the minimum population sample size) with the software CONTRIB, and demographic tests such as Fu's FS and Tajima's D were obtained with DnaSP v. 5.10. Analyses of the molecular variance (AMOVA) using haplotype frequencies were performed to examine the population structure, and their significance was tested running 16 000 permutations in Arlequin v. 3.557. Populations were grouped according to two different criteria; within four marine basins (Atlantic, Mediterranean, Indian and Pacific) and within North and South Hemispheres. In both grouping schemes, we had enough power in the permutation tests to reject the null hypotheses58. An AMOVA analysis was also run without grouping populations. For further analyses of population genetic structure, genetic distances (FST) between populations based on haplotype frequencies were assessed. The significance of the values was evaluated by performing 16,000 permutations with the Arlequin software. The measure of differentiation D proposed by Jost59 was also obtained using SPADE (available at http://chao.stat.nthu.edu.tw), with 1,000 bootstrap replicates used to estimate confidence intervals. A false discovery rate (FDR) correction was applied to the P-values (B-Y method as described by Narum60) to account for multiple tests. A multidimensional scaling analysis (MDS) was performed to graphically visualise interrelationships represented by the matrix of genetic distances derived from the FST values.

Finally, the potential effect of isolation by geographical distance was tested in Arlequin using the Mantel test procedure with 10,000 permutations. For this, we compared the correlation between the genetic distances matrix (FST/(1 − FST)) and a matrix of log-transformed geographical linear distances by sea.

Author Contributions

R.P.-P., M.R. and X.T. designed this study and collected most of the samples here analysed. V.A. sequenced all the samples and analysed part of the data. R.P.-P. and V.A. performed final analyses and prepared figures and tables. R.P.-P. and X.T. wrote the main manuscript and M.R. and V.A. revised it. All authors contributed with their ideas and reviewed the final version of the manuscript. All authors contributed to the revision of the manuscript after peer-review.