Microsatellites Reveal Genetic Homogeneity among Outbreak Populations of Crown-of-Thorns Starfish (Acanthaster cf. solaris) on Australia’s Great Barrier Reef

Abstract

Specific patterns in the initiation and spread of reef-wide outbreaks of crown-of-thorns starfish are important, both to understand potential causes (or triggers) of outbreaks and to develop more effective and highly targeted management and containment responses. Using analyses of genetic diversity and structure (based on 17 microsatellite loci), this study attempted to resolve the specific origin for recent outbreaks of crown-of-thorns on Australia’s Great Barrier Reef (GBR). We assessed the genetic structure amongst 2705 starfish collected from 13 coral reefs in four regions that spanned ~1000 km of the GBR. Our results indicate that populations sampled across the full length of the GBR are genetically homogeneous (G’_ST = −0.001; p = 0.948) with no apparent genetic structure between regions. Approximate Bayesian computational analyses suggest that all sampled populations had a common origin and that current outbreaking populations of crown-of-thorns starfish (CoTS) in the Swains are not independent of outbreak populations in the northern GBR. Despite hierarchical sampling and large numbers of CoTS genotyped from individual reefs and regions, limited genetic structure meant we were unable to determine a putative source population for the current outbreak of CoTS on the GBR. The very high genetic homogeneity of sampled populations and limited evidence of inbreeding indicate rapid expansion in population size from multiple, undifferentiated latent populations.

1. Introduction

Crown-of-thorns starfish (CoTS; Acanthaster spp.) naturally occur on coral reefs throughout the Indo-Pacific [1,2] While normally found at low densities [3], sporadic population outbreaks of CoTS cause significant localised coral loss and are a major contributor to the ongoing degradation of coral reefs throughout the Indo West-Pacific [4,5,6,7]. Numerous hypotheses have been put forward to explain the occurrence of CoTS outbreaks (reviewed by [1,2,8]), most of which incite an anthropogenic basis for the purportedly recent and increasing occurrence of outbreaks. While their inherent life-history characteristics (most notably their high fecundity [9]) predisposes CoTS to major fluctuations in abundance [10], there are two prominent theories proposed to trigger outbreaks; the larval survival hypothesis suggests that anthropogenic eutrophication of nearshore waters dramatically increases the survival of planktonic larvae [11,12,13], whereas the predator removal hypothesis postulates that overfishing of natural predators has allowed more CoTS to reach sexual maturity [1]. Tests of either hypothesis require improved knowledge of when and where outbreaks start and corresponding research on environmental conditions and population demographics of CoTS at these locations.

Outbreaks of CoTS on Australia’s Great Barrier Reef (GBR) were first documented in 1962 [14], though there are earlier reports of high densities of starfish (which may or may not have constituted an “outbreak”) on the GBR e.g., [15,16]. Since 1962, there have been three additional outbreak episodes on the GBR, starting in 1979, 1993 and 2009. However, each outbreak has followed a reasonably consistent pattern where primary outbreaks were first recorded on mid-shelf reefs between Lizard Island and Cairns (the ‘initiation box’; Fabricius et al. [13]), followed by a wave of ‘secondary outbreaks’ that tend to propagate southwards [17,18]. Biophysical modelling of larval dispersal patterns suggests that reefs within the initiation box are highly connected [19], thereby explaining why outbreaks that initiate in this region inevitably lead to reef-wide outbreaks. However, limited temporal and spatial resolution of monitoring e.g., [20], as well as inevitable delays in responding to new outbreaks mean that it is still unclear where exactly outbreaks arise. It is also unknown whether outbreaks start from a small cluster of reefs within this area or arise simultaneously on widely separated reefs [2]. Resolving the exact timing and location where outbreaks start is important to establish environmental triggers [13] or changes in population demographics [21] that cause outbreaks. This could also lead to improved management and containment strategies to stop outbreaks before they spread.

Although native to the GBR, the spatial distribution of CoTS following an outbreak closely resembles the spread of invasive species and infectious diseases [22]. Historical and observational data have previously identified invasion routes or disease vectors e.g., [23], but direct observations have proven ineffective for CoTS e.g., [20], where it is still unclear whether outbreaks start from a single reef or arise simultaneously from separate locations [24]. The rapid expansion of populations following biological invasions can, however, lead to distinct patterns of genetic structure and diversity. Genetic data have increasingly been used as an indirect method to describe the spread of invasive species [25,26] and infer the relationship between discrete populations or possible migration routes [27]. Elaborating on these patterns, model-based statistics can provide probabilistic estimations of the demographic and genetic history that are necessary to generate observed patterns of genetic structure e.g., [28,29]. Approximate Bayesian computation (ABC) approaches that incorporate the divergence and admixture of populations, as well as changes in population size and structure [30] can provide important information on the likely initiation and spread of species e.g., [31,32].

