Open Access

Deep genetic divergences among Indo-Pacific populations of the coral reef sponge Leucetta chagosensis (Leucettidae): Founder effects, vicariance, or both?

BMC Evolutionary Biology20088:24

DOI: 10.1186/1471-2148-8-24

Received: 07 February 2007

Accepted: 26 January 2008

Published: 26 January 2008



An increasing number of studies demonstrate that genetic differentiation and speciation in the sea occur over much smaller spatial scales than previously appreciated given the wide distribution range of many morphologically defined coral reef invertebrate species and the presumed dispersal-enhancing qualities of ocean currents. However, knowledge about the processes that lead to population divergence and speciation is often lacking despite being essential for the understanding, conservation, and management of marine biodiversity. Sponges, a highly diverse, ecologically and economically important reef-invertebrate taxon, exhibit spatial trends in the Indo-West Pacific that are not universally reflected in other marine phyla. So far, however, processes generating those unexpected patterns are not understood.


We unraveled the phylogeographic structure of the widespread Indo-Pacific coral reef sponge Leucetta chagosensis across its known geographic range using two nuclear markers: the rDNA internal transcribed spacers (ITS 1&2) and a fragment of the 28S gene, as well as the second intron of the ATP synthetase beta subunit-gene (ATPSb-iII). This enabled the detection of several deeply divergent clades congruent over both loci, one containing specimens from the Indian Ocean (Red Sea and Maldives), another one from the Philippines, and two other large and substructured NW Pacific and SW Pacific clades with an area of overlap in the Great Barrier Reef/Coral Sea. Reciprocally monophyletic populations were observed from the Philippines, Red Sea, Maldives, Japan, Samoa, and Polynesia, demonstrating long-standing isolation. Populations along the South Equatorial Current in the south-western Pacific showed isolation-by-distance effects. Overall, the results pointed towards stepping-stone dispersal with some putative long-distance exchange, consistent with expectations from low dispersal capabilities.


We argue that both founder and vicariance events during the late Pliocene and Pleistocene were responsible to varying degrees for generating the deep phylogeographic structure. This structure was perpetuated largely as a result of the life history of L. chagosensis, resulting in high levels of regional isolation. Reciprocally monophyletic populations constitute putative sibling (cryptic) species, while population para- and polyphyly may indicate incipient speciation processes. The genetic diversity and biodiversity of tropical Indo-Pacific sponges appears to be substantially underestimated since the high level of genetic divergence is not necessarily manifested at the morphological level.


Knowledge about the processes that generate and maintain marine biodiversity, and the evolutionary relationships and genetic variation of regional populations, along with assessments of the amount of demographic connection between these populations, are essential for understanding and effectively conserving and managing marine resources [1] such as the highly diverse coral reefs [2]. However, this information is currently lacking in the Indo-Pacific for most coral reef organisms other than fish and corals [3], which is surprising considering their diversity and the significant ecological and economic roles played by this ecosystem [4].

Because marine organisms frequently have high dispersal potential (e.g. [5]), their ranges have often been considered to be vast. However, populations have frequently been found to be highly genetically differentiated [6]. Sometimes they exhibit fine-scale endemism [7], and (cryptic) speciation is a common process in the sea [8]. Two main processes have been invoked as responsible for allopatric speciation in the tropical Indo-West Pacific (IWP): vicariance, where a species' previously coherent geographic range has become fragmented following the formation of a barrier to dispersal; or speciation through a founder effect, where a new population is established by a small number of individuals, often by long-distance dispersal, and subsequent restricted gene flow has led to speciation (reviewed in [9]). Both processes are probably the extremes of a continuum rather than being mutually exclusive [10], but the degree of their interplay remains poorly understood.

This study focuses on tropical marine sponges, a highly diverse, ecologically and economically important, but notoriously understudied, marine invertebrate taxon [11]. Biodiversity analyses of Australasian tropical sponges based on species occurrence data have, for example, revealed quite different trends from those of other marine phyla in the IWP [12]. Latitudinal gradients in sponge diversity were not evident, but various environmental factors were prominent at small spatial (α-diversity) scales. Patterns observed at larger spatial (γ and ε) scales of diversity have been ascribed to biogeographic factors and connectivity (reviewed in [13]). However, investigations based on morphometric data alone did not suffice to unravel the biogeographic factors and connectivity among populations.

The lemon-yellow calcareous sponge Leucetta chagosensis Dendy 1913 (Porifera: Calcarea: Leucettidae) served here as a model to test biogeographic and phylogeographic hypotheses, and to determine the degree of connectivity vs. isolation among populations, with a focus on the SW Pacific. L. chagosensis has a wide Indo-Pacific distribution, ranging from the northern Red Sea to the central Pacific (Moorea, Tuamotu) and from Okinawa (Japan) to Brisbane (Moreton Bay, Australia), and was considered to be a single species throughout this range [14]. However, the validity of this assumption has been challenged by molecular data [15]. L. chagosensis is relatively common and is typically found in shaded habitats below a water depth of 15 m. It is often the dominant macro-sponge in semi-cryptic habitats, e.g., on the Great Astrolabe Reef (Fiji) (Wörheide, pers. observ.). It is viviparous, and as such, is thought to have low dispersal capabilities, as is the case for many other sponges [16]. Initial attempts to use mitochondrial gene sequences (cytochrome oxidase II) to resolve phylogeographic patterns of L. chagosensis in the western Pacific were unsuccessful due to the low variability of this marker [17]. Such low intraspecific mtDNA variation is apparently a general mtDNA feature in sponges and anthozoan cnidarians, presumably due to the low rate of their mtDNA evolution [1820]. A subsequent phylogeographic study, focussing on the SW Pacific, used rDNA internal transcribed spacer (ITS) sequences that provided better resolution and uncovered a deep phylogeographic break on the Great Barrier Reef with distinct northern and southern clades. Each of these GBR clades was more closely related to the Indonesian clade than to each other, and a vicariance scenario of fragmentation and range expansion during and after the last glacial sea level low was suggested as a cause for this structure [15].

The aims of the present study were as follows: to unravel the phylogeographic structure of L. chagosensis over its full geographic range using two nuclear loci; to estimate the degree of connectivity among regional populations; to determine whether vicariance or founder dispersal was responsible for generating the observed patterns; and to test L. chagosensis' taxonomic status as a widespread Indo-Pacific species.



ITS, partial 28S sequences

Our data set of 176 individuals [see Additional file 1] includes samples covering the known geographic extent of L. chagosensis (Fig. 1). The first rDNA fragment used (referred to hereafter as ITS) includes the 3' end of the 18S gene, the full ITS1, the 5.8S gene, and the nearly complete ITS2. It was combined with a fragment from the C2-D2 region of the 28S gene. See Table 1 for general features of both fragments. Although previous studies did not detect any intragenomic polymorphisms (IGPs) in the ITS regions either by SSCP [15] or by the subcloning of amplicons [21], a few IGPs were detected here, both in the ITS and partial 28S gene. All polymorphic sequence types could be resolved into two sequence types per individual. Both rDNA fragments had a high GC content (Table 1, 53–61.9%), in the normal range for sponge rDNA [21, 22].
Figure 1

Geographical map showing the location of collection sites. Distribution of collection sites across the Indo-Pacific, covering all of the known geographic range of Leucetta chagosensis. Map was generated using Map-It [92].

Table 1

Features of sequenced regions.


aligned length (bp)

GC (%)

max. p-distance (%)

nucleotide diversity (π)






28S rDNA










Maximum aligned length (including indels), GC content, maximum uncorrected p-distance in percent, and nucleotide diversity of each sequenced fragment. ITS rDNA: ITS1 and ITS2 plus flanking regions; 28S rDNA: C2-D2 region of the 28S rDNA; ATPSb-iII: second intron of the ATP synthetase beta subunit gene.

rDNA sequence type network

The estimated statistical parsimony network, which had a maximum calculated connection limit of 14 steps (at 95% confidence), is shown in Figure 2. It was congruent in its general topology with a phylogeny estimated by Bayesian inference (BI), which had low branch support (similar to [15]) (not shown). The network showed the sequence types oriented along a central axis, with two main clusters near the top and bottom ends. The top cluster, centered around the most frequent sequence type S14, contained specimens from the Northern GBR, the Coral Sea, Guam, Taiwan, Papua New Guinea, and Sulawesi (shaded in blue). The bottom cluster, centered around the second most frequent sequence type S10, contained specimens from the Central and Southern GBR, the Coral Sea, and PNG (shaded in dark green). Sequence types from the Central South Pacific (shaded in lighter green) form an ambiguous loop in the network that could not be resolved with confidence. Leucetta villosa Wörheide & Hooper 1999, the putative GBR-sister taxon of Leucetta chagosensis, is located within this part of the network. Towards the top cluster, a number of deeply separated branches contain sequence types from the Philippines, the Red Sea, and one specimen from Vanuatu. The last branch, closest to the top cluster, contains sequence types from Japan, Indonesia, and PNG (shaded in lighter blue).
Figure 2

Unrooted network of 72 unique rDNA sequence types. Sequence type network constructed from the concatenated ITS and 28S (C2-D2 region) alignment using statistical parsimony. Numbers with the prefix "S" indicate sequence type number [compare with Additional Table 1] and small black dots denote hypothetical (not observed) intermediate sequence types. The size of the white circles indicates the relative frequency of the sequence type in the data set. Large clades are color coded according to their main geographic distribution. Blue: Northwest Pacific; green: Southern Pacific; orange: Philippines; red: Indian Ocean.

