Phylogeography of a pearl oyster (Pinctada maxima) across the Indo-Australian Archipelago: evidence of strong regional structure and population expansions but no phylogenetic breaks

Abstract

This study investigates the genetic structure and phylogeography of a broadcast spawning bivalve mollusc, Pinctada maxima, throughout the Indo-West Pacific and northern Australia. DNA sequence variation of the mitochondrial cytochrome oxidase subunit I (COI) gene was analysed in 367 individuals sampled from nine populations across the Indo-West Pacific. Hierarchical AMOVA indicated strong genetic structuring amongst populations (Φ_ST = 0.372, P < 0.001); however, sequence divergence between the 47 haplotypes detected was low (maximum 1.8% difference) and no deep phylogenetic divergence was observed. Results suggest the presence of genetic barriers isolating populations of the South China Sea and central Indonesian regions, which, in turn, show patterns of historical separation from northern Australian regions. In P. maxima, historical vicariance during Pleistocene low sea levels is likely to have restricted planktonic larval transport, causing genetic differentiation amongst populations. However, low genetic differentiation is observed where strong ocean currents are present and is most likely due to contemporary larval transport along these pathways. Geographical association with haplotype distributions may indicate signs of early lineage sorting arising from historical population separations, yet an absence of divergent phylogenetic clades related to geography could be the consequence of periodic pulses of high genetic exchange. We compare our results with previous microsatellite DNA analysis of these P. maxima populations, and discuss implications for future conservation management of this species.

biogeography, Coral Triangle, marine conservation, mtDNA, population genetics

Introduction

The Indo-Australian Archipelago is a globally significant reservoir of marine and coral reef biodiversity and is one of the world's most species-rich regions (e.g. Roberts et al., 2002). Marking the junction between both the Pacific and the Indian Ocean, and the tectonic convergence of the Australian and Eurasian continental plates, the region is characterized by geographical and oceanographic complexity (Gordon & Fine, 1996). As a possible reflection of these physical intricacies and their effect on the evolution and ecology of species, phylogenetic studies on marine and coral reef organisms of the Indo-Australian Archipelago have shown that this region is of unquestionable biogeographical and evolutionary importance for marine biodiversity (Benzie, 1998; Briggs, 1999, 2005; Barber et al., 2000; Bellwood & Wainwright, 2002). However, it is still uncertain whether broad-scale mechanisms affecting population dynamics in this region are pervasive across a collective of species, or if these processes impact the population structure of each organism differently as a result of varying ecologies and life histories (Carpenter et al., 2011).

The incorporation of molecular genetic analyses into population studies throughout the Indo-Australian Archipelago has revealed otherwise undetectable factors influencing population dynamics in this region (reviewed by Carpenter et al., 2011). For species that are broadly distributed across both Indian and Pacific Ocean basins, the Sumatra–Java island chain of the Indo-Australian Archipelago recurrently delineates a sharp intra-specific genetic break (e.g. Lavery, Moritz & 1996,Fielder, 1996; Benzie, 1998; Williams & Benzie, 1998). It is inferred that the Indo-Australian Archipelago has functioned as a restrictive barrier to gene flow over evolutionary timescales, particularly during historical periods of lowered sea level, driving a genetic divergence between Indian and Pacific Ocean populations. For some species, however, the Indo-Australian Archipelago has not been so impervious to historical dispersal, with several studies indicating an absence of deep phylogenetic divergence between regions, implying regular genetic exchange between ocean basins (Lessios, Kane & Robertson, 2003; Uthicke & Benzie, 2003).

Patterns of population differentiation across the Indo-Australian Archipelago are not only restricted to the axes of Indian–Pacific separation. During periods of lowered sea levels (most recently during the Late Pleistocene, ∼17 000 years ago), a complex bathymetry of expansive continental shelf, extended steep ridges and deep troughs within the archipelago itself has isolated ocean basins (e.g. Celebes and Sulu Seas) by exposing physical land barriers between populations that could once freely exchange migrants (see Voris, 2000). Several examples of sharp phylogenetic breaks have been observed between populations separated by small geographical distances, which have most likely been due to such separation of smaller ocean basins (Barber et al., 2002; Lourie, Green & Vincent, 2005; Timm & Kochzius, 2008). Population separation due to Pleistocene sea level fluctuations has also been cited as a major contributor to genetic disjunctions either side of the Torres Strait, where genetic exchange has been limited due to an exposed Sahul Shelf between northern Australia and Papua New Guinea (Chenoweth et al., 1998; Gopurenko & Hughes, 2002; Lukoschek, Waycott & Marsh, 2007). In addition to the effects historical isolation and separation of oceanic basins may have had on the formation of genetic structuring, species within this region are also impacted by oceanographic factors, such as immense passage of water conveyed by the Indonesian Throughflow current (ITF, Fig. 1), which moves water from the Pacific to Indian Oceans. Such factors can play a significant role in shaping population genetic patterns in the region by either facilitating or restricting passive larval transport between populations (Barber et al., 2002; Kochzius & Nuryanto, 2008).

Figure 1.