This study examines genetic diversity and structure, based on sampling of crown-of-thorns starfish during the current outbreak on the GBR. Over four thousand starfish were sampled at 13 reefs spanning ~1000 km. The spatial genetic structure of a CoTS outbreak will depend on the history of the source population(s), the size of the initial population(s), the dispersal of individuals that led to a primary outbreak and successive secondary outbreaks. While the demographic factors contributing to each stage of an outbreak are unclear, a recent review clarifies many aspects of their population dynamics and life-history characteristics [2]. A model-based approach would therefore allow us to take into account the stochasticity of these demographic processes and test multiple scenarios that would have generated the observed spatial population structure of CoTS on the GBR. Specifically, we test whether the outbreak was generated from a single source population in the ‘initiation box’ or multiple populations. If a small number of individuals from a single source population caused a localised primary outbreak, we would expect successive secondary outbreaks to be affected by a single bottleneck, be composed of highly related individuals and have low genetic diversity. If the primary outbreak originated from multiple populations, we would expect multiple bottlenecks, founding a widespread admixed population, and high genetic diversity. Using a model-based Bayesian approach, we explore the likelihood of competing outbreak scenarios against the observed spatial genetic structure of CoTS to determine the most parsimonious origin of primary outbreaks on the GBR.

2. Materials and Methods

2.1. Sample Collection

Crown-of-thorns starfish (Acanthaster cf. solaris) were collected between April 2013 and May 2015 from 13 reefs between Lizard Island (S 14.7; E 145.4) and the Swains reefs (S 22.3; E 1527) (Figure 1). Starfish were collected whilst snorkelling or SCUBA diving (depending on working depths). All starfish were kept alive in 500 L tanks connected to high flow-through seawater systems on live-aboard boats or at the Lizard Island research station for a maximum of 20 h before being processed. Starfish were placed on their aboral surface to remove tube feet for genetic analyses; 5–10 tube feet (depending on size) were removed using fine scissors to cut the tube feet close to their base. Multiple tube feet from each individual were placed together in 2-mL vials of 100% ethanol. Later, a single foot was taken, frozen and transported in dry ice for processing at King Abdullah University of Science and Technology (KAUST). All equipment used during sampling and processing was sterilised using a three-step rinse procedure involving bleach, water and ethanol.

2.2. Microsatellite Genotyping and Locus Characteristics

The genetic diversity amongst sample populations was assessed using 26 previously-described microsatellite loci [33,34,35,36]. Full details of microsatellite development and screening are available in Harrison et al. [36]. DNA extractions were performed from single tube feet following procedures described in the Nucleospin-96 Tissue kit (Macherey-Nagel, Germany), and microsatellites were amplified in four multiplex reactions of 6 or 7 loci. All PCR were performed using the QIAGEN Microsatellite Type-it kit (QIAGEN, Germany) and screened on an ABI 3370xl DNA Analyzer (Applied Biosystems). PCRs were repeated on 96 individuals to estimate locus-specific genotyping error. Individual genotypes were scored in Genemapper v4.0 (Applied Biosystems), and unique alleles were distinguished using marker specific bin sets in the R package ‘msatallele’ [37].

The numbers of genotyped individuals (N), number alleles (Na), observed heterozygosity (Ho) and expected heterozygosity (He) were estimated for each locus using GenAlEx v6.5 [38,39]. The exact test of Hardy–Weinberg and the score test for heterozygote deficiency were performed in Genepop on the web [40,41] alongside Weir and Cockerham’s estimate of F_IS [42]. The probability test for linkage disequilibrium was performed for each pair of loci in Genepop on the web based on 10,000 dememorisations, 5000 batches and 10,000 iterations. Significance levels of 0.05 were adjusted for a given false discovery rate of 10% to account for multiple testing [43]. The presence of null alleles was determined from 700 randomly selected individuals in Microchecker [44] based on the estimator of Brookfield [45].

2.3. Diversity Analysis of Sampled Populations