Analysis of Molecular Variance (AMOVA)

AMOVAs were carried out for two subsets of the data: subset 1 contained all samples from the total 15 geographically pooled populations in four regional groups (i.e. the whole data set), and subset 2 contained a more spatially restricted set of individuals from the SW Pacific only. For subset 1 (Table 2), the rDNA genetic variation was hierarchically structured, with about 41% distributed among populations within the groups, about 42% within populations, and only about 17% among the four groups. Fixation indices showed significant high genetic structuring at all hierarchical levels, with the highest structuring within populations and the lowest among groups. For subset 2, less than 3% of the variation was distributed between the two groups, with no significant differentiation. Slightly more than a third of the variation was distributed among populations within the groups, which were highly differentiated, and a little less than two thirds of the variation was distributed within highly differentiated populations.
Table 2

Results of AMOVA analysis (rDNA sequences).


Subset 1: 15 populations, 4 groups

Subset 2: SW Pacific only

Source of variation

Degrees of freedom

Variance components

% of variation

Fixation indices

Degrees of freedom

Variance components

% of variation

Fixation indices

Among groups




φct: 0.16893 *




φct: 0.02918

Among populations within groups




φsc: 0.49489 **




φ sc: 0.40156 **

Within populations




φst: 0.58022 **




φst: 0.41902 **

Hierarchical structuring of molecular variation (AMOVA) within and among a priori defined populations and groups based on rDNA sequences only. Subset 1: 15 populations grouped into four groups; Subset 2: 8 populations from the SW Pacific grouped into two groups. Statistical significance is based on 10,100 permutations in Arlequin: * P < 0.05; ** P < 0.001.

ATP-Synthetase beta subunit intron II (ATPSb-iII)

Sequences and alignment

Both alleles of ATPSb-iII, a phase-0 intron [see Additional file 2], could be resolved from 113 specimens of L. chagosensis. The length of the intron varied between about 750 bp in most individuals to 1089 bp in the population from Okinawa, Japan, due to a very long insertion. The maximum uncorrected p-distance among alleles was 9.57%, with a nucleotide diversity (π) of 0.03524 – about a four-fold increase in variation compared to the rDNA regions sequenced (see Table 1). No significant deviation from neutrality was found (Tajima's D: -1.20920, not significant) and no recombination was detected.

Several indels of various lengths were observed; some were minisatellite repeats restricted to certain populations, e.g. Vanuatu. Length variant heterozygotes differed either in the number of residues in a stretch of homomer thymidines starting at position 301 (max T10) or in a few private indels. After collapsing the 226 alleles of the 113 specimens sequenced, 89 unique alleles remained.

Phylogeny estimation

The estimated unrooted Bayesian phylogeny of the 89 unique alleles is shown in Figure 3 [see also Additional file 3]. An estimated phylogeny using the maximum likelihood optimality criterion showed the same topology, and a statistical parsimony network, constructed using TCS with the 95% connection limit, yielded several disconnected networks (not shown). In the Bayesian phylogeny, six main larger clades were present, all of which were well supported by posterior probabilities (PP, >95%) or bootstrap proportions in Maximum Likelihood analyses (BP, >70). The unrooted tree was characterized by one polytomy, in which the relationships among the four main clades could not be resolved. The general phylogeographic patterns were congruent with the rDNA network (Fig. 2). Several populations were reciprocally monophyletic (Red Sea, Maldives, Philippines, Japan, Samoa, Polynesia), while others were paraphyletic (e.g. Vanuatu) or polyphyletic (Papua New Guinea, Great Barrier Reef, Coral Sea). Alleles from the Philippines were most divergent.
Figure 3

Unrooted phylogeny of 89 unique ATPSb -iII alleles. Unrooted Bayesian phylogeny where highly supported branches with >95% posterior probability and >70% Maximum Likelihood non-parametric bootstrap proportions are indicated by a grey circle; numbers at branches are given if values were less than the above (main clades only). See Additional Figure 3 for detailed values. Numbers with an "A" prefix indicate allele number [compare with Additional Table 1]. Numbers in square brackets after allele numbers indicate how many individuals carried that allele, if more than one. Color code of larger clades is the same as in Figures 2 and 4.

An estimated phylogeny with all 226 alleles had the same topology, and showed that alleles from heterozygous individuals that were subcloned to resolve length variant alleles were found in the same larger clade and were not dispersed among distantly related clades (not shown).

Analysis of Molecular Variance (AMOVA)

AMOVAs were carried out for two subsets of the data, as defined above for the rDNA. For subset 1 (Table 3), genetic variation was, in part, hierarchically structured, with about 9% distributed among the four groups, about 57% among populations within the groups, and about 34% within populations. Fixation indices showed significantly high genetic divergence within populations and among populations within groups, with non-significant among-group differentiation. For subset 2, about 10% of the variation was distributed among the two groups, which were not significantly differentiated. About 46% of the variation was distributed among highly differentiated populations within groups, and about 44% within highly differentiated populations.
Table 3

Results of AMOVA analysis (ATPSb-iII sequences).


Subset 1: 15 populations, 4 groups

Subset 2: SW Pacific only

Source of variation

Degrees of freedom

Variance components

% of variation

Fixation indices

Degrees of freedom

Variance components

% of variation

Fixation indices

Among groups




φct: 0.08923




φct: 0.10153

Among populations within groups




φsc: 0.63065 **




φsc: 0.50910 **

Within populations




φst: 0.66360 **




φst: 0.55894 **

Hierarchical structuring of molecular variation (AMOVA) within and among a priori defined populations and groups based on ATPSb- iII sequences only. Subset 1: 15 populations grouped into four groups; Subset 2: 8 populations from the SW Pacific grouped into two groups. Statistical significance is based on 10,100 permutations in Arlequin: * P < 0.05; ** P < 0.001.

Pairwise Fixation Indices (FST)

All but two pairwise comparisons of genetic differentiation between eight populations in the SW Pacific based on FST values were statistically significant at the 99% confidence level (see Table 4). The lowest significant value observed was 0.1, between the population from the Central GBR and Papua New Guinea. Most other pairwise comparisons had values of more than 0.25, and the highest values observed were above 0.9 (0.94: Sunshine Coast/Brisbane vs. Northern GBR; 0.99: Sunshine Coast/Brisbane vs. Polynesia). Only the pairwise comparison between the population from the Capricorn Section of the southern GBR and the adjacent population from the Sunshine Coast and Brisbane indicated no significant genetic differentiation. Testing for a pattern of isolation-by-distance by performing a Mantel test between genetic (FST) and spatial distances did not reveal any significant relationships (r = 0.273, P = 0.0758) for the whole data set (subset 1: 15 populations), but an analysis of populations along the South Equatorial Current (subset 2: 8 populations) revealed a weak but significant correlation (r = 0.439, P = 0.0326) (not shown). However, some populations did not match this trend; as was the case for some of the pairwise comparisons among populations on the East-Australian coast: S'GBR-Central GBR, S'GBR-N'GBR, Brisbane-Central GBR, Brisbane-N'GBR (compare with Table 4). Calculations using log-transformed distances and log-transformed gene flow parameters (M) revealed significant negative correlations for the 15 population subset (r = -0.3927, P = 0.0144; Fig. 4a) and the eight population subset (r = -0.4177, P = 0.0327; Fig. 4b).
Figure 4

Graph of relationships between genetic and geographic distances. a) Relationships between Leucetta chagosensis log-transformed gene flow parameters (M) and log-transformed distances (km) of the 15 population subset. Thick lines indicate the Reduced Major Axis Regression. Mantel tests revealed a significant major axis regression slope of -0.3451 (P = 0.0144). b) The same relationships for the 8 population subset. Mantel test revealed a significant major axis regression slope of -0.4177 (P = 0.0332).

Table 4

Pairwise FST values of eight populations in the SW Pacific (ATPSb-iII alleles).


Northern GBR

Central GBR

