Population genetics and demography of the coral-killing cyanobacteriosponge, Terpios hoshinota, in the Indo-West Pacific

Sampling

All specimens were collected with permission and sampling permits where required (permits not required for sampling in Japan, Taiwan, China). A total of 234 specimens of Terpios hoshinota were examined in this study (Fig. 1, Table 1). We collected specimens at each location ranging from a minimum of 1 to a maximum of 41 specimens. This variation and non-uniformity in sampling was due to several factors: (1) variation in numbers of T. hoshinota at different sites, (2) field workers not identifying T. hoshinota in the field, and (3) sampling restrictions in some cases (e.g., the Great Barrier Reef). After collection, macro-invertebrates and commensals were carefully removed. All samples were preserved in 75% ethanol inside centrifuge tubes and stored at −20 °C until extraction.

For Guam, samples were collected with permission of the Department of Agriculture Division of Aquatic and Wildlife Resources (DAWR) Guam and certificate of origin from DAWR was obtained prior to shipping. Permit is not required to collect sponges in Green Island and in China. For Indonesia, permit was issued by the Indonesian Ministry of Research, Technology & Higher Education (KLN RISTEKDIKTI scheme, 2016-2018) No. 011/SP2H/LT/DRPM/IV/2017. For Malaysia, sampling permit was issued by the Marine Park and Resource Management Division, Department of Fisheries Malaysia (formerly Department of Marine Park Malaysia) with a Sampling permit: JTLM 630-7 JLD. 4(8). Samples from Australia were part of Census of Marine Life/Census of Coral Reefs, with permit No. G32313.1 issued by the Great Barrier Reef Marine Park Authority.

DNA extraction and PCR amplification