Allelic richness, observed and expected heterozygosities within populations were calculated in GenoDive 2.0 [46]. To assess population structure, we estimate global and pairwise genetic differences among sampled reefs within regions. Given the high level of heterozygosity in microsatellite markers, the degree of global genetic differentiation among populations was estimated from Hedrick’s standardised fixation index G’_ST [47]. Population structure among reefs was estimated from an analysis of molecular variance (AMOVA) with reefs nested in regions. Significance was assessed from 9999 permutations, and the standard deviation of genetic variance was obtained through jack-knifing over loci. The significance of population differentiation between all pairs of populations was assessed from 9999 permutations over populations. We tested for isolation by distance (IBD) by comparing the pairwise matrix of linearised genetic distance (Hedrick’s G’_ST) and geographic distance (km) between all sampled reefs. Statistical significance was assessed using a Mantel test with 9999 permutations, and significance values are shown.

2.4. Spatial Genetic Clustering

We applied a model-based Bayesian clustering method implemented in Structure v2.3.2 [48,49] to evaluate the most parsimonious allocation of samples to distinct genetic clusters following the method described by Evanno et al. [50]. If population outbreaks of CoTS stem from the successful fertilisation of multiple populations, we would expect the likelihood probability of the data to depart from its expected distribution of a single homogenous population. We performed three short runs for each number of populations K, from K = 1 to 10, and calculated the mean posterior probability for each value of K. Each run assumed population admixture for correlated allele frequencies and no a priori population assignment. The burn-in length was 30,000 followed by a Markov chain Monte Carlo (MCMC) length of 50,000 repetitions. The initial Dirichlet parameter for the degree of admixture ‘alpha’ was fixed to 1.0 in all simulations.

The statistical power to detect genetic structure amongst sampled populations was evaluated in the software Powsim v4.1 [51]. Three tests were performed to determine whether the number of sampled individuals, the number and diversity of loci could detect F_ST values of 0.0010 and 0.0005 and 0.0001. Each test assumed an effective population size of 5000 individuals and divergence times of 10, 5 and 1 generations, respectively. Due to software limitations in the number of alleles per locus, Apl19 was excluded from these analyses. The numbers of dememorisations, batches and iterations were set at 1000, 100 and 1000, respectively. A total of 1000 replicates were run for each test.

2.5. Approximate Bayesian Computation of Source Populations

We developed four scenarios to represent the possible origins and subsequent spread of CoTS outbreaks on the GBR and used an approximate Bayesian computation approach to determine the probability of each scenario to generate observed patterns of genetic diversity. All analyses were developed using DIYABC 2.1.0 [52,53]. Each scenario was defined to represent possible historical events and demographic changes, including possible population bottlenecks, the effective population sizes of ancestral populations, as well the number of individuals that would have contributed to primary outbreaks, the effective population size of sampled populations and possible divergence times. The four scenarios represent: (i) independent primary outbreaks in the northern and southern GBR, which led to separate secondary outbreaks in the northern and southern GBR (Figure 2a); (ii) a single primary outbreak, which led to secondary outbreaks with a common origin (Figure 2b); (iii) divergence of the northern and southern GBR populations followed by independent primary outbreaks in each of the five focal regions and subsequent secondary outbreaks (Figure 2c); and (iv) sequential secondary outbreaks starting from a common population at Lizard island (Figure 2d).

In order to make the ABC approach computationally feasible, we performed all tests on a subset of 68 randomly selected individuals from Lizard Islands, Cooktown, Cairns and Townsville and all 68 individuals from the Swains. All 17 loci were kept in the analysis. Parameters of mutation models of microsatellite loci were drawn from a set of 10,000 simulations generated for each putative scenario. The posterior probability value of each scenario was determined following 100,000 simulations of each scenario. The details of the ABC parameters and analyses are described in Figure S3.

2.6. Inbreeding

The mating of closely-related individuals, or closely-related groups of individuals, can result in an increase in the frequency of homozygote loci and, thus, heterozygote deficiency. The coefficient of inbreeding f is the probability that two alleles at any one locus are identical by descent (IBD) [54] Wright 1943) so that the progeny of close relatives will have many autozygous genotypes. We measured the multilocus heterozygosity (MLH) [55,56] to investigate the occurrence of inbred individuals in our sample. MLH is the proportion of loci at which an individual was heterozygous divided by the mean heterozygosity across all samples. If individuals are in fact inbred, all loci should exhibit signs of excess heterozygosity deficiency, and estimates of MHL from two random sets of loci should be similar. Furthermore, our sample should exhibit lower estimates of MLH than a random population generated from the observed allelic frequencies. We estimated MLH of CoTS sampled and tested the correlation between randomly selected sets of loci in the R package ‘Rhh’ [57]. Ninety-five percent confidence intervals around the mean were estimated from 999 iterations. Simulated genotypes were generated from observed allelic frequencies using R packages ‘gstudio’ [58] and ‘adegenet’ [59], and MLH between groups were compared used a paired t-test. Inbreeding should result in a significant over-representation of individuals with a low MLH when compared to randomly-assorted genotypes.