Map showing Pinctada maxima sampling locations throughout the Indo-Australian Archipelago. Shaded area indicates shoreline during Pleistocene low sea level stands (100 m below present-day level, following Voris, 2000). Black arrows indicate major ocean currents across the region, and dotted arrows show seasonally reversing currents (simplified from Wyrtki, 1961; Gordon & Fine, 1996). ITF, Indonesian Throughflow Current; NEC, North Equatorial Current; NECC, North Equatorial Counter Current; Sulu., Suluwesi; HA, Hainan Island; V, Vietnam; WP, West Papua; B, Bali; A, Aru; WA, Western Australia; TS, Torres Strait; EA, East Australia; SI, Solomon Islands.

In addition to its evolutionary and ecological significance, the highly productive marine ecosystems of the Indo-Australian Archipelago also provide a vital economic resource to surrounding nations through small-scale fisheries and tourism (Burke, Selig & Spalding, 2002). Unfortunately, the high incidence of unsustainable commercial practices and environmental degradation has placed this region under considerable threat of over-exploitation (Burke et al., 2002), resulting in the need to assign high priority to the implementation of conservation efforts and resource management practices throughout the region (Roberts et al., 2002). For commercially exploited marine organisms, identifying and understanding factors influencing population dynamics is a critical step towards their appropriate management and sustainable utilization; however, the evident complexity of population structuring throughout the Indo-Australian Archipelago requires continued investigation in order to better understand how historical and contemporary physical processes affect marine organisms in the region. Although understanding generalized phylogenetic patterns across this region is still considered to be in its infancy, continued molecular genetic investigations are required to provide support to broad-scale patterns that are beginning to emerge (Carpenter et al., 2011).

The silver-lipped pearl oyster, Pinctada maxima, is commercially valued for the production of ‘South Sea’ pearls and is widespread throughout the tropical and sub-tropical regions of the Indo-Australian Archipelago. An extended planktonic larval stage of 17–24 days (Rose & Baker, 1994) provides potential for P. maxima to broadly disperse via strong ocean currents; however, microsatellite DNA variation suggests dispersal and gene flow capabilities are limited, resulting in population genetic structuring throughout its distribution (Benzie & Smith-Keune, 2006; Lind et al., 2007). The phylogenetic history of P. maxima remains unknown. Deep phylogenetic divergence is important to identify for long-term genetic resource management strategies (Moritz, 1994), yet is difficult to detect using microsatellite markers because of size homoplasy and potentially high mutation rates. By adopting a phylogeographical approach, it is possible to ascertain whether the contemporary structure of P. maxima populations shown by microsatellites is reflective of a potential genetic footprint borne as a result of historical biogeographical influences, and will further strengthen understanding of historical and current population dynamics in this commercially significant species. Using sequence information from mitochondrial DNA, this study investigates population genetic aspects of P. maxima throughout a biogeographically complex and ecologically important region. Through the potential identification of genetically differentiated populations, this will allow more informed management of a commercially significant resource throughout the Indo-Australian Archipelago, and further increase understanding of the mechanisms of historical biogeographical influences on the genetic structure of marine organisms in the region.

Materials and Methods

Tissue sampling

Tissue samples were taken from 367 P. maxima individuals from nine populations throughout the Indo-Australian Archipelago, spanning China, Vietnam, Indonesia, tropical Australia, and the Solomon Islands (Fig. 1). Local divers collected live oysters from naturally occurring oyster beds, and foot, mantle or muscle tissue samples were excised and preserved in 70–80% ethanol.

DNA extraction, PCR conditions and haplotype detection

Preserved tissue was digested in a CTAB buffer with 20 mg mL⁻¹ proteinase K for 1–3 h at 55 °C, followed by a phenol/chloroform/isoamyl alcohol purification protocol to extract total genomic DNA (gDNA) (Sambrook, Fritsch & Maniatis, 1989). gDNA was quantified by comparison with DNA concentration standards after agarose gel electrophoresis using the ImageJ 1.33 software package (Wayne Rasband, 2004) and resuspended in ddH₂O to a concentration of 5 ng μL⁻¹.

Individual P. maxima were assessed for genetic variability using heteroduplex/temperature gradient gel electrophoresis (TGGE) analysis (Campbell et al., 1995; Elphinstone & Baverstock, 1997) to screen for sequence mutations of a 680-bp region of the mitochondrial cytochrome oxidase subunit I (COI) gene [using primers COIeF: 5′ATAATGATAGGAGGRTTTGG3′ and a GC-clamped reverse primer, COIeR-GC: 5′CGCCCCGCCGCGCCCCGCGCCCGTCCCGCCGCCCCCGCCCGCTCGTGTRCTACRTCCAT3′ (Arndt et al., 1996)]. The implementation of a GC-clamp to the reverse strand primer was to reduce complete denaturation of PCR fragments during TGGE, and can also improve the detection of single-base changes in DNA fragments across denaturing gradients (Sheffield et al., 1989). Polymerase chain reactions (PCR) were conducted in 15-μL volumes with reagent concentrations of 1× PCR buffer containing 1.5 mM MgCl₂ (QIAGEN), 0.2 mM each dNTP, 0.5 μM each primer, 0.03 U μL⁻¹Taq DNA polymerase (Qiagen) and 0.2–0.4 ng μL⁻¹ of DNA template. Thermocyler conditions were as follows: an initial denaturation step of 94 °C for 2 min, followed by 30 cycles of 94 °C for 1 min, 45 °C for 1 min, and 72 °C for 90 s; followed by a final elongation step of 72 °C for 5 min. Each individual was heteroduplexed (HD) with a reference sample (see below) and subjected to TGGE on 5% polyacrylamide gels (8 M urea, 1× ME buffer) following the HD/TGGE procedures outlined by Campbell et al. (1995). A melt transition of ∼48 °C was determined from a perpendicular gel with a temperature gradient of 20–60 °C (300 V, 1.5 h), and subsequent parallel gels were run across a temperature gradient of 28–58 °C (300 V, 2.5 h). Three candidate reference samples were tested on a subsample of ten individuals to determine which reference could best resolve heteroduplex bands. Haplotype 1 was the best resolving reference and thus used for all subsequent HD/TGGE. To aid haplotype scoring accuracy, positive (i.e. reference-known haplotype HD) and negative (i.e. reference-reference HD) control samples were run on each gel. HD banding patterns were visualized through silver staining.