Genomic DNA was extracted from 25 mg of clean sponge tissue using the DNAeasy-Tissue Kit (Qiagen) according to the manufacturer’s instructions. Specimens used for COI sequencing were from Green Island (n = 12), Lanyu (3), Wanlitung (1), Guam (2), Maldives (5), Miyako (2), Bise (6), Shiraho (1), Yakomo (2), Lizard Island (1), and Xisha (3). A fragment containing 650-bp COI of T. hoshinota was amplified by the universal primer set (Folmer et al., 1994) LCO 1490: 5′-GGTCAACAAATCATAAAGATATTGG-3′) and HCO 2198: 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′). The thermal cycles for PCR consisted of 1 cycle at 94 °C for 2 min, followed by 35 cycles at 94 °C for 50 s, 40 °C for 55 s, 72 °C for 60 s, and a final cycle of 72 °C for 7 min. The amplified products were directly sequenced in both directions (Mission Biotech Taiwan, http://www.missionbio.com.tw).

In addition, universal primers for amplifying the entire nuclear ribosomal internal transcribed spacer (ITS), (ITS4: 5′-TCCTCCGCTTATTGATATGC-3′; EF3RCNL: 5′-CAAACTTGGTCATTTAGAGGA-3′) (Lord et al., 2002), were used to design specific ITS 2 primers for T. hoshinota as below. The number of samples used for cloning were; Green Island (n = 9) LiuMaGou (2), ChaiKou (1), BiTouXi (2), ZhaoRi hot spring (2), YouZi Lake (1) and Shiraho, Japan (1). DNA samples were amplified by PCR with 1 cycle at 94 °C for 3 min, 94 °C for 30 s, 54 °C for 90 s, and 72 °C for 60 s at 35 cycles, followed by a final cycle of 72 °C for 10 min. The length of the PCR products (∼650 bp) was then confirmed by gel electrophoresis (1% agarose). PCR products were then cloned by using an Invitrogen TA cloning kit (Thermo Fisher Scientific, https://www.thermofisher.com/tw/zt/home/life-science/cloning/ta-cloning-kits.html) following the manufacturers protocol. Treated plasmids containing Ampicillin (50 µL), IPTG (0.5 mM), X-Gal (50 µg/ml) were then added into Luria Bertani agar plates for cloning. Positive clones were selected during the transformation step and a minimum of 20 clones were picked up for each set of sponge samples, followed by PCR with universal primer (ITS4, EF3RCNL). A total of 111 positive clones were selected and amplified for the completed ITS region. Cloned ITS DNA sequences were submitted to GenBank using the blastn option, and sequences that matched 100% to non-poriferan species were excluded from the following analyses. Of the 111 ITS clones, 72 sequences matched sponge species (Hymeniacidon heliophila, Accession: AB250764) with 80 to 89% similarity values (Fig. S2). Primers for highly variant regions of internal transcribed spacer 2 (ITS2) were designed by Geneious Pro 4.5 “Primer3’ (Untergasser et al., 2012) module with default value and rectified its structure to prevent hairpin structures or self-annealing problem by Primer Select (Lasergene). The designed primer set was ITS2Fwd: 5′- TTTGCAGACGGACAGCCTCA -3′, ITS2Rev: 5′- TTTTACTGTGCACCCCTCTCAG-3′). Temperature Gradient PCR test was performed to find the best annealing temperature for ITS2 Primer, temperature starting from 47 °C to 65 °C were set as temperature gradient scope coming by 10 steps of temperature gradient with each gradient steps up at 2 °C. Gel electrophoresis images showed that the best annealing temperature was 53 °C. Optimum PCR condition of ITS2 primers are: 1 Cycle at 94 °C for 3 min, 94 °C for 30 s with 53 °C for 30 s and 72 °C for 30 s at 35 cycles, then finished by 1 cycle of 72 °C for 10 min. These new primers were then used for subsequent PCR amplification and direct sequencing of the remaining T. hoshinota specimens (Genomics Taiwan, https://en.genomics.com.tw).

Data analyses

Due to unequal sampling numbers across different reefs, we collapsed specimen groupings (Table 1, Tables S2, S3) based on marine ecoregions (Spalding et al., 2007) to obtain better resolution on the genetic diversity of T. hoshinota in the IWP. DNA sequences were cleaned visually based on quality and accuracy using Geneious R10 (https://www.geneious.com). Sequences were trimmed and aligned by the MAFFT algorithm (Katoh & Standley, 2013) for downstream analyses. All COI sequences obtained in this study were submitted to GenBank with accession numbers MZ914515 –MZ914532 and OK576574 –OK576593. Previously reported COI haplotypes available from GenBank were downloaded and added to our dataset (KJ008098, MN507872, MN507873, KP764915, KP764916, MN507874 –MN507878). Genetic diversity indices nucleotide diversity (π), number of haplotypes (h), and haplotype diversity (h_d) were calculated using DnaSP V.6.12.03 (Rozas et al., 2017). Mean uncorrected pairwise genetic distances (p- distance) were calculated with MEGA 7 (Kumar, Stecher & Tamura, 2015). Relationships among sequenced samples were analyzed using a statistical parsimony network algorithm implemented in PopART V.1.7 (Leigh & Bryant, 2015).

ITS2 chromatograms were analyzed following the procedure for nuclear sequence markers described in (Fontaneto, Flot & Tang, 2015; Fontaneto, Flot & Tang, 2015). Sequences were phased using the statistical website phylogeny.fr. Phasing was trivial for homozygous sponges (that had no double peaks) and sponges that had only one double peak (i.e., comprising two ITS1 sequence types that differed at a single position). Sponges that had only a few double peaks, as expected from the superposition of two ITS2 sequence types of equal lengths differing at several positions, were phased using the programs PHASE V.2.1.1 (Stephens, Smith & Donnelly, 2001) and SeqPHASE (Flot, 2010) where the allele combinations with probability scores above 0.8 were chosen for the analyses. Using the ITS2 marker, 234 samples of T. hoshinota were successfully sequenced. There were 191 homozygous and 43 heterozygous individuals, and hence a total of 277 sequences corresponding to 23 haplotypes (128 bp long) was obtained (Table 1). All sequences were submitted to GenBank with accession numbers MZ468904 –MZ469137. Haplotype networks were constructed with a statistical parsimony network algorithm (Templeton, Crandall & Sing, 1992) by PopART V.1.7 software based on marine ecoregions (Spalding et al., 2007). Genetic diversity indices mean uncorrected pairwise genetic distances (p- distance) were calculated as described in the previous section (COI sequence analysis). Best-fit models for nucleotide substitution were performed by jModelTest V.2.1.10 (Darriba et al., 2020), with the most appropriate evolutionary model TrN for the nucleotide substitution model of the ITS2 sequences. To test past demographic events, we used Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) in DnaSP with 10,000 simulations under the selective neutrality model. Mismatch distribution graphs were plotted in DnaSP V.6.12.03 to reveal evidence of spatial range expansion or stable population (Tajima, 1989). Demographic changes were also examined by the Ramos-Onsins’ R₂ in DnaSP V.6.12.03 with 10,000 permutations, which is more powerful for detecting past demographic events than Harpending’s raggedness index (Harpending et al., 1993) or Fu’s Fs when sample size is relatively small (Ramos-Onsins & Rozas, 2002). R₂ can detect recent strong population growth, and a non-significant value shows no deviation to the null hypothesis and thus represents a constant population size (Ramos-Onsins & Rozas, 2002). Population differentiation was examined using the Weir and Cockerham F_ST estimator (Weir & Cockerham, 1984) by Arlequin 3.5 (Excoffier & Lischer, 2010), and significance level (p = 0.001) was assessed using 10,000 permutation replicates. Global F_ST and average F_ST per population following with UPGMA clustering of population were obtained from GraphPad Prism 9. Isolation-by-distance analysis (Slatkin, 1987; Slatkin, 1993) was performed through linearized pairwise F_ST estimates (F_ST /1- F_ST), and correlated geographical distances of T. hoshinota populations based on marine ecoregions using the standard Mantel test implemented in GENODIVE V.3.05 (Meirmans, 2020). The linear regression of linearized pairwise F_ST values and log transformed geographical distances were plotted by GraphPad Prism 9.

Mitochondrial COI genetic diversity and statistical parsimony network

The COI sequences obtained from 37 Terpios hoshinota specimens across 10 populations in the IWP were nearly identical, with only 3 phylogenetically informative sites identified. By adding 10 sequences from GenBank, the COI sequences for T. hoshinota formed a single statistical parsimony network comprising five haplotypes, with 1 to 2 mutational steps evident among sequences (Fig. 1A) having a mean nucleotide diversity (π) of 0.00185, haplotypic diversity (h_d) of 0.416, and p-distance of 1.09 (Table S3), respectively. C2 was the dominant haplotype (n = 35, 74.47%) in specimens from the major marine ecoregions, except for specimens from the Torres Strait (Northern GBR). C3 was ranked as the second-most dominant haplotype (n = 9, 14.89%) in samples from Guam, Indonesia, and southern Kuroshio. The remaining three COI haplotypes were each identified from single individual specimens from Xisha, the GBR, and the Maldives. Due to the small sample sizes and relatively low genetic variability, COI was not utilized in the historical demography and population genetic analyses.

ITS2 genetic diversity, demography, population differentiation, and statistical parsimony network

The number of samples (n), nucleotide diversity (π), p-distance, haplotype (h), and haplotype diversity (h_d) of ITS2 in Terpios hoshinota based on marine ecoregions are shown in Table 1. Identical ITS2 sequences in two samples from the Torres Strait, northern GBR, were excluded from the comparisons of π, h_d, p- distance, and following demography and F_ST -statistics. The mean p-distances ranged from 0.269 in the Maldives to 3.132 in the South China Sea Oceanic Islands. Similar patterns were seen in nucleotide diversity, with 0.0021 in the Maldives and 0.00245 in the South China Sea Oceanic Islands. The number of haplotypes ranged from 1 in the Torres Strait, northern GBR, to 12 in South Kuroshio, with haplotype diversity ranging from 0.117 in the Maldives to 1 in Northeast Sulawesi.

The historical demography of Terpios hoshinota inferred by Tajima’s D was positive for most populations except those from Sunda Shelf, Sulawesi Sea, and the Maldives (Table 1). The results of Fu’s Fs test also showed positive values for most populations except those from South Kuroshio, Southern Java, Sulawesi Sea, northeast Sulawesi, and the Maldives. Only the Sulawesi Sea and Maldives had both negative and statistically significant Tajima’s D and Fu’s Fs, indicating an excess of rare haplotypes over what would be expected under neutrality and rejecting the null hypothesis of constant population size at mutation-drift equilibrium; i.e., Sulawesi Sea and Maldives populations might be experiencing population expansions or be under positive selection. The remaining populations were not statistically diverged from neutrality, suggesting their population sizes were at mutation-drift equilibrium and none significantly deviated from neutrality (Ramos-Onsins’ R₂, Table 1). Four types of mismatch distribution models were observed in Terpios hoshinota across the IWP: unimodal, skewed unimodal, multimodal, and bimodal (Table 1, Fig. S1). While Lesser Sunda showed a unimodal distribution that suggested an association with a sudden expansion model, both the Sulawesi Sea and the Maldives showed skewed unimodal distributions that suggested recent expansions or bottleneck events. The remaining populations had either multimodal or bimodal distributions that suggested constant population sizes or the presence of secondary contacts of two distinct lineages (Alvarado Bremer et al., 2005; Jenkins, Castilho & Stevens, 2018).

Two hundred and seventy-seven ITS2 sequences of T. hoshinota were categorized into 23 haplotypes across the IWP (Fig. 2), forming a single statistical parsimony network comprising 1 to 5 mutational steps evident among specimens (Fig. 1B). H3 was the most dominant haplotype (n = 133, 48.01%) present in all sampling sites except the Maldives, GBR, Xisha, and Yakomo (Fig. 2). H6 was the second-most dominant haplotype (n = 35, 12.64%), and present only in sites near the equator or in the southern hemisphere, including the Sulawesi Sea, Southern Java, the Maldives, and GBR. H2 was the third-most dominant haplotype (n = 31, 11.91%), present in samples from Guam, South China Sea Oceanic Islands, Sunda Shelf, and Southern Java. The rest of the haplotypes were present either in low frequencies or restricted to single populations; for example, H22 and H23 occurred in Xisha, China, H19 in southern Java, and H18 in Papua, Indonesia (Figs. 1B, 2).

Population differentiation based on marine ecoregions was examined by pairwise F_ST -statistics (Fig. 3, Table S4). Eighteen of 55 (32.73%) pairwise F_ST values were significantly different (p <0.001), with negative F_ST values in South Kuroshio - Sulawesi Sea (−0.003), South Java - Northeast Sulawesi (−0.1017), and Maldives –GBR (−0.286). Negative F_ST values indicate more variation within than between populations or results from uneven sample sizes (Holsinger & Weir, 2009; Gerlach et al., 2010). The F_ST heatmap (Fig. 3) clearly showed that the Maldives, located in the central Indian Ocean, had significant pairwise F_ST values (except one with Northeast Sulawesi) in comparison with all other populations in the marine ecoregions of the eastern Indian Ocean, Coral Triangle, and West Pacific, with the largest F_ST (0.6687) shown between the Maldives and Sunda Shelf. Testing for a pattern of isolation-by-distance by performing a Mantel test between genetic (F_ST) and spatial distances revealed a significant relationship (r² = 0.2144, p < 0.05) for T. hoshinota across the IWP (Fig. 4).

Variability of COI and ITS2 in Terpios hoshinota

COI sequences obtained from this study showed extremely low intraspecific variation in T. hoshinota. A slow evolution pattern of mitochondrial genomes has been reported in lower metazoans, including certain lineages of poriferans and anthozoans in general (Shearer et al., 2002; Tseng, Wallace & Chen, 2005; Wörheide, 2006; Huang et al., 2008; Chen, Kuo & Chen, 2009; Pöppe et al., 2010; Sperling et al., 2012). In the 650-bp COI fragment, only 5 haplotypes with a low p- distance of 1.09 and low nucleotide diversity (π = 0.00185) were identified from 47 T. hoshinota specimens collected across the IWP. Slow mitochondrial gene evolution was previously also found in the coralline demosponge Astrosclera willeyana across the Indo-Pacific, which contained only three COI haplotypes with a maximum p-distance of 0.418 and low nucleotide diversity (π = 0.00049) (Wörheide, 2006). As well, in yellow candle sponges (Aplysina spp.) from the Caribbean, there was no interspecific variation in the COI fragment and only 6 nucleotide differences among species observed in the complete mitochondrial genome (19,620-bp) (Sperling et al., 2012). Conversely, in recent studies on some deep-sea sponges such as Poecillastra laminaris, COI shows significant genetic differences among provinces and geomorphic features in New Zealand (Zeng et al., 2019).

In contrast to COI, ITS2 sequences, although only 128 bp in length, had much higher genetic variation (mean p-distance of 1.717 and π of 0.0164) and presented 23 haplotypes across the IWP region, suggesting that ITS2 could provide more information for intraspecific inferences in T. hoshinota. However, heterogeneous levels of ITS intra-genomic polymorphisms (IGPs) in a diverse range of marine sponge taxa have been documented, and it has been recommended that caution must be taken to survey the different levels of IGPs before applying ITS in phylogenetic studies below the species level (Wörheide, 2006; Wörheide, Epp & Macis, 2008). In this study, we developed ITS2 through a series of sub-cloning evaluations that enabled us to perform PCR and directly sequence the entire region. Relatively low IGPs were detected in ITS2. Following the procedure described in Fontaneto, Flot & Tang (2015), ITS2 could be used as a potential marker to infer the historical demography and population genetics of T. hoshinota across the IWP.

Potential dispersion and outbreak mechanisms: long-range dispersal, local outbreaks, or both?

COI and ITS2 produced similar statistical parsimony networks (SPN) showing a dominant haplotype comprised of several mutational steps evident among sequences from specimens. Although the number of specimens and sampling locations were not the same for the two markers, our results indicated that, by count, ITS2 (π = 0.0164) evolved 10 times faster than COI (π = 0.00185). Both markers resolved most of the representative haplotypes from the distributional range of Terpios hoshinota in the IWP, finding 5 COI haplotypes out of 47 samples (10%) and 23 ITS2 haplotypes out of 234 samples (9%). In the COI SPN, haplotype C2 was dominant, corresponding well with the most dominant ITS2 haplotype H3 from Southern Kuroshio, H2 from Mariana Islands, H6 from Maldives, and H22 and H23 from Xisha, South China Sea Islands. Interestingly, C2 also corresponded to the COI haplotype previously identified from the Spermonde Archipelago, Indonesia (vander Ent, Hoeksema & de Voogd, 2016), and the Kimberley inshore reef, Western Australia (Fromont, Richards & Wilson, 2019), suggesting that C2/H3 is the major haplotype in shaping the demography, population genetic structure, and probably outbreaks of T. hoshinota in the IWP.

Since its first discovery in Guam (Bryan, 1973), occurrences and outbreaks of Terpios hoshinota have been reported along the South Kuroshio marine ecoregion towards the South China Sea Oceanic Islands, Torres Strait Northern GBR, marine ecoregions in the Coral Triangle including those of Malaysia and Indonesia, and more recently into the Indian Ocean (Table 1, Table S2). By superimposing the main ocean currents onto the sampling sites and ecoregions (Fig. 2), long-range dispersal, or the “recent expansion, perhaps from Guam” hypothesis could, at least in part, explain how T. hoshinota might have dispersed in the IWP. In this hypothesis, T. hoshinota is represented by its two most “common” ITS2 haplotypes, H3 and H2, which spread from Guam along the North Equatorial Current (NEC) to the west and the Kuroshio Current to the north to invade reefs in Japan and Taiwan (Liao et al., 2007; Chen, Kuo & Chen, 2009; reviewed in Reimer et al., 2011). These two ITS2 haplotypes also invaded the South China Sea Oceanic Islands along a branch of the Kuroshio Current, entering the South China Sea through the Luzon Strait and then traveling clockwise on the South China Sea Surface Current (SCSSC) to cause outbreaks in Itu Aba, Spratlys Yang et al., 2018, and the Sunda Shelf in the Gulf of Thailand (Fig. 2). On the other hand, the other branch of the NEC, the Mindanao Current, turns southward and transported T. hoshinota to the reefs around the Coral Triangle. The Mindanao Current also accounts for the trans-equatorial transportation of T. hoshinota from the Coral Triangle to invade Kimberley inshore reefs, Western Australia (sharing the C2/H3 haplotype) (Fromont, Richards & Wilson, 2019). Trans-equatorial transportation of T. hoshinota is further supported by long-term benthic surveys and examination of museum collections, which indicate that T. hoshinota may have arrived recently in the Kimberley, possibly after 2014 (Fromont, Richards & Wilson, 2019). Similar observations were also seen with regards to the T. hoshinota outbreaks at Green Island reefs that appeared in 2007 but were not recorded in ReefCheck surveys between 1999 and 2006 (Liao et al., 2007; Chen, Kuo & Chen, 2009).

Results from pairwise F_ST and isolation-by-distance (IBD) is consistent with the “recent expansion, perhaps from Guam” hypothesis. Only 32.73% of pairwise F_ST values were significant, suggesting that there were recent expansion events or moderate gene flow existed among marine ecoregions in the West Pacific, Coral Triangle, and eastern Indian Ocean, but not the central Indian Ocean. The F_ST heatmap clearly showed that the Maldives had much higher F_ST values compared to counterparts in other marine ecoregions (see Fig. 3 and Table S4). In addition, IBD, a pattern of population differentiation based on the stepping-stone model of population structure in which genetic differences between populations increase with geographic scale, was also significant. Trans-oceanic genetic divergence across the Pacific and Indian Oceans has been documented in many tropical marine taxa (Crandall et al., 2008; Wörheide, Epp & Macis, 2008; Bowen et al., 2013; Ma et al., 2018). A phylogeographic study on the lemon-yellow calcareous sponge, Leucetta chagosensis, showed that while populations from the Philippines, the Red Sea, the Maldives, Japan, Samoa, and Polynesia were reciprocally monophyletic, suggesting long-standing isolation, the populations along the South Equatorial Current in the south-western Pacific, on the contrary, showed isolation-by-distance effects. It was concluded that the dispersal pattern of L. chagosensis was stepping-stone dispersal with some putative long-distance exchanges (Wörheide, Epp & Macis, 2008).

Despite the possibility of ocean currents mediating the dispersal of Teripos hoshinota in the IWP, the alternative hypothesis of “local outbreak” should also be considered, wherein the sponge may appear in a new environment through occasional or rare events. Once T. hoshinota arrives on a new reef, it might remain undetected at low abundances for a long period of time until an environmental disturbance triggers a new outbreak. Larval dispersal of T. hoshinota is short, and its pseudoblastula larvae have been shown to crawl out from the oscula and settle rapidly within 1 day, subsequently initiating metamorphosis on substrates consisting of dead coral skeletons and a glass petri dish but not on living corals (Wang et al., 2012a; Wang et al., 2012b; Hsu, Wang & Chen, 2013; Nozawa, Huang & Hirose, 2016). The absence of planktonic larvae and having quick settlement runs against basic assumptions of the long-range dispersal hypothesis. However, rafting with floating materials (e.g., pumice), as seen in other benthic organisms, is consistent with long-range dispersal (Jokiel, 1984; Jokiel, 1990; Bryan et al., 2012), and rapid movement caused by typhoons (Nozawa, Huang & Hirose, 2016) could be another alternative way for T. hoshinota to travel long distances with ocean currents.

Moderate gene flow with IBD and a demographic history of mutation-drift equilibrium without significant expansion supports this hypothesis with respect to T. hoshinota dispersal and outbreaks in the IWP. All the indices (Tajima’s D, Fu’s F_S statistics, mismatch modelling) used to examine the demographic history of T. hoshinota based on marine ecoregions showed that most populations were consistent with the model of drift–mutation equilibrium with no evidence of selection, except for populations in the Sulawesi Sea and the Maldives (significantly negative Tajima D, Fu’s Fs), and skewed unimodal mismatch distribution may suggest sudden population expansions with bottlenecks in these populations. However, Ramos-Onsins’ R₂, which is more powerful for detecting past demographic events, especially in populations with small sample sizes, remained non-significant in all populations, suggesting that T. hoshinota might have existed in the sampled reefs for a long period of time in relatively low abundances, followed by local outbreaks through asexual propagation and also consequently noted during reef surveys. Previous studies have shown that T. hoshinota can increase its coverage through asexual propagation at 1.5–2 times the speed of other sponges ((Rützler & Muzik, 1993), Rossi et al., 2015), and spread twice as fast as scleractinian corals (Elliott et al., 2016a; Elliott et al., 2016b). Therefore, local outbreaks such as those in Japan and Taiwan might have resulted from asexual propagation of the dominant haplotypes within a short period of time, leading to low genetic diversity without any significance in demographic indices.

Limitation of conventional molecular markers

Contrasting interpretation showing a moderate gene flow detected among the marine ecoregions with a clear pattern of isolation-by-distance (IBD) and demography analyses hinting that T. hoshinota populations were in the status of mutation-drift equilibrium without significant expansion might be due to the characteristics of molecular markers. While mtCOI is slow-evolving, ITS2 is too short. Also, both contain few phylogenetic informative sites that provided limited information for genetic analyses. Terpios hoshinota is a cyanobacteriosponge containing relatively high density of symbiotic cyanobacteria (Hirose & Murakami, 2011; Tang et al., 2011) that retards the development of Terpio s-specific gene loci. Future applications of the single-cell capture method (Tang et al., 2009) and next generation sequencing are expected to obtain in-depth genomic information for genetic research in T. hoshinota.

Population genetic implications for anthropogenic and natural disturbances causing T. hoshinota outbreaks

Occurrences of Terpios hoshinota have been reported from the West Pacific from the early 1970s (Bryan, 1973; Plucer-Rosario, 1987), and now extend into the Indian Ocean (Elliott et al., 2016a; Elliott et al., 2016b; Montano et al., 2014; Fromont, Richards & Wilson, 2019). Nevertheless, the ecological mechanisms of T. hoshinota outbreaks remain uncertain. The outbreaks are suspected to be related to large-scale disturbances caused by coastal development (Rützler & Muzik, 1993; Chen, Kuo & Chen, 2009) or volcanic eruptions (Schils, 2012), which could have resulted in increased turbidity, inshore nutrient run-off, and trace metals such as iron from volcanic ash being present in seawater. These inputs accelerate the propagation of symbiotic cyanobacteria and growth of sponge hosts, resulting in the spread of T. hoshinota. Once a population outbreak is established, it can occupy reefs for decades, only disappearing possibly due to strong waves caused by storms, but still having the ability to revert to outbreak proportions at a later time (Reimer et al., 2011; Schils, 2012; Yomogida et al., 2017). As such, occurrence and outbreak through “long-distance dispersal” and “local outbreak” are non-mutually exclusive. Meaning both of these two mechanisms may have resulted in the outbreaks of T. hoshinota in the past 20 years. It is important to continue long-term ecological studies using newly developed genomic applications to unveil the outbreak mechanisms of Terpios hoshinota, particularly given the increasing impacts of climate change.