3. Results

A total of 4082 individual CoTS were collected between April 2013 and May 2015 from 13 reefs between Lizard Island and the Swains reefs (Figure 1). All individuals were genotyped at 26 microsatellite loci, and the final dataset was curated to focus on five focal regions (

Table 1. Sampling locations of 2705 crown-of-thorns starfish collected on the Great Barrier Reef that were then genotyped with 17 microsatellite loci. Table includes collection period, the number of genotyped individuals (N), allelic richness (AR) and the observed (Ho) and expected heterozygosity of each sampled population (Hs).

Region	Sampling Site	Latitude	Longitude	Date of Collection	N	AR	Ho	Hs
Lizard Is	Lizard Island reefs	−14.6916	145.4479	Oct 2013–Feb 2015	385	5.6	0.673	0.685
Lizard Is	MacGillivray Reef	−14.6524	145.4892	Oct 2013	192	5.4	0.675	0.685
Lizard Is	Nth Direction	−14.7445	145.5399	Oct 2013	247	5.6	0.674	0.690
Lizard Is	Sth Direction	−14.8562	145.4825	Oct 2013	182	5.3	0.672	0.684
Cooktown	Emily Reef	−15.6325	145.6511	Feb 2014	278	5.5	0.687	0.688
Cooktown	Endeavour Reef	−15.7823	145.5847	Feb 2014	208	5.6	0.683	0.692
Cooktown	Pickersgill Reef	−15.8838	145.5640	Feb 2014	151	5.6	0.670	0.693
Cooktown	Spitfire Reef	−16.1148	145.6424	Oct 2013	155	5.4	0.673	0.688
Cairns	Arlington Reef	−16.7749	145.9767	Sept 2014	262	5.4	0.670	0.685
Cairns	Hedley Reef	−17.2474	146.4637	Sept 2014	275	5.4	0.666	0.684
Cairns	McCulloch Reef	−17.2996	146.4257	Sept 2014	265	5.4	0.673	0.686
Townsville	Townsville reefs	−18.4303	146.8209	Nov 2014	137	5.5	0.651	0.682
Swains	Swains reefs	−22.3112	152.6720	May 2015	68	5.5	0.685	0.694

). Individuals removed from the data included 642 individuals sampled prior to September 2013 or less than 80 mm in length and 121 individuals that were sampled from reefs where less than 50 individuals were sampled.

Poor DNA quality or the possible presence of PCR inhibitors in DNA extractions led to a large amount of missing data and genotyping error. A training dataset [36] identified four loci (Apl21, Apl36, Apl39 and Maki12) with over 3% genotyping error and five loci (Apl01, Apl25, Apl27, Apl29 and Apl37) with over 30% missing data, which were all discarded prior to analysis. In addition, 600 individuals with four or more missing loci were discarded from further analyses along with 14 duplicate genotypes associated with the presence of missing data.

Of an initial 4082 samples, our final data included 2705 unique individuals from 13 reefs from Lizard Island (n = 1002), Cooktown (n = 692), Cairns (n = 802), Townsville (n = 137) and the Swains (n = 68). Individuals ranged from 90 mm to 510 mm (mean = 282 mm ± 66 mm SD) with no statistical difference in the size of individuals between regions. Amongst all genotyped individuals, 1925 were identified as either mature male or mature female with 1.34 males for every female. This ratio was used to parameterise the ABC models.

3.1. Microsatellite Data