Representative samples (2–4 per haplotype) from each putative haplotype resolved by TGGE were amplified with COIeF and COIeR (no GC-clamp: 5′GCTCGTGTRCTACRTCCAT3′) primers using the same PCR conditions as above. PCR products were then purified by centrifugation through Sephadex™ (Sigma-Aldrich) columns before being sequenced by Macrogen Inc. (Seoul, Korea). Where uncertainty in scoring TGGE haplotypes by eye was encountered (e.g. abnormal sample migration due to an air-bubble or imperfection in a gel), both samples were re-amplified and sent for sequencing. In all instances of uncertain scoring, sequence data confirmed both samples to be the same haplotype as indicated on the gel.

Statistical analyses

Sequence data were proofread and aligned by eye using SEQUENCHER v4.5 (Gene Codes Corporation) before being subjected to further phylogenetic analyses. Haplotype diversity and nucleotide diversity (Nei, 1987) were estimated for each population using ARLEQUIN v3.1 (Excoffier, Laval & Schneider, 2005). To test the selective neutrality of mutational differences observed across populations, and to investigate for indications of rapid population expansions, Tajima's D-test statistic and Fu's F-test were generated using ARLEQUIN v3.1. Further investigation for genetic signals of population equilibrium was performed using Harpending's raggedness index (HRI), and the sum of squared deviations (SSD) between the observed and expected distribution of pairwise sequence mismatches, calculated using ARLEQUIN v3.1. A hierarchical analysis of molecular variance (AMOVA) was undertaken using ARLEQUIN v3.1 to determine whether spatial partitioning of genetic variation was present on a regional (Φ_CT) and/or population level (Φ_ST) (based on 1000 permutations). Pairwise Φ_ST values were also estimated between all populations, and the significance of population differentiation was determined after 1000 permutations. Using pairwise Φ_ST values together with geographical location of each population, Monmonier's maximum difference algorithm (Monmonier, 1973) was implemented using BARRIER v2.2 (Manni, Guerard & Heyer, 2004) to provide a computational geometric approach toward identifying putative genetic boundaries across the distribution of P. maxima. In association with spatial information defined by population locations, the algorithm defines a path of the greatest rate of change in a given distance measure (in this case pairwise Φ_ST) based on Delaunay triangulation and Voronoi tessellation (Manni et al., 2004). A Mantel test to determine the statistical significance of the relationship of pairwise Φ_ST and the geographical distance between two populations was implemented with GENALEX v6 (Peakall & Smouse, 2006), where geographical distance between sample locations was measured as the shortest distance by water using GOOGLE EARTH (Google Inc.). According to likelihood ratio tests implemented in MODELTEST (Posada & Crandall, 1998), the sequence evolution model of Tamura & Nei (1993) with a predicted gamma distribution of 0.069 (TrN + G) best fitted the data set and was used in AMOVA and subsequent phylogenetic analyses. To investigate the phylogenetic relationship amongst haplotypes, minimum spanning trees based on Tamura & Nei's (1993) genetic distance were generated in ARLEQUIN v3.1, and bootstrap maximum parsimony and neighbour-joining trees (based on 10 000 bootstrap replicates) were generated in PAUP* v4.0b10 (D. L. Swofford, Florida State University).

Results

Genetic diversity and sequence phylogeny

From 367 individual pearl oyster samples, TGGE and sequencing identified 47 haplotypes (GenBank Accession numbers: JQ990784–JQ990830, see Table 1). The 635-bp mtDNA COI fragment analysed exhibited 45 variable sites (7%), of which 14 (2%) were parsimony informative. The most frequently observed haplotype was present in samples from every population except the Solomon Islands, whilst the second most abundant haplotype was far more restricted in its distribution, observed predominantly throughout south-east Asian locations (Table 1). Populations from central Indonesia (i.e. Bali, West Papua) exhibited generally higher genetic diversity at both haplotype and nucleotide levels compared with other populations, with the exception of eastern Australia and Solomon Islands (Table 2). Samples from the Solomon Islands showed greatest haplotype and nucleotide diversity, although this result must be interpreted with caution given only five individuals were sampled from this location. A mean pairwise sequence divergence of 1.0% (maximum 1.8%) indicated that haplotypes were closely related, which is visually evident from a minimum spanning haplotype network (Fig. 2). This was supported by phylogenetic analyses, with maximum parsimony and neighbour-joining trees giving low bootstrap support (< 50%) for the presence of any deep phylogenetic lineages or divergence amongst sequences (data not shown).