Capricorn Section (S'GBR)

Sunshine Coast & Brisbane

Queensland Plateau, Coral Sea

Papua New Guinea (PNG)


Northern GBR



Central GBR




Capricorn Section (S'GBR)





Sunshine Coast & Brisbane



0.08 ns



Queensland Plateau, Coral Sea







Papua New Guinea (PNG),
















Moorea & Tuamotu (Polynesia)








Pairwise FST values of eight populations in the SW Pacific (ATPSb-iII alleles). Bold values indicate a significance level of <0.01, * a significance level of <0.05 after 10,100 permutations in Arlequin 3.1. ns = non-significant. Populations correspond to subset 2 of Tables 2 and 3.

Combined rDNA and ATPSb-intron analyses

Phylogeny estimation

All three fragments – ITS, the C2-D2 regions of the rDNA, and ATPSb-iII – were obtained for 105 specimens [see Additional file 1]. Among those, 92 unique genotypes were found. The estimated Bayesian phylogeny is shown in Figure 5a. The general topology was congruent with the phylogeny estimated from the ATPSb-iII-only alignment. Several deeply diverging major clades were detected, highly supported by posterior probabilities and maximum likelihood bootstrap values, allowing for an increased resolution compared to the ATPSb-iII-only phylogeny and solving the polytomy among the four larger clades with their main geographic distribution in the Indo-Northwest Pacific (compared with Fig. 3). A closer relationship of the Indian Ocean clades (Fig. 5b, red) to the remaining ones from the Indo-Northwest Pacific (Fig. 5b, blue) was revealed in an unrooted phylogeny, which also suggested a closer relationship of the Philippine clade with the clades from the southern Pacific [see also Fig. 5a and Additional file 4]. The main clade geographically restricted to the southern Pacific (Fig. 5, green) was, itself, substructured into two sister-groups: the south-central Pacific sensu stricto (Fiji, Samoa, Polynesia) and the south-western Pacific, restricted to PNG and Australia's East coast. Several regional populations showed reciprocal monophyly (e.g., the Red Sea, Maldives, Japan, Philippines, Samoa, Polynesia), while others were para- or polyphyletic (e.g., Vanuatu and PNG/GBR). The geographic distribution of two of the main clades (green and blue in Fig. 5a) in the area of the GBR and PNG are also illustrated in Fig. 5c.
Figure 5

Phylogeny of 92 genotypes from the concatenated alignment. a) Bayesian phylogeny of 92 unique genotypes from the concatenated rDNA and ATPSb-iII alignment rooted with the populations from the Indian Ocean. Highly supported branches with >95% posterior probability from Bayesian analysis and >70% Maximum Likelihood non-parametric bootstrap proportions are indicated by a grey circle; numbers at branches are given if the values were less than the above. Terminal branches are collapsed for clarity of presentation. See Additional Figure 4 for a fully resolved phylogeny and detailed support values. Numbers with a "G" prefix indicate genotype number (compare with Additional Table 1). Larger clades are color-coded as in Figures 2 and 3. b) Unrooted phylogeny with the same color-coding of clades as in a) showing only the main branches. c) The map displays the area of apparent overlap of two of the deepest diverging clades that have their main geographic distribution in the NW Pacific (blue) and the Southern Pacific (green) in the area of Papua New Guinea and the Great Barrier Reef (GBR). Map also shows subdivision of GBR populations that follows the sections of the Great Barrier Reef Marine Park, as defined for population genetic analyses (see text for details).

Neighbor-Net analyses

Neighbor-Net analyses of the two separate rDNA and ATPSb-iII partitions of the 92 genotype data set were carried out to explore ambiguities in the data and to evaluate the degree of congruence among both loci [see Additional file 5]. Some ambiguities were detected in the ATPSb-iII Neighbor-Net, especially among genotypes from the southern Pacific (small loops in the network). The rDNA Neighbor-Net showed a much higher degree of incompatible splits (i.e., larger loops). Both were largely congruent with the phylogenies estimated using other methods (Figs. 2 and 3).

A comparison of the two Neighbor Nets revealed that the position of several genotypes was different in the two networks (highlighted in red color). G9 (L. villosa) was found among southern GBR genotypes in the ATPSb-iII Neighbor-Net, whereas they were found among South-Central Pacific genotypes in the rDNA Neighbor-Net, where they were closely related to those from Vanuatu (compare also with Fig. 2). Similarly, G46 (Bali) grouped with genotypes from the north-western Pacific in the ATPSb-iII Neighbor-Net, but was closest to a second genotype from Bali (G47) in the rDNA Neighbor-Net. Also, the position of G14 (Hook Reef, GBR) was different in both networks.

Estimation of migration rates and directions

Migration rates and directions were estimated for two subsets: 1) a three large-population case, where populations from the East coast of Australia, PNG, and the South Pacific were pooled based on geography (Table 5), and a spatially more restricted six population case, also based on geography, containing five populations from Australia plus the one from PNG (Table 6). Populations contained the same specimens as those used for the AMOVA analysis. Very low migration among populations was estimated for subset 1 (Nem: 0.05 to 0.31, Table 5), but for the more detailed subset 2, higher asymmetrical migration was estimated within the GBR (mean Nem up to 3.03, Table 6). Outside the GBR, the only population that was estimated to receive more than one effective immigrant per generation was PNG (Nem up to 1.32 from the Central GBR, Table 6).
Table 5

Migration rates among three pooled populations from Papua-New Guinea (PNG), Australia, and the SW Pacific.


to Australia

to PNG

to the S'Pacific


Θ 0.01483

0.31 (109.68)

0.05 (16.26)


0.12 (31.81)

Θ 0.01127

0.05 (16.68)


0.05 (13.04)

0.05 (16.2)

Θ 0.01126