Data presented here contained 17 polymorphic loci with 1.7% missing data (

Table 2. Characteristics of 17 microsatellite used in the analysis of genetic diversity in crown-of-thorns starfish, Acanthaster cf. solaris. The number of genotyped individuals (N), number of alleles (Na), observed (H_O) and expected (H_E) heterozygosity were measured from 2705 unique individual from the Great Barrier Reef. The exact test of Hardy–Weinberg (HWE-p) and a specific test for heterozygosity deficiency (HWE-h) were measured alongside Weir and Cockerham’s measure of F_IS (F). Significant values following correction for multiple comparisons are shown in bold.

Locus name	N	Na	Ho	He	Missing Data (%)	HWE-p	HWE-h	F
AP1	2705	4	0.648	0.657	0.0	0.399	0.159	0.0091
AP12QS	2703	6	0.490	0.514	0.1	0.002	0.001	0.0321
AP654	2689	11	0.685	0.702	0.6	0.389	0.047	0.0089
AP9	2702	23	0.845	0.843	0.1	0.471	0.360	0.0009
Hisayo01	2688	32	0.804	0.829	0.6	0.010	0.003	0.0179
Maki03	2701	28	0.775	0.792	0.1	0.001	0.005	0.0123
Sayo03	2659	9	0.492	0.519	1.7	0.010	0.006	0.0273
AP11QS	2703	11	0.751	0.762	0.1	0.171	0.119	0.0056
AP30QS	2705	9	0.617	0.625	0.0	0.010	0.001	0.0635
AyU03	2704	29	0.885	0.876	0.0	0.276	0.917	−0.0041
Yukina06	2681	32	0.878	0.917	0.9	0.000	0.000	0.0287
AP5QS	2705	7	0.590	0.600	0.0	0.000	0.014	0.0211
Apl07	2440	7	0.390	0.409	9.8	0.040	0.009	0.0269
Apl19	2474	62	0.949	0.965	8.5	0.644	0.010	0.0066
Etsuko01	2657	10	0.772	0.777	1.8	0.001	0.183	0.0045
Sayo01	2704	7	0.214	0.220	0.0	0.095	0.035	0.0148
Apl02	2577	10	0.653	0.669	4.7	0.335	0.032	0.0108

). The mean number of alleles per locus was 17.4 and ranged from four to 62 alleles. Similarly, the average observed heterozygosity was 0.67 and ranged from 0.21 to 0.95. Each locus was tested for departure from Hardy-Weinberg equilibrium, and nine loci showed non-random association of alleles after correction for multiple testing. A specific test for heterozygous deficiency highlighted 12 loci with a higher than expected frequency of heterozygotes after correction for multiple testing. Estimates of Weir and Cockerham’s F_IS were positively skewed with an average F_IS of 0.017 ± 0.004 SE across all loci, which suggests some degree of inbreeding in these populations. Amongst 136 pairwise comparisons, two locus pairs measured linkage disequilibrium after correction for multiple testing, and no locus showed evidence of null alleles.

3.2. Spatial Patterns of Genetic Diversity

The mean allelic richness (AR) and genetic diversity (H_S) within sampled reefs were high (AR = 5.5; H_S = 0.687) and consistent among reefs (Table 1). There was no evidence for genetic differentiation amongst sampled reefs with global estimates of G’_ST = −0.001 (p = 0.948) and no genetic variance amongst regions (F_CT = 0.000 ± 0.000 SD) or amongst populations nested in regions (F_SC = 0.000 ± 0.000 SD;

Table 3. Analysis of molecular variance amongst crown-of-thorns starfish collected from 12 reefs in 4 regions of the Great Barrier Reef. The significance of F-statistics was measured from 9999 permutations.

Source of Variation	Nested in	% var 1	F-stat	F-value	SD 2	CI 2.5%	CI 97.5%	p-Value 3
Within Individual		0.98	Fit	0.020	0.004	0.012	0.028	-
Among Individual	Population	0.02	Fis	0.020	0.004	0.012	0.028	0.000
Among Reefs	Region	0	Fsc	0.000	0.000	0.000	0.000	0.903
Among Region		0	Fct	0.000	0.000	0.000	0.000	0.830

¹ Per cent of total genetic variance; ² Standard deviations of F-statistics were obtained through jack-knifing over loci; ³ 95% confidence intervals of F-statistics were obtained through bootstrapping over loci.

). Pairwise genetic differences amongst populations were not significantly different from zero (

Table 4. Pairwise differentiation between crown-of-thorns starfish collected from 13 reefs of the Great Barrier Reef (lower diagonal). The significance of p-values (upper diagonal) was measured from 9999 permutations.

	Lizard Island	MacGillivray	North Direction	South Direction	Emily	Endeavour	Pickersgill	Spitfire	Arlington	Hedley	McCulloch	Townsville	Swains
Lizard Island		0.210	0.958	0.391	0.721	0.990	0.707	0.302	0.752	0.333	0.826	0.730	0.952
MacGillivray	0.000		0.380	0.749	0.380	0.657	0.988	0.407	0.358	0.210	0.764	0.585	0.634
North Direction	0.000	0.000		0.316	0.662	0.992	0.784	0.158	0.331	0.480	0.700	0.954	0.949
Soutth Direction	0.000	0.000	0.000		0.827	0.764	0.248	0.380	0.622	0.960	0.821	0.980	0.857
Emily	0.000	0.000	0.000	0.000		0.912	0.678	0.570	0.672	0.425	0.845	0.948	0.854
Endeavour	−0.001	0.000	−0.001	0.000	−0.001		0.876	0.551	0.637	0.850	0.889	0.991	1.000
Pickersgill	0.000	−0.001	0.000	0.000	0.000	−0.001		0.349	0.553	0.096	0.484	0.601	0.745
Spitfire	0.000	0.000	0.001	0.000	0.000	0.000	0.000		0.060	0.202	0.712	0.703	0.874
Arlington	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.001		0.717	0.865	0.854	0.739
Hedley	0.000	0.000	0.000	−0.001	0.000	0.000	0.001	0.000	0.000		0.892	0.948	0.952
McCulloch	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000		0.822	0.895
Townsville	0.000	0.000	−0.001	−0.001	−0.001	−0.001	0.000	0.000	−0.001	−0.001	0.000		0.986
Swains	−0.001	0.000	−0.001	−0.001	−0.001	−0.002	−0.001	−0.001	−0.001	−0.001	−0.001	−0.002

), and showed no evidence of isolation by distance (Mantel test: p = 0.9).