Figure 2.

Pinctada maxima COI haplotype network. Circles represent different haplotypes, with relative size proportionate to its observed frequency across populations of the Indo-Australian Archipelago. Each haplotype is separated by a single mutational change, with crossed connections indicating the number of mutational differences between linked haploypes; colours indicate proportional distribution of a haplotype amongst each population. Minor text labels indicate haplotype name.

Table 1.

Frequencies of mtDNA cytochrome oxidase I haplotypes observed in Pinctada maxima populations across the Indo-Australian Archipelago; GenBank accession numbers refer to haplotype sequence database entries available at http://www.ncbi.nlm.nih.gov/GenBank/

Haplotype	GenBank accession no.	Solomon Islands	East Australia	Torres Strait	West Australia	Aru	Bali	West Papua	Vietnam	Hainan Island	Total
1	JQ990784	–	7	11	36	10	12	4	14	25	119
2	JQ990785	–	–	1	–	–	18	23	2	–	44
3	JQ990786	–	–	25	8	2	–	1	–	–	36
4	JQ990787	–	–	–	4	26	–	–	–	–	30
5	JQ990788	–	–	–	–	4	5	14	–	–	23
6	JQ990789	–	4	9	–	–	3	2	–	–	18
7	JQ990790	–	–	3	2	4	1	–	–	–	10
8	JQ990791	–	–	5	1	2	–	–	–	–	8
9	JQ990792	–	–	–	2	4	–	–	–	–	6
10	JQ990793	–	–	–	–	3	–	3	–	–	6
11	JQ990794	–	–	–	–	–	–	5	–	–	5
12	JQ990795	–	–	–	–	–	–	–	5	4	9
13	JQ990796	–	–	–	–	–	4	–	–	–	4
14	JQ990797	–	–	–	3	–	–	–	–	–	3
15	JQ990798	–	–	–	–	–	2	–	–	–	2
16	JQ990799	–	–	–	–	–	2	–	–	–	2
17	JQ990800	–	–	–	–	–	1	1	–	–	2
18	JQ990801	–	–	–	–	–	–	2	–	–	2
19	JQ990802	–	–	–	–	–	–	2	–	–	2
20	JQ990803	–	–	–	–	–	–	2	–	–	2
21	JQ990804	2	–	–	–	–	–	–	–	–	2
22	JQ990805	–	2	–	–	–	–	–	–	–	2
23	JQ990806	–	2	–	–	–	–	–	–	–	2
24	JQ990807	–	2	–	–	–	–	–	–	–	2
25	JQ990808	–	1	–	–	–	–	–	1	–	2
26	JQ990809	–	–	–	–	–	–	–	2	1	3
27	JQ990810	–	–	–	–	–	1	–	–	–	1
28	JQ990811	–	–	–	–	–	1	–	–	–	1
29	JQ990812	–	–	–	–	–	1	–	–	–	1
30	JQ990813	–	–	–	–	–	1	–	–	–	1
31	JQ990814	–	–	–	–	–	–	1	–	–	1
32	JQ990815	–	–	–	–	–	–	1	–	–	1
33	JQ990816	–	–	1	–	–	–	–	–	–	1
34	JQ990817	1	–	–	–	–	–	–	–	–	1
35	JQ990818	1	–	–	–	–	–	–	–	–	1
36	JQ990819	1	–	–	–	–	–	–	–	–	1
37	JQ990820	–	1	–	–	–	–	–	–	–	1
38	JQ990821	–	–	–	–	–	–	–	1	–	1
39	JQ990822	–	–	–	–	–	–	–	1	–	1
40	JQ990823	–	–	–	–	–	–	–	1	–	1
41	JQ990824	–	–	–	–	–	–	–	1	–	1
42	JQ990825	–	–	–	–	–	–	–	–	1	1
43	JQ990826	–	–	–	–	–	–	–	–	1	1
44	JQ990827	–	–	–	–	–	–	–	–	1	1
45	JQ990828	–	–	–	–	–	–	–	–	1	1
46	JQ990829	–	–	–	–	–	–	–	–	1	1
47	JQ990830	–	–	–	–	–	–	–	–	1	1
All		5	19	55	56	55	52	61	28	36	367

Table 2.

Genetic diversity statistics and tests of selective neutrality for Pinctada maxima populations across the Indo-Australian Archipelago

Population	Location	N	Nh	Np	h ± SE	π ± SE	Fu's Fs	Tajima's D	SSD	HRI
1. Hainan Island	19.66°N, 112.00°E	36	9	8	0.51 ± 0.10	0.001 ± 0.001	−6.16**	−1.72*	0.002	0.073
2. Vietnam	9.08°N, 105.25°E	28	9	8	0.73 ± 0.08	0.002 ± 0.001	−3.97**	−1.03	0.013	0.101
3. West Papua	1.13°N, 130.54°E	61	13	17	0.80 ± 0.04	0.004 ± 0.002	−2.60	−0.86	0.012	0.038
4. Bali	8.32°S, 114.92°E	52	13	16	0.82 ± 0.04	0.004 ± 0.002	−3.08	−0.83	0.026	0.078
5. Aru	6.43°S, 134.63°E	55	8	13	0.74 ± 0.05	0.004 ± 0.002	1.41	0.30	0.072**	0.219**
6. Western Australia	19.29°S, 119.75°E	56	7	13	0.57 ± 0.07	0.003 ± 0.001	0.01	−1.13	0.440**	0.360
7. Torres Strait	10.97°S, 143.18°E	55	7	9	0.73 ± 0.04	0.003 ± 0.002	0.23	−0.13	0.029	0.105
8. Eastern Australia	24.20°S, 152.84°E	19	7	8	0.82 ± 0.06	0.002 ± 0.002	−2.37**	−0.94	0.006	0.089
9. Solomon Islands	8.30°S, 158.60°E	5	4	8	0.90 ± 0.16	0.006 ± 0.004	−0.26	−1.17	0.132	0.270