Θ (Effective population size [Ne] * mutation rate [μ]), migration rates (M), and the number of effective migrants per generation (Nem) among three larger populations in the SW Pacific as estimated by MIGRATE. Populations were pooled: Australia; Papua New Guinea (PNG); Southern Pacific (S'Pacific). For immigration occurring between populations, the upper row is the receiving population and the left column is the broadcasting population. Θ values are given in the diagonal. The first number in the other cells is Nem, derived by calculating (Θ * M)/4, and the numbers in parentheses are migration rates (M). The number of alleles used were (ATPSb-iII/ITS) for Australia: 116/186; PNG: 26/38; and S'Pacific: 34/46, for a total of 176/270.

Table 6

Migration rates among six populations along the Northeast coast of Australia and Papua-New Guinea (PNG).


to N' GBR

to Cen. GBR

to Cap. GBR

to BNE

to Coral Sea

to PNG


Θ 0.00259

3.03 (598.86)

0.37 (274.15)

0.20 (285.07)

0.37 (320.57)

1.15 (231.89)

Cen. GBR

0.16 (254.62)

Θ 0.02026

0.31 (231.56)

0.20 (276.95)

0.31 (269.76)

1.32 (265.87)

Cap. GBR

0.15 (237.15)

1.61 (317.63)

Θ 0.00537

0.29 (409.75)

0.27 (234.11)

1.09 (219.98)


0.13 (202.29)

1.33 (261.62)

0.51 (381.39)

Θ 0.00285

0.25 (215.01)

1.02 (205.91)

Coral Sea

0.14 (220.66)

1.46 (288.97)

0.25 (184.13)

0.19 (262.04)

Θ 0.00458

1.21 (243.72)


0.13 (198.65)

1.45 (285.66)

0.25 (188.54)

0.19 (273.07)

0.28 (248.76)

Θ 0.01982

Θ (Effective population size [Ne] * mutation rate [μ]), migration rates (M), and number of effective migrants per generation (Nem) among six populations along the Northeast Australian coast and Papua New Guinea (PNG), as estimated by MIGRATE. Populations were: N'GBR (Northern Great Barrier Reef [GBR]); Cen. GBR (Central GBR); Cap. GBR (Capricorn Section GBR); BNE (Sunshine Coast & Brisbane); Coral Sea (Osprey, Bougainville and Holmes Reefs); and PNG (Papua New Guinea). For immigration occurring between populations, the upper row is the receiving population and the left column is the broadcasting population. Θ values are given in the diagonal. The first number in the other cells is Nem, derived by calculating (Θ * M)/4, and the numbers in parentheses are migration rates (M). The numbers of alleles used were (ATPSb-iII/ITS) for N'GBR: 22/23; Cen. GBR: 18/16; Cap. GBR: 42/33; BNE: 16/8; Coral Sea: 26/27; and PNG: 28/20, for a total of 152/127.


In this study, the phylogeographic structure of the widely distributed coral reef sponge Leucetta chagosensis was investigated throughout its known Indo-Pacific range, using two unlinked nuclear DNA markers, the rDNA ITS, a fragment of the 28S gene (C2-D2 region), and a novel marker for sponge evolutionary studies, the second intron of the ATP synthetase beta subunit gene (ATPSb-iII). Deep genetic divergences, substantial phylogeographic structure, and substantial amounts of regional isolation with low amounts of migration among regional populations were uncovered, and were congruent across both loci. Based on the present data, we argue that life history traits (low dispersal capabilities) in combination with historical factors (tectonics, sea level fluctuations) are responsible for the diversification of the study taxon over space and time.

Connectivity, or lack thereof

A relationship is expected to exist between the amount of genetic differentiation of allopatric marine populations and their dispersal abilities [23], with low dispersing taxa showing higher genetic structuring and a significant correlation between genetic and geographic distance (isolation-by-distance, IBD) compared to those capable of wide dispersal [24]. In this study, genetic variation was significantly structured among and within populations, but no IBD was detected across the whole set of populations. This suggests that, overall, the colonization process of L. chagosensis was discontinuous, with founder effects occurring when new populations were established [10], and that the very low migration rates were responsible for the lack of genetic cohesiveness. A weak, but significantly positive correlation between FST and geographic distance was detected among populations connected along the South Equatorial Current (SEC), consistent with expectations based on low dispersal capabilities; some populations on the Australian east coast were also highly differentiated despite their close geographic proximity, probably because there was an area of overlap of two of the deepest diverging clades (see below).

Occasional gene flow in a stepping-stone model [25, 26], in which dispersal and genetic exchange occur only between adjacent populations, was suggested by the significant negative slopes of the RMA regression of the logarithmic transformed values (distances and M) for both the whole set of populations and the ones along the SEC. However, the negative slopes (respectively -0.3451 and -0.4177) were shallower than expected for a strict stepping-stone model (-0.5 for a two-dimensional array, [26]) and might indicate some long-distance dispersal [27]. Unfortunately, low geographic coverage of L. chagosensis samples from the South Pacific (e.g. Fiji, Samoa, etc.) prohibited further analysis and inference at this stage.

The low migration rates, below one effective migrant per generation, frequently observed here are less than the minimum required (i.e., 1 Nem) to compensate for genetic drift, and thwart continued genetic divergence of populations [28]. For the GBR, these values are lower than those reported for GBR invertebrates [29, 30]. In a conservation context, where some authors have argued for the use of a threshold of 10 effective migrants per generation as a minimum to counteract local extinction risks [31], the overall low number of migrants encountered here is cause for concern because isolated populations are clearly prone to local extinction. Such a high degree of regional isolation certainly needs to be taken into consideration for the effective design of marine protected areas [1].

Results from the present study suggest that the colonization of new distant habitats by L. chagosensis and the genetic cohesiveness among them cannot be sustained by the dispersal of sexual propagules alone. Occasional long-distance dispersal, either by asexual fragments [32, 33] or by rafting [34] might play a role. To date, there is no direct evidence for long-range dispersal in the studied taxon, but the discovery of a budding specimen of L. chagosensis (Fig. 6; G. Wörheide pers. obs. at Ribbon Reef No. 10 in February 2006) provides evidence that asexual dispersal does occur naturally. Such asexual fragments (buds) constitute fully functional mini-sponges that are probably capable of surviving for a considerable amount of time in the plankton (or rafting on pumice, for example) before potentially colonizing distant reef habitats.
Figure 6

Leucetta chagosensis in its natural habitat. Leucetta chagosensis in its natural habitat at Ribbon Reef #10, Great Barrier Reef, February 2006. Left: A small specimen with one central osculum. Right: An individual producing asexual buds separating from the main animal. The size of both specimens is approximately 4 cm.

Genealogical patterns

The structure of the rDNA sequence type network, with the two most common sequence types at opposite ends of the network ("star-bursts") connected by a longer branch, is reminiscent of a "dumbbell" shaped network [35], which apparently indicates that two long-separated populations each underwent a (recent) expansion. This interpretation would support a scenario where populations were separated during low sea level stands, and subsequently expanded with rising sea levels from refuges on the Queensland Plateau and the shallower shelf south of the GBR [15]. This hypothesis, originally formulated by Davies [36], is also supported by a number of other genetic studies on fish and corals (e.g. [3739]). However, such a vicariance scenario is complicated by the intermediate allopatric populations in the rDNA network, and was not supported by a recent study of calcareous sponges from the same calcarean family Leucettidae, Pericharax heteroraphis, using the same molecular markers (ITS rDNA, ATPSb-iII) [40]. Here, Bentlage and Wörheide observed much less variation in both markers and no phylogeographic (genealogical) structure on the GBR, and attributed this to the retention of ancestral ITS rDNA polymorphisms and a relatively recent expansion after a population bottleneck. This suggests considerable amounts of idiosyncrasy in each species' (demographic) history.

Using the rDNA internal transcribed spacers (ITS) for phylogenetic inferences poses some difficulties. The internal transcribed spacers (ITS) separate the 18S and 28S rDNA genes in the tandemly repeated rDNA cistron [41]. While variation among the multiple ITS copies is normally homogenized by a process called 'concerted evolution' [42], intragenomic polymorphisms (IGP) [21] do occur, e.g., if concerted evolution is slow [43]. The occurrence of potentially paralogous ITS copies can then confound phylogenetic inferences [21, 44]. The few rDNA IGPs detected in this study could all be resolved into two different sequence types, and were always found to be closely related and clade-specific, most likely representing orthologs. Nonetheless, reconstruction of sequence types from rDNA IGPs by haplotype inference should always be preceded by subcloning of rDNA amplicons of several individuals to check the extent of intragenomic polymorphism [21]. We consider the risks of analysing paralogous rDNA sequence types to be minimal in this case, but the inclusion of additional and independent loci is necessary to untangle the true demographic history from locus-specific forces because only genealogical concordance across multiple unlinked loci can elucidate whether phylogeographic breaks are caused by stochasticity or by real barriers to gene flow [45].

To this end, our analysis of ATPSb-iII sequences provided pivotal insight. This new marker showed about five-fold higher substitution rates than the rDNA, and higher phylogeographic resolution was achieved with congruent deeper coalescent patterns across loci. Based on the data presented here, we argue that the high concordance between the two loci is due to population substructuring and putative reproductive isolation, and not to the linked inheritance of both loci (for which no data currently exist).

The large divergence of the reciprocally monophyletic genotypes from the Philippines, with their divergence especially pronounced among the intron alleles, points to a long-standing isolation of the Philippine population. The Philippine islands are an area of complex island-geomorphology that have undergone complicated tectonics, with some islands on crustal fragments originating far south-east of their present position [46]. Several of today's islands have emerged since the late Miocene (about 10 million years ago), and the sampled area at the northern reef of Bohol probably belonged to Proto-Mindanao during the Pliocene [47]. Initial colonization could have happened by long-distance dispersal soon after reef habitats were established. This could then have been followed by persistence when subsequent changes in shoreline constellations due to cycling sea levels reinforced the genetic isolation of this region [29], resulting in the deeply divergent and reciprocally monophyletic clade discovered here. Such apparent founder speciation has been previously described in the terrestrial fauna of the Philippines [47], but might be more widespread in marine biota than previously appreciated [10]. This could also be responsible for generating and maintaining the other reciprocally monophyletic populations discovered here. However, for more in-depth interpretations of the remarkable divergence of the Philippine population, a broader geographic coverage of the whole island group would be required.

Substantial amounts of local endemism (reciprocal monophyly) were also uncovered, pointing to considerable amounts of regional isolation similar to what has been found in some other reef invertebrates [7]. This is also consistent with expectations based on low dispersal capabilities. However, some populations were found to be para- or polyphyletic. Such a pattern of population poly- or paraphyly may be caused by either incomplete lineage sorting or recurrent (recent and historic) gene flow among previously separated clades [35], or both of these factors. These two processes of shared alleles among distinct clades are hard to distinguish on the basis of genetic data alone [48]. Avise [35] proposed that deep coalescence of polyphyletically distributed lineages would point to incomplete lineage sorting, whereas shallow genealogical patterns would be caused by recent gene flow. For L. chagosensis, all polyphyletic populations (e.g. PNG, GBR) showed deep coalescence, as expected for a scenario of incomplete lineage sorting. However, the previous differentiation of geographically restricted reciprocally monophyletic lineages with subsequent re-dispersal into each others' ranges are equally likely to have caused the apparent overlap in geographic distribution of the two deepest diverging clades in the area of PNG and East Australia. Based on the presence of putative hybrids (see below), we favor the latter scenario. However, more extensive geographic sampling from those critical regions would be beneficial to allow additional conclusions.

Under the phylogenetic species concept [49, 50], all reciprocally monophyletic groups uncovered here would represent different recognizable taxa and ESUs (Evolutionary Significant Units [51]). The degree of reproductive isolation among those ESUs cannot be discerned at present, i.e., whether they are from a single, still-interbreeding widespread species, or instead constitute a reproductively isolated sibling (and morphologically cryptic) species. Support for the latter came from the near-fixation of alternative ATPSb-iII alleles in some populations, e.g., between the populations from Brisbane and Polynesia. The occurrence of a different combination of ATPSb-iII alleles and rDNA sequence types in some individuals, uncovered by their different clade-affiliations, is an indication of occurrences of hybridization and incomplete reproductive isolation. One of those hybrids is L. villosa, was described as a closely related sister-species of L. chagosensis based solely on morphological data [14]. This is evidence that those hybrids are also manifested as (morphologically) distinguishable phenotypes.

The temporal aspect

Resolving the temporal frame of clade divergences in L. chagosensis remains a challenge due to the paucity of unequivocally identifiable fossil remains. We also lack information about mutation rates in this taxon. Only from the study of Wörheide et al. [21], who analysed divergences in rDNA ITS sequences in several reciprocally monophyletic populations of Prosuberites laughlini (Demospongiae: Hadromerida) across the Isthmus of Panama, can we derive a mutation rate of roughly 1% ITS rDNA sequence divergence per million years for this demosponge. While the validity of applying this mutation rate estimate to a phylogenetically distant taxon is untested, it does currently represent the only option for a rough estimation of divergence times among L. chagosensis clades. Thus, the deepest divergences among L. chagosensis ITS rDNA sequence types would have happened about 2 million years ago in the late Pliocene, and subsequent divergences would have occurred in the Pleistocene (see Table 1). The late Pliocene, an era with fluctuating sea levels [52], and the resulting differences in shoreline distributions and current (dispersal) patterns, has been estimated to be the time of vicariant speciation in a number of marine taxa, as has the Pleistocene [29].


Deep phylogeographic structure was uncovered in this study, congruent across the two nuclear markers used. There was low gene flow among most regional populations, some of which showed isolation-by-distance effects. Overall, dispersal was indicated in a stepping-stone model with some long-distance genetic exchange – all consistent with expectations from the allegedly low dispersal capability of the taxon. Reciprocally monophyletic populations have long been isolated and constitute putative sibling species, but the degree of reproductive isolation between para- and polyphyletic populations remains to be investigated in more detail. Hybrids displaying mixed genotypes from deeply divergent clades were rare, but one (L. villosa) was morphologically distinguishable. L. chagosensis might therefore be regarded as a complex of closely related (cryptic) species rather than a single widespread species.

However, we cannot conclusively determine at present whether the observed structure was generated solely by vicariance resulting from glacio-eustatic changes during the late Pliocene and Pleistocene, or was simply a result of general dispersal limitation with stochastic long-distance dispersal and founder speciation. However, both processes supposedly contributed with varying degrees to the generation and/or maintenance of the observed phylogeographic structure in each region.

The low migration rates observed here are not sufficient to counteract the continued genetic divergence and regional isolation of populations. Geographically restricted and deeply divergent genealogical lineages of taxa with low dispersal capabilities are therefore prone to local extinction, since immigrants from more viable, distant habitats cannot easily replenish populations depleted of genetic diversity. These results highlight the need for comprehensive biodiversity estimates (i.e., based on true genetic diversity) when considering an effective design and implementation of marine reserves, in order to conserve the complete diversity to promote reef resilience.


Sample collection

Leucetta chagosensis specimens were collected by SCUBA diving on subtidal reef slopes. Specimens from the East Coast of Australia, the Coral Sea, Vanuatu, the Red Sea and Japan were collected by Gert Wörheide (GW); other specimens were received from various colleagues. Due to this fact, and because this taxon can be quite rare (e.g. in the Gulf of Aqaba and the Red Sea), sample sizes were limited in some localities, varying from two in Fiji to 10 in the Philippines individuals per site. Voucher sample numbers and localities are given in Additional file 1. Immediately after collection on deck, macroinvertebrate commensals were carefully removed, if present, and parts of the choanosome were cut into small pieces. Sometimes, small specimens were preserved whole. Samples were preserved in >90% ethanol or in silica-gel after being cut into small pieces, and were stored at -20°C or room temperature until extraction (see also [15, 53]). The external surface was avoided to minimize potential contamination. Voucher samples have been deposited at the Queensland Museum, Brisbane, Australia (indicated by the prefix "QMG" in Additional Table 1 see Additional File 1) or are held by G. Wörheide.

DNA extraction and PCR amplification

Genomic DNA was extracted from the choanosomal tissue using the DNeasy-Tissue Kit (Qiagen) according to the manufacturer's instructions. Internal transcribed spacer (ITS) and partial 28S (C2-D2 region [54]) rDNA sequences were generated according to previously published protocols [15, 21, 55].

The second intron of the ATP synthetase beta subunit gene (ATPSb-iII) was amplified using the following strategy. Initially, degenerate primers (ATPSbf1: 5'-CGT GAG GGH AAY GAT TTH TAC CAT GAG ATG AT-3'; ATPSbr1: 5'-CGG GCA CGG GCR CCD GGN GGT TCG TTC AT-3'; [56]) were used for PCR amplification of four samples from disparate geographic areas (Red Sea, Taiwan, GBR, Tuamotu). BLAST searches [57] were run to confirm that the amplified sequences were of poriferan origin. Exon/intron boundaries were determined by comparison with ATPS beta gene sequences from Genbank. The intron was determined to start about 2.3 kb downstream from the 5'-end of Exon 1 in Drosophila melanogaster (Position 3416; GenBank accession no. X86015). From an initial alignment, nested species-specific primers were designed (ATPSbf2: 5'-TTG TCT TGG ACA AGG AGG GG-3'; ATPSbr2: 5'-TCG TTC ATT TGA CCG TAC AC-3'), which were still located in the exon but were about 10–15 bp closer to the exon/intron boundary. These primers were then used for subsequent PCR amplification of the larger data set.

PCR reactions consisted of 2 mM MgCl2, 16 mM (NH4)2SO4, 67 mM Tris-Cl (pH 8.8 at 25°C), and 0.01% Tween-20. Each 25 μ l reaction included 1.0 μ l of DNeasy-extracted DNA of various concentration, 0.625 μl dNTPs (10 mM each), 1.5 mM of each primer, and 0.25 U Bio-Taq DNA polymerase (Bioline, Luckenwalde). PCR cycling conditions included an initial denaturation of 2 min at 94°C; 38 cycles of 94°C for 20 s, 54°C for 60 s, and 72°C for 50 s; and a final extension step at 72°C for 10 min. Amplicons were purified from Agarose gels using a silica-based method [58]. PCR products were sequenced directly with ABI BigDye Terminator Cycle Sequencing (Version 3.1) on an ABI 3100, and some PCR products were subcloned into pGEM-T (Promega) according to the manufacturer's instructions. A minimum of three positive clones were sequenced in both directions using M13 vector primers to determine if more than two alleles were present, indicating paralogous copies. Amplicons were directly sequenced as above. Because the specific primers were very close to the exon/intron boundaries, the first 42 bases of the 5'- end and the first 53 bases of the 3'- end of the intron were not included in the final alignment.

Sequence assembly and alignment

Double-stranded sequences of novel sequences generated in this study were assembled with CodonCode Aligner [59], automatically aligned using ClustalW [60], and manually inspected and optimized using Se-Al v2.0a11 [61]. Existing rDNA ITS sequence types of Leucetta chagosensis from previous studies [15, 21] (Genbank accession nos. AF458852–AF458870) were added to the alignment. Polymorphic sites in heterozygotes were detected using the 'find mutations' option in the CodonCode Aligner, and coded using IUPAC [62] ambiguous DNA characters. Alleles of heterozygotes that showed no length variation within individuals were resolved manually using a parsimony approach, in an attempt to minimize the number of alleles observed in the data set [63]. Intragenomic polymorphisms in rDNA were resolved into two different sequence types per individual using the same approach, primarily to enable estimation of migration rates in MIGRATE (see below). Due to the low intraspecific diversity of the C2-D2 fragment, rDNA alignments were concatenated for subsequent phylogeographic analyses. All sequences from this study were submitted to the EMBL Nucleotide Sequence Database and can be accessed under the following accession numbers: [EMBL: AM850255–AM850677].

Intron heterozygous indel detection

We attempted to resolve length variant heterozygotes (LVH), which often cause problems when analyzing intron sequences [64], using Single-Stranded Conformational Polymorphism (SSCP) analyses [65], as described in [40]. However, the intron was too long for alleles to be resolved using this method. Therefore, to obtain an overview of the distribution of LVH alleles on the estimated phylogenies, i.e. whether one heterozygote individual harbours alleles distributed in different larger clades on the estimated phylogeny, amplicons of several such heterozygotes from each larger clade were subcloned as above, sequenced, and both alleles were included in the alignment. Due to restricted resources, it was not possible to subclone and sequence both alleles of all heterozygous individuals.

Indels and minisatellite repeats were treated as missing data, and were not recoded for phylogenetic analyses because they represented autapomorphies for regional populations, and as such did not contribute to resolving phylogeographic structure.

Phylogenetic analyses

All sequences were assembled and aligned as described above. Uncorrected p-distances were calculated using Paup*4 [66]. Nucleotide diversities and GC contents were calculated, and tests of neutrality using Tajima's D [67] were carried out in DNAsp v4.1 [68]. Nuclear loci are subject to recombination, potentially confounding phylogeny estimation [69]. Consequently, recombination rates and events were estimated using RDP 2.0 [70], employing the default options with all methods implemented in the program suite (RDP, GENECONV, Bootscan, Chimera, SciScan). The program MODELTEST version 3.7 [71] was used to find the model of DNA substitution that best fit the rDNA and intron data. The best-fit model selected for each sequenced region by hierarchical Akaike information criterion (AIC) [72] (ITS rDNA: TIM+G, partial 28S rDNA: HKY+I, intron: TrN+G) was used for a subsequent Bayesian phylogeny inference (BI) using MrBayes 3.1.2 [73] with default priors and mixed models.

Phylogenies were estimated from the combined ITS/partial 28S rDNA alignment, the intron alignment, and a combined ITS/28S/intron alignment containing only specimens from which all three fragments could be obtained. For the latter, consensus sequences of ATPSb-iII alleles were used and concatenated with both rDNA fragments into one alignment. Alignments were collapsed to contain only unique sequence types/alleles in COLLAPSE 1.2 [74] for phylogeny estimation, and intron sequences were additionally analyzed using the full number of detected alleles.

Three independent runs with one cold and seven heated Markov chains each per analysis were performed simultaneously in MrBayes 3.1.2 until the average standard deviation of split frequencies between the three runs dropped below 0.005 (lowered from the default value of 0.01 to improve chain convergence). Analyses were carried out with the MPI-enabled parallel version of MrBayes [75] on a 64-node Linux-cluster at the Gesellschaft für wissenschaftliche Datenverarbeitung Göttingen (GWDG), using one processor for each of the 24 Markov chains per analysis. Batch files are available upon request. All MrBayes analyses were run at least twice to check for the consistency of results. The second analysis was allowed to run longer without a stop value, until the wall time of 48 hours was reached, to check if convergence of chains could be further improved. Trees were sampled every 1000th cycle, with a burn-in of 25% of sampled trees. The appropriateness of burn-in values was determined after the graphical display of likelihood values, and the convergence of chains was evaluated using AWTY [76]. The remaining trees after burn-in were used to generate a 50% majority rule consensus tree, where posterior probabilities for internal branches were indicated by their sample frequency.

For comparison, Maximum Likelihood bootstrap analyses were conducted with GARLI 0.94 [77] using a heuristic search with the default option, i.e. under the GTR model of nucleotide substitution, with gamma distributed rate heterogeneity, a proportion of invariant sites, and 100 bootstrap replicates. Phylogenies estimated from the intron-only and combined rDNA/intron data set were rooted with the Indian Ocean clades, as they represent a different biogeographic region (e.g. [78]). Outgroup rooting with intron sequences from a closely related species was not possible because such sequences were unalignable due to high evolutionary rates, and because the same intron in Pericharax heteroraphis, for example, is about 60% shorter [40].

Due to the low divergence of rDNA sequences, we estimated a statistical parsimony sequence type network [79] using the program TCS 1.21 [80] with the default options. The maximum number of mutational steps that constitute a parsimonious connection between two sequence types was calculated with 95% confidence. An attempt at constructing a nested clade design [81] was not successful because no unambiguous nested design could be constructed due to a central loop that could not be resolved with confidence.

To obtain an overview of the level of ambiguity and phylogeographic congruency within and among loci and to visually explore the different signals contained in the data, the concatenated alignment was divided into rDNA and ATPSb-iII partitions. Each partition contained samples from 92 specimens from which both loci could be sequenced. Both loci were analyzed separately using Neighbor-Net [82], as implemented in SplitsTree4 [83], using uncorrected p-distances.

In brief, Neighbor-Net is a distance-based method that provides an effective tool to visualize and detect conflicting signals or alternative phylogenetic histories. First, a collection of weighted splits is constructed from a distance matrix, and then these splits are represented using a splits graph, where a totally compatible collection of splits would be precisely represented as a tree, but incompatible splits as cycles or boxes. Such incompatibilities might represent, among other things, hybridization, recombination, gene duplication, or in general, conflicting signals in the data. In general, the more (and larger) cycles/boxes connecting operational taxonomic units (OTUs) in a split graph, the more incompatibilities of splits and incongruences exist in the data.

Population genetic analyses

Analysis of molecular variance (AMOVA) [84] was conducted to estimate the significance of population structure at several hierarchical levels. Pairwise fixation indices (Fst) [85] and their significance were calculated for intron alleles only to estimate population differentiation because true heterozygotes cannot be distinguished from intragenomic variation in the rDNA cistron (see [86]), however, an AMOVA was carried out using the rDNA sequence types. Calculations were carried out in ARLEQUIN 3.1 [87] using 10,100 random permutations for significance tests. All analyses were run at least twice to check for the consistency of results. Two sets of analyses were run for each data set separately, because Arlequin does not allow for unequal sample sizes among loci. First, sample localities were pooled into 15 geographic populations. Pooling of sample localities into populations on the Great Barrier Reef (GBR) followed the sections of the Great Barrier Reef Marine Park: N'GBR, Central GBR, Capricorn (GBR), Brisbane, Queensland Plateau; Taiwan, Guam, Okinawa, Philippines, Indonesia; Maldives, Red Sea; and PNG, Samoa/Fiji/Vanuatu, Polynesia. The 15 populations were grouped into four regional groups (separated by semicolons in the previous list: Australia, NW Pacific, S'Pacific, and Indian Ocean). Due to low sample sizes from some archipelagos in the SW Pacific (Fiji, Samoa) in the intron data set, those localities were pooled with samples from Vanuatu as one population in order to increase the statistical power. To enable comparison between the two loci, the same geographic population structure was employed for the rDNA data set, where more sequences were available from those archipelagos. A second analysis included only eight populations from the SW Pacific (Northern GBR, Central GBR, Capricorn, Brisbane, Queensland Plateau, PNG, Samoa/Fiji/Vanuatu, and Polynesia), grouped into two groups (Australia, S'Pacific), where sampling was more comprehensive both in terms of geography and sample sizes. A Mantel test was carried out using the Isolation by Distance Web Service [88] to test for correlations between spatial and genetic (Fst) distances [89], also using gene flow M calculated as (1/FST-1)/4 [26]. The significance of the slope of the reduced major axis (RMA) regression was assessed by 30,000 randomizations.

Because of the acknowledged difficulties of using FST's to estimate gene flow [90], we also took a coalescent-based approach to estimate the parameters of populations [91], such as theta (θ = 4Ne μ; where Ne is the effective number of individuals and μ is the mutation rate in mutations per generation), migration rates (M), and the number of effective migrants per generation (Nem). We used the coalescent-based Markov-Chain Monte Carlo method implemented in MIGRATE 2.1.7, which explicitly takes into account historical processes and asymmetrical migration/gene flow [91]. The parallelized version of MIGRATE was compiled to run on the LINUX-cluster of the GWDG (see above), requesting several processors for each run depending on the numbers of replicates. Several independent runs were conducted to check for convergence, the consistency of results, and the shape of posterior distributions. The Bayesian search strategy was optimized for the following settings (parm files available on request): recorded genealogies [a]: 100,000; increment (record every x genealogy [b]: 500; and the number of concurrent chains [c]: 4, resulting in 2*108 visited (sampled) genealogies [a*b*c] with 10,000 discarded trees per chain (burn-in). For MIGRATE analyses, each diploid sponge individual was represented by two alleles/sequence types in the input file, whether it was homozygous or heterozygous. rDNA sequences showing IGPs were resolved into two sequence types using the parsimony approach outlined above, since MIGRATE does not allow ambiguously (IUPAC) coded nucleotide sites.

Due to the large geographic distances between some populations without intermediate sample localities, (e.g. Red Sea-Maldives-IWP) and the low sample sizes of some, we focussed our attention on estimating the migration rates among populations along the East-Australian coast, Papua New Guinea, and the southern Pacific only.



We are sincerely grateful to all colleagues who helped to obtain samples from remote places, among them Lori and Pat Colin, John N.A. Hooper, Tomoki Kase, Scott Nichols, Gustav Paulay, Olaf Schlegel, Peter Schupp, Nicole de Voogd, and Nerida Wilson. John N.A. Hooper's funding from Natural Products Discovery at Griffith University and the Queensland Museum (all Brisbane, Australia) enabled biodiversity surveys of the Great Barrier Reef (GBR) and Coral Sea Territories of Australia. The Great Barrier Reef Marine Park Authority is gratefully acknowledged for permitting the fieldwork on the GBR (Permit nos G98/142, G98/022, G00/638, G06/16547.1). We are also grateful to the Egyptian Environmental Affairs Agency (EEAA), especially Mohammed Fouda, for permitting fieldwork in Egypt and to Alexander Keck and Christian Alter for their support during fieldwork in Egypt. We are grateful for critical comments by and inspiring discussions with members of the Molecular Geobiology Lab of the Courant Research Center Geobiology, Göttingen, as well as John N.A. Hooper, John Benzie, and four anonymous reviewers who significantly improved earlier drafts of the manuscript. The use of MIGRATE was stimulated during GW's participation in the Workshop on Molecular Evolution in Woods Hole (USA) in 2005; GW would like to thank Peter Beerli for many subsequent helpful tips on using his program. This study was supported by funds from the Australian Biological Resources Study (ABRS), the German Research Foundation (DFG), and the European Marie-Curie project HOTSPOTS (contract MEST-CT-2005-020561). All experiments carried out in this study were in compliance with German and Australian laws.

Authors’ Affiliations

Courant Research Center Geobiology, Georg-August-Universität Göttingen
Institut für Biochemie und Biologie, Evolutionsbiologie/Spezielle Zoologie, Universität Potsdam


  1. Palumbi SR: Marine Reserves and Ocean Neighborhoods: The Spatial Scale of Marine Populations and Their Management. Annual Review of Environment and Resources. 2004, 29: 31-68. 10.1146/ Article
  2. Hellberg M: Footprints on water: the genetic wake of dispersal among reefs. Coral Reefs. 2007, 26: 463-473. 10.1007/s00338-007-0205-2.View Article
  3. Rocha L, Craig M, Bowen B: Phylogeography and the conservation of coral reef fishes. Coral Reefs. 2007, 26: 501-512. 10.1007/s00338-007-0261-7.View Article
  4. Moberg F, Folke C: Ecological goods and services of coral reef ecosystems. Ecological Economics. 1999, 29: 215-233. 10.1016/S0921-8009(99)00009-9.View Article
  5. Avise JC: Conservation genetics in the marine realm. Journal of Heredity. 1998, 89: 377-382. 10.1093/jhered/89.5.377.View Article
  6. Palumbi SR, Hedgecock D: The life of the sea: Implications of marine population biology to conservation policy. Marine Conservation Biology: The Science of Maintaining the Sea's Biodiversity. Edited by: Norse EA, Crowder LB. 2005, Washington: Island Press, 33-46.
  7. Meyer CP, Geller JB, Paulay G: Fine scale endemism on coral reefs: archipelagic differentiation in turbinid gastropods. Evolution. 2005, 59: 113-125.View ArticlePubMed
  8. Knowlton N: Sibling species in the sea. Annual Review of Ecology and Systematics. 1993, 24: 189-216. 10.1146/ Article
  9. Heads M: Towards a panbiogeography of the seas. Biological Journal of the Linnean Society. 2005, 84: 675-723. 10.1111/j.1095-8312.2005.00466.x.View Article
  10. Paulay G, Meyer C: Diversification in the tropical Pacific: comparison between marine and terrestrial systems and the importance of founder speciation. Integrative and Comparative Biology. 2002, 42: 922-934. 10.1093/icb/42.5.922.View ArticlePubMed
  11. Kelly-Borges M, Valentine C: The sponges of the tropical island region of Oceania: A taxonomic status review. Marine and Coastal Biodiversity in the Tropical Island Pacific Region. Vol. 1 Species Systematics and Information Management Priorities. Edited by: Maragos JE, Peterson MNA, Eldredge LG, Bardach JE, Takeuchi HE. 1995, Honolulu: Program on Environment, East-West Center, 83-120.
  12. Hooper JNA, Kennedy JA, Quinn RJ: Biodiversity 'hotspots', patterns of richness and endemism, and taxonomic affinities of tropical Australian Sponges (Porifera). Biodiversity and Conservation. 2002, 11: 851-885. 10.1023/A:1015370312077.View Article
  13. Wörheide G, Solé-Cava AM, Hooper JNA: Biodiversity, molecular ecology and phylogeography of marine sponges: patterns, implications and outlooks. Integrative and Comparative Biology. 2005, 45: 377-385. 10.1093/icb/45.2.377.View ArticlePubMed
  14. Wörheide G, Hooper JNA: Calcarea from the Great Barrier Reef. 1: Cryptic Calcinea from Heron Island and Wistari Reef (Capricorn-Bunker Group). Memoirs of the Queensland Museum. 1999, 43: 859-891.
  15. Wörheide G, Hooper JNA, Degnan BM: Phylogeography of western Pacific Leucetta 'chagosensis' (Porifera: Calcarea) from ribosomal DNA sequences: implications for population history and conservation of the Great Barrier Reef World Heritage Area (Australia). Molecular Ecology. 2002, 11: 1753-1768. 10.1046/j.1365-294X.2002.01570.x.View ArticlePubMed
  16. Maldonado M: The ecology of the sponge larva. Canadian Journal of Zoology. 2006, 84: 175-194. 10.1139/z05-177.View Article
  17. Wörheide G, Degnan BM, Hooper JNA: Population phylogenetics of the common coral reef sponges Leucetta spp. and Pericharax spp. (Porifera: Calcarea) from the Great Barrier Reef and Vanuatu. Abstracts, 9th International Coral Reef Symposium, Bali, October 2000. 2000, 23-
  18. Shearer T, van Oppen MJH, Romano SL, Wörheide G: Slow mitochondrial DNA sequence evolution in the Anthozoa (Cnidaria). Molecular Ecology. 2002, 11: 2475-2487. 10.1046/j.1365-294X.2002.01652.x.View ArticlePubMed
  19. Wörheide G: Low variation in partial cytochrome oxidase subunit I (COI) mitochondrial sequences in the coralline demosponge Astrosclera willeyana across the Indo-Pacific. Marine Biology. 2006, 148: 907-912. 10.1007/s00227-005-0134-y.View Article
  20. Hellberg ME: No variation and low synonymous substitution rates in coral mtDNA despite high nuclear variation. BMC Evolutionary Biology. 2006, 6: 24-10.1186/1471-2148-6-24.PubMed CentralView ArticlePubMed
  21. Wörheide G, Nichols S, Goldberg J: Intragenomic variation of the rDNA internal transcribed spacers in sponges (Phylum Porifera): implications for phylogenetic studies. Molecular Phylogenetics and Evolution. 2004, 33: 816-830. 10.1016/j.ympev.2004.07.005.View ArticlePubMed
  22. Erpenbeck D, Cleary DFR, Voigt O, Nichols SA, Degnan BM, Hooper JNA, Wörheide G: Analysis of evolutionary, biogeographical and taxonomic patterns of nucleotide composition in demosponge rRNA. Journal of the Marine Biological Society of the United Kingdom. 2007, 86: 1607-1614.
  23. Grantham BA, Eckert GL, Shanks AL: Dispersal potential of marine invertebrates in diverse habitats. Ecological Applications. 2003, 13: S108-S116. 10.1890/1051-0761(2003)013[0108:DPOMII]2.0.CO;2.View Article
  24. Kirkendale LA, Meyer CP: Phylogeography of the Patelloida profunda group (Gastropoda: Lottidae): diversification in a dispersal-driven marine system. Molecular Ecology. 2004, 13: 2749-2762. 10.1111/j.1365-294X.2004.02284.x.View ArticlePubMed
  25. Kimura M, Weiss GH: The Stepping Stone Model of Population Structure and the Decrease of Genetic Correlation with Distance. Genetics. 1964, 49: 561-576.PubMed CentralPubMed
  26. Slatkin M: Isolation by Distance in Equilibrium and Non-Equilibrium Populations. Evolution. 1993, 47: 264-279. 10.2307/2410134.View Article
  27. Hellberg ME: Stepping-stone gene flow in the solitary coral Balanophyllia elegans: equilibrium and nonequilibrium at different spatial scales. Marine Biology. 1995, 123: 573-581. 10.1007/BF00349236.View Article
  28. Mills L, Scott , Allendorf FW: The One-Migrant-per-Generation Rule in Conservation and Management. Conservation Biology. 1996, 10: 1509-1518. 10.1046/j.1523-1739.1996.10061509.x.View Article
  29. Benzie JAH: Genetic structure of coral reef organisms: ghosts of dispersal past. American Zoologist. 1999, 39: 131-145.View Article
  30. Ayre DA, Hughes TP: Climate change, genotypic diversity and gene flow in reef-building corals. Ecology Letters. 2004, 7: 273-278q. 10.1111/j.1461-0248.2004.00585.x.View Article
  31. Vucetich JA, Waite TA: Is one migrant per generation sufficient for the genetic management of fluctuating populations?. Animal Conservation. 2000, 3: 261-266. 10.1111/j.1469-1795.2000.tb00111.x.View Article
  32. Battershill CN, Bergquist PR: The influence of storms on asexual reproduction, recruitment, and survivorship of sponges. New Perspectives in Sponge Biology. Edited by: Rützler K. 1990, Washington D.C.: Smithsonian Institution Press, 397-403.
  33. Zilberberg C, Sole-Cava AM, Klautau M: The extent of asexual reproduction in sponges of the genus Chondrilla (Demospongiae: Chondrosida) from the Caribbean and the Brazilian coasts. Journal of Experimental Marine Biology and Ecology. 2006, 336: 211-220. 10.1016/j.jembe.2006.05.010.View Article
  34. Thiel M, Haye PA: The ecology of rafting in the marine environment. III. Biogeographical and evolutionary consequences. Oceanography and Marine Biology: An Annual Review. Edited by: Gibson RN, Atkinson RJA, Gordon JDM. 2006, 44: 323-429.
  35. Avise J: Phylogeography: The history and formation of species. 2000, Cambridge MA: Harvard University Press
  36. Davies P: Evolution of the Great Barrier Reef. Australian Geologist. 1994, 92: 21-24.
  37. van Herwerden L, Doherty PJ: Contrasting genetic structures across two hybrid zones of a tropical reef fish, Acanthochromis polyacanthus (Bleeker 1855). Journal of Evolutionary Biology. 2006, 19: 239-252. 10.1111/j.1420-9101.2005.00969.x.View ArticlePubMed
  38. Doherty PJ, Planes S, Mather P: Gene Flow and Larval Duration in Seven Species of Fish from the Great Barrier Reef. Ecology. 1995, 76: 2373-2391. 10.2307/2265814.View Article
  39. Smith-Keune C, van Oppen M: Genetic structure of a reef-building coral from thermally distinct environments on the Great Barrier Reef. Coral Reefs. 2006, 1-10.
  40. Bentlage B, Wörheide G: Low genetic structuring among Pericharax heteroraphis (Porifera: Calcarea) populations from the Great Barrier Reef (Australia), revealed by analysis of rDNA and nuclear intron sequences. Coral Reefs. 2007, 26: 807-816. 10.1007/s00338-007-0267-1.View Article
  41. Coleman AW: ITS2 is a double-edged tool for eukaryote evolutionary comparisons. Trends in Genetics. 2003, 19: 370-375. 10.1016/S0168-9525(03)00118-5.View ArticlePubMed
  42. Elder JF, Turner BJ: Concerted evolution of repetive DNA sequences in eukaryotes. The Quarterly Review of Biology. 1995, 70: 297-320. 10.1086/419073.View ArticlePubMed
  43. Wei XX, Wang XQ, Dong DY: Marked intragenomic heterogeneity and geographical differentiation of nrDNA ITS in Larix potaninii (Pinaceae). Journal of Molecular Evolution. 2003, 57: 623-635. 10.1007/s00239-003-2512-8.View ArticlePubMed
  44. Leo NP, Barker SC: Intragenomic variation in ITS2 rDNA in the louse of humans, Pediculus humanus: ITS2 is not a suitable marker for population studies in this species. Insect Molecular Biology. 2002, 11: 651-657. 10.1046/j.1365-2583.2002.00367.x.View ArticlePubMed
  45. Kuo C-H, Avise JC: Phylogeographic breaks in low-dispersal species: the emergence of concordance across gene trees. Genetica. 2005, 124: 179-186. 10.1007/s10709-005-2095-y.View ArticlePubMed
  46. Hall R: The plate tectonics of Cenozoic SE Asia and the distribution of land and sea. Biogeography and geological evolution of SE Asia. Edited by: Hall R, Holloway JD. 1998, Leiden: Blackhuis Publishers, 99-132.
  47. Steppan SJ, Zawadski C, Heaney LR: Molecular phylogeny of the endemic Philippine rodent Apomys (Muridae) and the dynamics of diversification in an oceanic archipelago. Biological Journal of the Linnean Society. 2003, 80: 699-715. 10.1111/j.1095-8312.2003.00274.x.View Article
  48. Bull V, Beltran M, Jiggins CD, McMillan WO, Bermingham E, Mallet J: Polyphyly and gene flow between non-sibling Heliconius species. BMC Biology. 2006, 4: 11-10.1186/1741-7007-4-11.PubMed CentralView ArticlePubMed
  49. Knowlton N, Weight LA: Species of marine invertebrates: a comparison of the biological and phylogenetic species concepts. Species. The units of biodiversity. Edited by: Claridge MF, Dawah HA, Wilson MR. 1997, London: Chapman & Hall, 199-219.
  50. Baum D: Phylogenetic species concepts. Trends In Ecology & Evolution. 1992, 7: 1-2. 10.1016/0169-5347(92)90187-G.View Article
  51. Moritz C: Defining "evolutionary significant units" for conservation. Trends in Ecology and Evolution. 1994, 9: 373-375. 10.1016/0169-5347(94)90057-4.View ArticlePubMed
  52. Haq BU, Hardenbol J, Vail PR: Chronology of Fluctuating Sea Levels since the Triassic. Science. 1987, 235: 1156-1167. 10.1126/science.235.4793.1156.View ArticlePubMed
  53. Wörheide G: The reef cave dwelling ultraconservative coralline demosponge Astrosclera willeyana Lister 1900 from the Indo-Pacific. Micromorphology, ultrastructure, biocalcification, isotope record, taxonomy, biogeography, phylogeny. Facies. 1998, 38: 1-88. 10.1007/BF02537358.View Article
  54. Wuyts J, De Rijk P, Van de Peer Y, Winkelmans T, De Wachter R: The European Large Subunit Ribosomal RNA Database. Nucleic Acids Research. 2001, 29: 175-177. 10.1093/nar/29.1.175.PubMed CentralView ArticlePubMed
  55. Usher KM, Sutton DC, Toze S, Kuo J, Fromont J: Biogeography and phylogeny of Chondrilla species (Demospongiae) in Australia. Marine Ecology Progress Series. 2004, 270: 117-127. 10.3354/meps270117.View Article
  56. Jarman SN, Ward RD, Elliot NG: Oligonucleotide primers for PCR amplification of coelomate introns. Marine Biotechnology. 2002, 4: 347-355. 10.1007/s10126-002-0029-6.View ArticlePubMed
  57. BLAST: Basic Local Alignment and Search Tool. [http://​www.​ncbi.​nlm.​nih.​gov/​BLAST/​]
  58. Boyle JS, Lew AM: An inexpensive alternative to glassmilk for DNA purification. Trends in Genetics. 1995, 11: 8-10.1016/S0168-9525(00)88977-5.View ArticlePubMed
  59. Codon Code Corp. [http://​www.​codoncode.​com]
  60. Thompson JD, Higgins DG, Gibson TJ: CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, positions-specific gap penalties and weight matrix choice. Nucleic Acids Research. 1994, 22: 4673-4680. 10.1093/nar/22.22.4673.PubMed CentralView ArticlePubMed
  61. Rambaut A: Se-Al: Sequence Alignment Editor. 1996, [http://​tree.​bio.​ed.​ac.​uk/​software/​seal/​]
  62. International Union of Pure and Applied Chemistry. [http://​www.​iupac.​org]
  63. Clark AG: Inference of haplotypes from PCR-amplified samples of diploid populations. Molecular Biology and Evolution. 1990, 7: 111-122.PubMed
  64. Creer S, Pook CE, Malhotra A, Thorpe RS: Optimal intron analyses in the trimeresurus radiation of Asian pitvipers. Systematic Biology. 2006, 55: 57-72. 10.1080/10635150500431213.View ArticlePubMed
  65. Orti G, Hare MP, Avise JC: Detection and isolation of nuclear haplotypes by PCR-SSCP. Molecular Ecology. 1997, 6: 575-580. 10.1046/j.1365-294X.1997.00212.x.View ArticlePubMed
  66. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods) Version 4. [http://​paup.​csit.​fsu.​edu/​]
  67. Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123: 585-595.PubMed CentralPubMed
  68. Rozas J, Sanchez-DelBarrio JC, Messeguer X, Rozas R: DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics. 2003, 19: 2496-2497. 10.1093/bioinformatics/btg359.View ArticlePubMed
  69. Posada D, Crandall KA, Holmes EC: Recombination in Evolutionary Genomics. Annual Review of Genetics. 2002, 36: 75-97. 10.1146/annurev.genet.36.040202.111115.View ArticlePubMed
  70. Recombination Detection Program (RDP). [http://​darwin.​uvigo.​es/​rdp/​rdp.​html]
  71. Posada D, Crandall KA: MODELTEST: Testing the model of DNA substitution. Bioinformatics. 1998, 14: 817-818. 10.1093/bioinformatics/14.9.817.View ArticlePubMed
  72. Posada D, Buckley TR: Model selection and model averaging in phylogenetics: advantages of akaike information criterion and bayesian approaches over likelihood ratio tests. Systematic Biology. 2004, 53: 793-808. 10.1080/10635150490522304.View ArticlePubMed
  73. Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.View ArticlePubMed
  74. Collapse: Describing haplotypes from sequence alignments. [http://​darwin.​uvigo.​es/​software/​collapse.​html]
  75. Altekar G, Dwarkadas S, Huelsenbeck JP, Ronquist F: Parallel Metropolis coupled Markov chain Monte Carlo for Bayesian phylogenetic inference. Bioinformatics. 2004, 20: 407-415. 10.1093/bioinformatics/btg427.View ArticlePubMed
  76. AWTY: A system for graphical exploration of MCMC convergence in Bayesian phylogenetic inference. [http://​ceb.​csit.​fsu.​edu/​awty]
  77. Zwickl DJ: Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets under the maximum likelihood criterion. The University of Texas at Austin. PhD thesis. 2006, [http://​www.​zo.​utexas.​edu/​faculty/​antisense/​garli/​Garli.​html]
  78. Benzie JAH: Major genetic differences between crown-of-thorns starfish (Acanthaster planci) populations in the Indian and Pacific Oceans. Evolution. 1999, 53: 1782-1795. 10.2307/2640440.View Article
  79. Templeton AR, Crandall KA, Sing CF: A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data: III. Cladogram estimation. Genetics. 1992, 132: 619-633.PubMed CentralPubMed
  80. Clement M, Posada D, Crandall K: TCS: a computer program to estimate gene genealogies. Molecular Ecology. 2000, 9: 1657-1660. 10.1046/j.1365-294x.2000.01020.x.View ArticlePubMed
  81. Templeton AR: Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history. Molecular Ecology. 1998, 7: 381-397. 10.1046/j.1365-294x.1998.00308.x.View ArticlePubMed
  82. Bryant D, Moulton V: Neighbor-Net: An Agglomerative Method for the Construction of Phylogenetic Networks. Molecular Biology and Evolution. 2004, 21: 255-265. 10.1093/molbev/msh018.View ArticlePubMed
  83. Huson DH, Bryant D: Application of Phylogenetic Networks in Evolutionary Studies. Molecular Biology and Evolution. 2006, 23: 254-267. 10.1093/molbev/msj030.View ArticlePubMed
  84. Excoffier L, Smouse P, Quattro J: Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics. 1992, 131: 479-491.PubMed CentralPubMed
  85. Wright S: The interpretation of population structure by F-statistics with special regard to systems of mating. Evolution. 1965, 19: 395-420. 10.2307/2406450.View Article
  86. van Oppen MJH, Wörheide G, Takabayashi M: Nuclear markers in evolutionary and population genetic studies of scleractinian corals and sponges. Proceedings of the 9th International Coral Reef Symposium, Bali. Edited by: Moosa KM, Soemodihardjo S, Soegiarto A, Romimohtarto K, Nontji A, Soekarno Suharsono. 2002, Jakarta: Ministry for Environment, Indonesian Institute of Sciences, International Society for Reef Studies, 1: 131-138.
  87. Excoffier L, Laval G, Schneider S: Arlequin ver. 3.0: An integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online. 2005, 1: 47-50.PubMed Central
  88. Jensen JL, Bohonak AJ, Kelley ST: Isolation by distance, web service. BMC Genetics. 2005, 6: 13-10.1186/1471-2156-6-13. v.3.14, [http://​ibdws.​sdsu.​edu/​]PubMed CentralView ArticlePubMed
  89. Wright S: Isolation by distance. Genetics. 1943, 28: 114-138.PubMed CentralPubMed
  90. Bossart JL, Prowell DP: Genetic estimates of population structure and gene flow: limitations, lessons and new directions. Trends in Ecology & Evolution. 1998, 13: 202-206. 10.1016/S0169-5347(97)01284-6.View Article
  91. Beerli P, Felsenstein J: Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. Proceedings of the National Academy of Science of the USA. 2001, 98: 4563-4568. 10.1073/pnas.081068098.View Article
  92. Map-It. [http://​woodshole.​er.​usgs.​gov/​mapit/​]


© Wörheide et al; licensee BioMed Central Ltd. 2008

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://​creativecommons.​org/​licenses/​by/​2.​0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.