There was, however, evidence of low, but significant heterozygote deficiency within the sample. Individual heterozygosity was lower than would be expected in a large, randomly-mating population with a global estimate of G_IS = 0.022 (p = 0.001) and amongst individuals within individuals reefs (F_IS = 0.021 ± 0.004 SD). Such patterns are commonly associated with evidence of the mixing of different source populations (Wahlund effect), inbreeding due to the mating of close relatives or the non-random sampling of a limited number of familial pools. The level of heterozygote deficiency also varied amongst reefs, ranging from −0.004 to 0.048 and an average G_IS of 0.021 ± 0.012 SD amongst reefs.

Using the whole sample of individuals without prior information of sampling location, the population could not be partitioned into independent populations, indicating that the sample has a common origin. The largest mean log-likelihood values of the data were for K = 1 population and decreased with increasing values of K (Figures S1 and S2). The result provided by the analysis of spatial genetic clustering could not unambiguously detect separate groups of individuals in the sample, indicating a homogenous population, and does not support a Wahlund effect as a source of heterozygote deficiency. The power of the tests suggest that the number of sampled individuals, the number of loci and the allelic diversity at these loci were sufficient to detect genetic structure amongst sampled populations. At a level of differentiation of F_ST = 0.0005, the power to detect genetic heterogeneity with 95% confidence was 100%. At the lowest level of differentiation (F_ST = 0.0001), the power was reduced to 53.7%. Analyses were repeated after excluding nine loci that did not meet equilibrium assumptions (AP12QS, Hisayo01, Maki03, Sayo03, AP30QS, Yukina06, AP5QS, Apl07, Etsuko01). Results from these runs were not different from runs with the full dataset.

3.3. Inbreeding

We measured low, but significant levels of heterozygote deficiency in most screened loci, which could indicate the presence of inbred individuals in our sample. However, the mean multilocus heterozygosity correlation was −0.003 (CI 95% −0.025 to 0.019) indicating that the heterozygote deficiency was not significantly different across loci. To confirm that inbreeding was not the cause of heterozygote deficiency in sampled CoTS, we compared the MLH between sampled and simulated genotypes. Amongst all CoTS sampled in the GBR, the average MLH was 0.29 ± 0.11 SD, which was not significantly different from simulated genotypes (0.26 ± 0.10 SD, t = −9.41, df = 2704, p-value = 1).

3.4. ABC Framework