n, Number of individuals sampled; N_h, number of haplotypes; N_p, number of polymorphic sites per population; h, haplotype diversity; π, nucleotide diversity; SSD, sum of squared differences of the mismatch distribution; HRI, Harpending's raggedness index (significance as indicated).

**

P < 0.01; *0.01 < P < 0.05.

Demographic history and population structure

Significantly negative values of Fu's F_s statistic were observed in the peripherally located populations of Hainan Island, Vietnam, and eastern Australia, indicating historical demographic expansion events in these regions (Table 2). Tajima's D also indicated a significant deviation from selective neutrality in the Hainan Island population; however, this result can also be interpreted as a consequence of demographic changes (Table 2). A significantly negative Tajima's D statistic can also indicate population size expansion (by highlighting an excess of low-frequency polymorphisms) (Tajima, 1989a, b), which, in the case of the Hainan Island population, seems to be a likely explanation due to its location at the periphery of P. maxima's natural distribution. This is supported by non-significant sum of squared deviations of the mismatch distributions (SSD) in all populations except Aru and West Australia, rejecting the null hypothesis of demographic expansion in the latter two populations (Table 2). A significant partitioning of genetic variation amongst populations was revealed by AMOVA when no regional grouping was inferred (Φ_ST = 0.372, P < 0.001), highlighting genetic structuring throughout the natural range of P. maxima on a population level. On a regional scale, however, hierarchical AMOVA showed that genetic structuring was more complex than simple broad-scale explanations of genetic differentiation (e.g. differentiation across Wallace's Line or Indo-Pacific separation). The partitioning of genetic variation was best explained when populations were grouped into regions coinciding with northern South-East Asia, central Indonesia, and north-west Australia (including Aru) with eastern Australia (Table 3). The realization of strong genetic structuring throughout the region was further confirmed when pairwise Φ_ST values amongst populations showed significant differentiation in all pairwise comparisons except between Hainan Island and Vietnam, and between Bali and West Papuan populations, indicating genetic exchange between these population pairs is enough to prevent genetic differentiation (Table 4). A statistical relationship of geographical distance between populations and pairwise Φ_ST was not observed (Mantel's test, P = 0.515), indicating that an isolation-by-distance (IBD) model of genetic differentiation was not applicable on a broad scale throughout the range of P. maxima, and that other restrictions/barriers to historical gene flow may be in effect across its distribution. The approximate geographical locations of several putative genetic barriers based on pairwise Φ_ST values were identified using Monmonier's algorithm, and are shown in Figure 3. Interestingly, the geographical area bound by putative barriers A and B bears a striking similarity to the path of the ITF (Gordon & Fine, 1996) (Fig. 1), which indicates that oceanographic factors may have a significant influence in shaping the genetic structure of this species.

Figure 3.

Location of likely biogeographical barriers identified with Monmonier's algorithm overlayed on a map of the Indo-Australian Archipelago. Thick black lines indicate paths of the greatest rate of change in pairwise Φ_ST genetic distances (i.e. possible barriers), thin lines show Voronoi tessellation. 1, Hainan Island; 2, Vietnam; 3, West Papua; 4, Bali; 5, Aru; 6, Western Australia; 7, Torres Strait; 8, East Australia; 9, Solomon Islands.

Table 3.

Hierarchical analysis of molecular variance (AMOVA) to partition mtDNA COI gene variation amongst possible regional groupings of Pinctada maxima populations

Hierarchical grouping	Amongst population		Amongst groups	P-value
	ΦST	ΦCT	% Variation
East–West partition [HA – V – WP – B – WA] [A – TS – EA – SI]	0.393*	0.072	7.23	0.30
Sahul Shelf partition (Torres Strait) [HA – V – WP – B – A – WA] [TS – EA – SI]	0.369*	0.007	0.00	0.45
Wallace's Line partition [HA – V – B] [WP – A – WA – TS – EA – SI]	0.329*	0.065	6.54	0.68
Regional partitioning (based on Monmonnier's algorithm) [HA – V] [WP – B] [A – WA – TS – EA] [SI]	0.436*	0.362	36.22*	< 0.001
No groups	0.372*	–	–	< 0.001

HA, Hainan Island; V, Vietnam; WP, West Papua; B, Bali; WA, West Australia; A, Aru; TS, Torres Strait; EA, East Australia; SI, Solomon Islands.

*

P < 0.001.

Table 4.

Pairwise Φ_ST estimates between Pinctada maxima populations, based on the genetic distance of Tamura & Nei (1993) with gamma correction

Population	Hainan Island	Vietnam	West Papua	Bali	Aru	West Australia	Torres Strait	East Australia	Solomon Islands
1. Hainan Island	–
2. Vietnam	0.035	–
3. West Papua	0.541	0.434	–
4. Bali	0.484	0.363	0.011	–
5. Aru	0.241	0.210	0.447	0.407	–
6. Western Australia	0.059	0.061	0.474	0.410	0.157	–
7. Torres Strait	0.271	0.256	0.538	0.492	0.241	0.110	–
8. Eastern Australia	0.166	0.148	0.492	0.445	0.131	0.107	0.192	–
9. Solomon Islands	0.812	0.684	0.283	0.312	0.547	0.683	0.719	0.693	–

Bold values indicate significant differences [P < 0.05 after false discovery rate correction (Benjamin & Hochberg, 1995)]; underlined values become non-significant after correction.

Discussion

Genetic diversity and range-wide gene flow patterns

Based on frequency and sequence analyses of the mitochondrial COI gene, there exists strong regional genetic structure throughout the natural distribution of P. maxima. Moderate to high levels of genetic diversity were observed in P. maxima populations, which is within the range of variability seen across the Indo-Pacific in other marine species for the same gene region (Barber et al., 2002; Benzie et al., 2002; Gopurenko & Hughes, 2002; Uthicke & Benzie, 2003; de Bruyn et al., 2005; Kochzius & Nuryanto, 2008). Note that it is possible that some haplotypes were not detected using the TGGE screening approach; however, this is most likely to occur with haplotype sequences that are highly similar to the reference (i.e. differing by only a few base-pair mutations). Given we were able to successfully discriminate many haplotypes that differed by only a single based-pair mutation (Fig. 2), we are confident that if any haplotypes were missed it would be of minor frequency and likely to have minimal influence on our analyses. Mitochondrial diversity also followed patterns of diversity observed with microsatellite markers, in that P. maxima populations located near range peripheries generally tend to be less genetically diverse (Lind et al., 2007). Previous studies have identified population structure in P. maxima using microsatellite DNA markers (Benzie & Smith-Keune, 2006; Lind et al., 2007), although mtDNA data presented here indicate a much stronger genetic partitioning than previously realized. Nevertheless, this confirms that despite an extended planktonic larval stage (17–24 days) providing the potential to broadly disperse, gene flow throughout the distribution of P. maxima is restricted.

A stepping-stone, IBD model of gene flow implies that more geographically distant populations are less likely to exchange genetic migrants and are therefore more likely to exhibit genetic differences arising from the random processes of genetic drift (Kimura & Weiss, 1964). Such a significant correlation between genetic differentiation and geographical distance of populations (reflecting patterns of IBD) was not observed in P. maxima across the scale of its distribution. Interestingly, this contrasts with patterns observed in microsatellite data generated on the same populations studied here, which showed that nuclear DNA markers exhibit an IBD pattern of differentiation within the South-East Asia region (Lind et al., 2007). In addition, significant genetic differentiation was detected amongst Indonesian populations using mtDNA, which was not detected with nuclear microsatellite markers (Lind et al., 2007). Given that maternally inherited mtDNA is more prone to population perturbations and the influence of genetic drift as a consequence of its smaller effective size (one-quarter the effective size of the nuclear genome), significant differentiation and the failure of IBD patterns from mtDNA may indicate a more substantial imprint of complex oceanographic influences or historical factors such as biogeographical barriers and repeated population expansions/contractions on the genetic structure of P. maxima. Indeed, analysis of gene flow patterns through pairwise Φ_ST values shows strong differentiation amongst some populations (e.g. Torres Strait vs. Bali, West Papua; Hainan Island; Vietnam vs. Solomon Islands) yet little genetic difference is detected between others (e.g. Hainan and Vietnam; Bali and West Papua) (Table 4), suggesting a more complex pattern of regional genetic structure is present. In addition to a smaller effective size, the effect of protandry, i.e. the sequential hermaphrodism from male to female in P. maxima, may also have important contributions towards the exaggerated genetic structure observed in maternally inherited mtDNA markers compared with the previously mentioned studies based on nuclear microsatellite markers. As highlighted by Arnaud-Haond et al. (2003), protandrous hermaphodism in pearl oysters may result in a strong male-biased effective sex ratio, which can further reduce the mtDNA effective population size compared with sex-separated or simultaneous hermaphrodite species. Male-biased sex ratios have been recorded in younger, smaller demographics of Western Australian P. maxima populations [up to dorso-ventral margin (DVM) size of 170 mm], whereas sex ratios amongst larger animals (170–200 mm DVM) are biased towards females (Hart & Joll, 2006). The lack of IBD patterns in mtDNA markers may also be explained by this factor, given that P. maxima migrants may have a greater chance to reproduce as a male than to reach the age to produce female gametes, thus contributing additional stochasticity to effective mtDNA genetic exchange. Apart from the Indonesian populations mentioned previously, however, patterns of genetic differentiation are broadly concordant with those elucidated in P. maxima using microsatellite markers (Benzie & Smith-Keune, 2006; Lind et al., 2007).

Influence of prevailing ocean currents on gene flow patterns

Given that P. maxima is a sessile benthic organism once its planktonic larvae have settled, passive transport via ocean currents is effectively the sole mechanism allowing gene flow across broad geographical regions, and therefore should play a significant role in the genetic structuring of this species. However, throughout Indo-Pacific marine systems the planktonic larval duration or dispersal capability of an organism has been shown to be an unreliable predictor of population connectivity and genetic structuring patterns (Barber et al., 2002; Ovenden et al., 2004; Bay, Crozier & Caley, 2006). Based on particle dispersion models across the north-west shelf of Australia, P. maxima larvae will passively travel up to 60 km from their origin of spawning, although are most likely to settle within ∼30 km due to the prevailing ocean current and tidal conditions in this region (Condie et al., 2006). This level of movement is perhaps less than expected given P. maxima's extended planktonic larval phase of 17–24 days (Rose & Baker, 1994); however, it appears that single generation dispersal of this magnitude is still sufficiently large to maintain high gene flow over evolutionary time amongst populations spanning thousands of kilometres along the north-west Australian coastline (Johnson & Joll, 1993; Benzie, Smith & Sugama, 2003; Benzie & Smith-Keune, 2006).

A high connectivity between Hainan Island and Vietnam populations is suggested by hierarchical AMOVA and pairwise Φ_ST values, yet limited genetic exchange between more southern regions (Table 3). Strong ocean currents in the western South China Sea, particularly along the coastline of Vietnam (Wyrtki, 1961) (Fig. 1), suggest that a high connectivity via passive larval transport on ocean currents would be likely, and is a plausible explanation for the genetic similarity observed between Vietnam and Hainan Island. Additionally, with approximately 10 million m³ s⁻¹ of water flowing via the ITF from the Pacific Ocean towards the Indian Ocean, mostly through the Makassar Strait (Gordon & Fine, 1996; Gordon, Susanto & Vranes, 2003), passive larval transport via the ITF provides an obvious predictor for high gene flow through the constricted Indonesian seaways. Interestingly, Monmonior's algorithm predicts two possible genetic barriers in locations which isolate the South China Sea and surround an area of close resemblance to the path of the ITF (A and B, Fig. 3), with the populations sampled within each of these areas (Hainan Island–Vietnam and Bali–West Papua, respectively) showing a high genetic similarity (Table 4). Evidence for ITF-mediated gene flow patterns has been observed in several other Indo-Pacific invertebrate and fish species (e.g.Barber et al., 2002; Kochzius & Nuryanto, 2008; Timm & Kochzius, 2008) and is consistent with results presented here. The ITF's deflection away from Aru and Western Australia towards the Indian Ocean may also explain the genetic dissimilarity between these populations and those from central Indonesia. The role of the ITF is a significant factor in shaping population genetic patterns of marine species within Indonesia, and mtDNA evidence from this study suggests that the ITF also has a prominent influence on the genetic structuring of P. maxima.

Genetic patterns observed in the Solomon Islands are also intriguing, indicating a large genetic divergence in populations from this region, with haplotypes more closely related to those found predominantly in West Papua/Bali (Fig. 2). It has been observed in two giant clams species (Tridacna maxima and T. gigas) that populations from the Solomon Islands are more genetically similar to populations from the Philippines than those from the Great Barrier Reef, Australia, and highlight the significance of historical dispersal patterns rather than present-day ocean circulation in the formation of genetic patterns across the west Pacific (Benzie &,Williams, 1995, 1997). Although data observed in this study suggest a similar pattern occurring in P. maxima, such conclusions must be drawn with appropriate caution, given only five individuals were sampled from this region.

Genetic impact of Pleistocene ocean basin isolation and exposure of continental land masses

In addition to oceanographic influences, episodes of lowered sea level have contributed to phylogeographical patterns across the Indo-Australian Archipelago through the formation of physical land barriers across the region, blocking passages of gene flow and causing allopatric differentiation between previously (or presently) connected populations. During periods of Pleistocene polar glaciation (most recently ∼17 000 BP), lowered sea levels of up to 120 m below present-day levels left the Sunda Shelf, Sahul Shelf and other shallow seafloor regions exposed, causing separation of ocean basins across the Indo-Malay region (Voris, 2000). Signatures of historical vicariance between ocean basins has persisted in several present-day populations across the Indo-Malay region and northern Australia, where deep genetic divergence has been observed in marine invertebrates (Barber et al., 2002; Gopurenko & Hughes, 2002; Barber, Erdmann & Palumbi, 2006) and fish species (Chenoweth et al., 1998; Lourie et al., 2005; Timm & Kochzius, 2008). A historically isolated South China Sea basin caused by an exposed Sunda Shelf to the west and the Philippines to the east may have contributed to the genetic differences observed in P. maxima from this region compared with other populations across the Indo-Malay archipelago, yet sustained an adequate gene flow to maintain homogeneity between Hainan Island and Vietnam.

The expansive, shallow continental regions of the Sunda Shelf across the Gulf of Thailand, and the Sahul Shelf between Australia and Papua New Guinea (Fig. 1) are also prominent features in marine biogeography within the Indo-Australian Archipelago, and have probably played a significant role in shaping population genetic patterns in P. maxima. Rapid re-colonization of exposed land regions with rising sea levels (particularly the Sunda Shelf, Hanebuth, Stattegger & Grootes, 2000) could heighten genetic differentiation of populations in these regions through founder effects, and is considered a significant contributor to genetic patterns in several marine species throughout the Indo-Malay region (e.g. Arnaud, Bonhomme & Borsa, 1999; Nelson et al., 2000; Lourie et al., 2005; Mahidol et al., 2007). Genetic differentiation from (re)colonization events are particularly relevant to broadcast spawning marine bivalves, where a small number of effective breeders can contribute large proportions of offspring within a generation with a ‘sweepstakes’-like chance of reproductive success (Hedgecock, 1994; Hedrick, 2005), increasing the likelihood of genetic drift in small or newly colonized populations. In P. maxima, Fu's F-test and non-significant SSD estimates (Table 2) indicate that genetic signatures of demographic and/or range expansions are present in the northern populations of Hainan Island and Vietnam, which is consistent with expectations based on Sunda Shelf re-colonization. Fu's F_s also indicates significant signatures of demographic expansion in the peripheral East Australia population (Table 2). Historical population expansions within the Indo-Malay region have also been suggested from genetic patterns in a tropical abalone species (Haliotis asinina, Imron et al., 2007), as well as in mudcrabs (Gopurenko, Hughes & Keenan, 1999) and sea cucumbers (Holothuria nobilis, Uthicke & Benzie, 2003). It must be noted, however, that signatures of range expansions in P. maxima may also be simply due to the peripheral locations of these populations (especially eastern Australia and Hainan Island), which may have seen repeated population expansion and contraction in response to historical fluctuations in environmental conditions in range peripheries.

No deep phylogenetic divergence yet geographical haplotype associations

Phylogenetics studies have shown a strong influence of regional vicariance on the formation of divergent genetic clades across regions of the Indo-Australian Archipelago separated by only hundreds of kilometres (Barber et al., 2002; Kochzius & Nuryanto, 2008; Timm & Kochzius, 2008), with a realization that genetic connectivity can be significantly restricted despite a potential to broadly disperse. In P. maxima, however, the presence of deep phylogenetic divergence was not detected. A lack of deep genetic divergence and the close relationship of haplotype sequences (1% mean divergence between 47 haplotypes) in P. maxima may suggest that its high dispersal potential is occasionally realized and long-distance dispersal events have periodically occurred. Alternatively, phylogenetic patterns seen in P. maxima are also in agreement with patterns observed over much broader scales. It is often observed that populations from north-west Australia are more phylogenetically related to Asian/Pacific clades than those from other Indian Ocean regions, most likely because of positive influences on population connectivity caused by the IFT (Williams & Benzie, 1998; Benzie et al., 2002; Uthicke & Benzie, 2003; Bay et al., 2004). This is likely to also be the case in P. maxima.

The observation of several small clusters of closely related haplotypes found only in Hainan Island, Solomon Islands, and to a lesser extent in Bali and West Papua (Fig. 2), could also be indications that preliminary lineage sorting has occurred across the distribution of P. maxima (Avise, 1994).

A comprehensive review of phylogenetic studies conducted within the Coral Triangle (the marine regions surrounding much of Indonesia, Malaysia, the Philippines, Brunei, Timor L'Este, Papua New Guinea, and the Solomon Islands) has highlighted that several broad phylogenetic patterns amongst various invertebrate and vertebrate marine species are beginning to emerge (Carpenter et al., 2011). The authors suggest broad management units that coincide with concordant phylogenetic patterns or major ocean currents could be suitable if additional studies continue to corroborate with existing studies. A lack of deep phylogenetic divergence in P. maxima is consistent with genetic patterns observed in some species reviewed by Carpenter et al. (2011), which is amongst several other broad-scale phylogenetic patterns throughout the region. According to the authors just mentioned, this suggests that additional demographic data or habitat data would be necessary to support conservation management of this species. The importance of such approaches is highlighted by recent surveys in the Solomon Islands, which reveal depleted P. maxima populations have struggled to recover from previous commercial overexploitation (Hawes et al., 2011).

Conclusions

Pinctada maxima shows strong genetic structure throughout its natural distribution yet low sequence divergence amongst COI mtDNA haplotypes. Historical biogeographical barriers to gene flow are likely to be a main cause of gene frequency-based differentiation in this species, although high genetic connectivity in some regions appears to coincide with present-day patterns of major ocean currents. The lack of deep genetic divergence in P. maxima could be a consequence of periodic pulses of high genetic exchange. Our study did not detect the presence of significant evolutionary divergence across a broad geographical region, suggesting broad-scale conservation management strategies may be appropriate for this species. Further study on other mtDNA or nuclear gene regions would provide greater support to this. Additional investigation on whether genetic differentiation observed in this study (as well as others) is indicative of phenotypic or physiological differences amongst different populations would be necessary for comprehensive conservation management.

Acknowledgements

This work was funded by an Australian Research Council Linkage-Research Grant (LP-0560298). We are grateful to G. McGlauchlin (Western Australia), D. Williams (Coral Sea Pearls, Eastern Australia), T. Nguyen (Vietnam), P. Southgate and A. Wang (Hainan Island), K. Takami (Kazu Pearls, Torres Strait), and Atlas South Sea Pearl Ltd./P.T. Cendana Indopearls (Indonesia) for kindly providing tissue samples from the locations indicated. We thank two anonymous reviewers for helpful comments.