The most likely scenario supported by the ABC analyses was Scenario 4 (Figure 2d), whereby all sampled populations originated sequentially from a single common primary outbreak (Figure S4). However, the outcomes of this scenario could not be distinguished from Scenario 2, whereby sampled populations originated from a single common primary outbreak at the same time point (Figure 2b). Both scenarios were more representative of the observed data than alternative scenarios that incorporated a divergence between the Swains and the northern populations (Figure 2a,c). In all simulated scenarios, ABC analyses were particularly sensitive to the effective population size for the ancestral population (N_anc). The posterior distribution of N_anc was however consistent amongst scenarios, and estimates suggest that the effective population size that would have led to a primary outbreak would be ~5000 individuals (

Table 5. Prior distribution of historic and genetic data associated with evolutionary scenarios in the ABC analysis.

Parameter	Definition	Distribution (Interval, Mean, SD)
Nanc	Effective population size of an ancestral population	Normal (103 to 104, 5 × 103, 103)
Nprim	Effective population size of a primary outbreak	Normal (103 to 104, 4 × 103, 8 × 102)
Nsamp	Effective population size of a sampled population	Normal (5 × 104 to 2 × 105, 8 × 104, 104)
t1	Divergence time of independent outbreaks	Uniform (1 to 10)
db	Foundation time of a primary outbreak	Uniform (1 to 10)
t2	Divergence time of northern and southern populations	Uniform (102 to 1.5 × 103)
Mean μ	Mean mutation rate	Gamma (10−4 to 10−3, 10−4, 2)
Mean P	Mean of the geometric distribution of the number of repeats	Gamma (10−2 to 10−0, 5.5 × 10−1, 3)
Mean μSNI	Mean single nucleotide insertion/deletion mutation rate	Gamma (10−8 to 5 × 10−5, 1.5 × 10−5, 3)

; Figure S3).

4. Discussion

Assessing a species’ dispersal ability and capacity to colonise new habitats is critical for our understanding of their population biology and ecology [60,61], particularly where species have ecological or economic importance. CoTS outbreaks represent the single most important biological disturbance on coral reefs throughout the Indo West-Pacific [62] and often account for up to 50% of coral loss recorded on coral reefs over the last few decades [5,6,7]. For these starfish, knowledge of population structure and the movement of individuals among reefs can greatly influence management decisions e.g., [19], leading to improved detection and understanding of the patterns of outbreaks, as well as prioritisation of reefs for direct intervention (culling) in an attempt to contain outbreaks. In the current study, we investigated the genetic diversity and structure of CoTS on the GBR with the specific intention of identifying the origin and the direction of subsequent spread for current outbreaks apparent at reefs between Cooktown and Townsville, as well as at Swains reefs, in the southernmost portion of the GBR. This would have further important management ramifications, whereby containment of future outbreaks would be most effective by concentrating monitoring and control on the reef(s) where outbreaks initiate. Our results indicate that populations sampled across the full length of the GBR are genetically homogeneous, highly diverse and have no apparent genetic structure. Furthermore, model-based Bayesian analyses showed that the current outbreaking population of CoTS in the Swains is not independent of the outbreak populations in the northern GBR, but share a common origin.

We found no evidence of genetic structure amongst CoTS genotyped from 13 reefs and five regions spanning over 1000 km along the GBR. While the genetic diversity in our sample was high, there was no variance in diversity among reefs. Power analysis confirms that our data were sufficient to detect even very low degrees of genetic differentiation (F_ST > 0.0005) had there been any structure in the sampled population. The results indicate that CoTS from Lizard Island to the Swains are genetically homogenous. There was a small but significant deficiency in heterozygosity amongst individuals that could have arisen from: (i) the successful reproduction among differentiated populations (Wahlund effect); (ii) the mating of close relatives (inbreeding) or (iii) the sampling of close relatives. We could not distinguish independent clusters in our sample, and therefore, a Wahlund effect is unlikely (sensu [63]). Moreover, multilocus genotypes did not provide consistent evidence of inbreeding. We did find evidence that some individuals were highly related, which would result in a small but significant deficit in heterozygosity (data not shown). The strength and number of relationships between individuals from the same or different reefs could not, however, distinguish between alternative explanations for deficiencies in heterozygosity. Furthermore, a spatial autocorrelation analysis did not identify any relationship between the relatedness of individuals and their spatial distribution. Previous sampling during non-outbreak periods reported small but significant genetic structure among latent CoTS populations on the GBR (F_ST = 0.003, [64]), even though they also found limited structure when sampling each outbreak population.

Our ability to determine a putative source population for the current outbreak of CoTS on the GBR is significantly constrained by the limited genetic structure (sensu [65,66]), despite hierarchical sampling and large numbers of CoTS genotyped from individual reefs and regions. Very high levels of homogeneity across the 2705 individual starfish points to the recent and very rapid increase in population size of CoTS across the GBR, which has arisen from either a single source population or multiple undifferentiated populations. A single origin would be expected to generate some inbreeding and highly related individuals, which the data do not support. A more parsimonious explanation, therefore, is that the current outbreak arose almost simultaneously across a number of reefs, but with largely undifferentiated latent populations. This is consistent with reports during the last major outbreak in 1994, whereby increasing CoTS densities occurred almost simultaneously at Lizard Island and several nearby reefs (including Linnet Reef, North Direction Island and Rocky Islet) before being reported on reefs to the south [20]. High levels of gene flow amongst these closely-positioned populations would have likely resulted in admixture and high levels of genetic homogeneity, as recorded in this study. However, given the lack of genetic structure between regions, we are unable to unequivocally state whether the current outbreak did or did not originate in the northernmost sector of the initiation box, nor can we establish the directionality in the spread of individuals between regions. It is also likely that the latent population of CoTS is sufficiently large and sufficiently diverse to prevent genetic drift from occurring within regions and therefore maintaining the genetic homogeneity in the latent population. We cannot, therefore, dismiss the possibility that outbreaks arise almost simultaneously and/or independently across the entire area of the initiation box, as suggested by Fabricius et al. [13].

While our genetic analyses did not resolve competing hypotheses about the initiation and spread of CoTS outbreaks on the GBR, these data could be used in conjunction with demographic information and/or fine-scale monitoring data to better resolve the patterns of initiation and spread. Extensive and intensive monitoring of CoTS populations was undertaken across the northern GBR throughout the period of this study, recording the extent and severity of outbreaks at every reef between Lizard Island and Cairns (P. Doherty, unpublished data), as well as documenting the size-structure of CoTS populations at select reefs throughout this range [67]. However, these surveys were undertaken (in 2014 to 2015) only after outbreaks had become well established throughout the entire area, such that sequential (cf. simultaneous) occurrence of outbreaks will only be apparent based on spatial variation in the size and abundance of CoTS. In the future, systematic and intensive monitoring should be undertaken across a range of reefs within the initiation box to unequivocally establish the sequence and inter-dependence of outbreaks within this area [2]. It is also possible (but not certain) that next generation sequencing might reveal greater genetic structure among existing samples and thereby provide meaningful differences among sub-populations to explicitly test for directionality in spread. The recent compilation and publication of an entire mitochondrial genome for Acanthaster cf. solaris collected from Japan [68] certainly paves the way for much more detailed studies of population genetics for CoTS. Ongoing genetic sampling and re-analyses of existing genetic samples from the GBR are underway.

An unexpected outcome of this study was that outbreak populations of CoTS in the southern GBR (Swains reefs) were not significantly differentiated and have a similar origin to outbreak populations sampled in the northern and central GBR during 2014/2015. For the most part, CoTS outbreaks in the Swains have been thought to occur independently of outbreaks in the northern GBR and have an altogether different origin (e.g., [19]), though the appearance of high CoTS densities at Swains reefs in 2014 is consistent with continual and progressive southerly spread of the outbreak that started at and near Lizard Island in 1993/1994 [69]. Model-based Bayesian analyses resulted in a higher posterior probability of selected scenarios that considered a common origin for outbreaks in both the Swains reefs and northern GBR, as opposed to two independent primary outbreaks (sensu [19]). Both scenarios that represent a single common primary outbreak indicate a strong goodness-of-fit to our genetic data, though we could not distinguish whether secondary outbreaks originated from a single time point or sequentially. It is possible that the limited time (in generations) elapsed between primary and secondary outbreaks, and even between successive waves of outbreaks, would result in minimal observable genetic differences.

5. Conclusions

Resolving the specific location where reef-wide outbreaks of CoTS actually originate on the GBR remains a high priority, both to understand potential causes or triggers of outbreaks and to develop more effective and highly targeted management responses [2]. However, our capacity to establish the origin of the current outbreak was significantly constrained by the limited genetic structure apparent based on 17 microsatellite loci. Very high homogeneity observed within the current outbreak population, with limited evidence of inbreeding, suggests that rapid expansion in population size most likely arose from multiple and undifferentiated latent populations. Indeed, our data suggest that CoTS outbreaks may have occurred almost simultaneously and independently across the entire area of the ‘initiation box’, from Cairns to Lizard Island. The priority, therefore, is to undertake intensive sampling on reefs throughout the ‘initiation box’ in the lead up to the next outbreak of CoTS on the GBR. However, ongoing sampling should be combined with testing of alternative molecular markers (e.g., SNPs) and the application of new techniques to sample larvae [70] and newly settled starfish [71] to further resolve the origin of reef wide outbreaks.