Integrating phylogeographic and ecological niche approaches to delimitating cryptic lineages in the blue–green damselfish (Chromis viridis)

Sample collection and DNA extraction

We collected 252 C. viridis between 2007 and 2017 from 15 locations across its Indo-Pacific distribution (Table 1; Fig. 1). We collected specimens by hand net and clove oil, either by scuba diving or snorkeling and preserved tissue samples (fin clips, a piece of muscle, or both) in 95% alcohol and stored at 4 °C. Large fish were released following fin clipping, but individuals too small for fin clipping were euthanized with clove oil and preserved whole in 95% ethanol. The sampling depth of all specimens used in present study was less than 10 m. This experiment doesn’t involve animal experiment and the field sampling process complied with the regulation drafted by the Animal Care and Use Committee of National Sun Yat-sen University. Additionally, the permission to perform research activities in Indonesia was issued by the Indonesian Government and the Ministry of Research and Technology under research permit No. 272/SIP/FRP/SM/VII/2013.

We extracted genomic DNA from tissue fragments using Geneaid Tissue Genomic DNA mini Kit (Geneaid Biotech, New Taipei City, Taiwan) following manufacturers’ protocol and eluted extracted DNA in TE buffer and stored at −20 °C until amplification by polymerase chain reaction (PCR).

Amplification of genetic markers

We amplified a portion of the mitochondrial cytochrome b (Cytb) gene using universal primers GluDG-L and H16460 (Palumbi, 1996). PCRs reactions were 30 μL in volume, containing 10–40 ng template DNA, three μL 10X buffer, 0.2 mM dNTPs, 1.5 mM MgCl₂, 10 mM each primer, and 0.2 units of Taq polymerase (MDbio, Taipei). The thermocycling profile consisted of initial denaturation at 94 °C for 2 min, followed by 41 cycles of denaturation at 94 °C for 30 s, annealing at 57 °C for 30 s, and extension at 72 °C for 40 s, concluding with a final extension at 72 °C for 2 min.

Because mitochondrial and nuclear genes can have very different histories, we also amplified the nuclear recombination-activating (Rag2) gene for a subset of samples that represented distinct clades recovered in the mtDNA, using primers RAG2F and RAG2R (Westneat & Alfaro, 2005). PCR reactions were as above, except that we used a five mM MgCl concentration and the following thermocycling parameters: 39 cycles of denaturation at 94 °C for 30 s, annealing at 56 °C for 30 s, and extension at 72 °C for 40 s, and a final extension at 72 °C for 2 min.

The nucleotide sequences of the PCR products of both loci were determined using ABI 3730XL automated sequencer (Carlsbad, CA, USA) by Genomics (https://www.genomics.com.tw). Nucleotide sequences were assembled and edited using the SEQUENCHER version 4.2 software (Gene Code, Ann Arbor, MI, USA). Sequences were uploaded to GenBank under the accession number MH743228–MH743691.

Phylogenetic analyses

Prior to any analyses, we aligned DNA sequences from each gene region in Clustal W (Thompson, Higgins & Gibson, 1994) and exported these sequences to MEGA 6 (Tamura et al., 2013) to visually inspect all alignments for accuracy. To quantify genetic diversity measures, we calculated standard genetic diversity indexes including haplotype diversity (h) and nucleotide diversity (π) in Arlequin 3.5 (Excoffier & Lischer, 2010). We then inferred phylogenetic trees from each individual locus using maximum likelihood (ML) and Bayesian inference performed on the CIPRES Science Gateway (Miller, Pfeiffer & Schwartz, 2010), and in BEAST 2.4.5 (Bouckaert et al., 2014). We conducted the ML analyses in RAxML version 8.1.24 (Stamatakis, 2014) using the GTR+G substitution model which selected by MEGA 6 (Tamura et al., 2013) as best-fit substitution model. For Bayesian inference, we used MrBayes version 3.2.2 (Ronquist et al., 2012), implementing two parallel runs of four simultaneous Markov chains for 10 million generations, sampling every 1,000 generations and using the default parameters. Run parameters employed unlinked substitution models, rated heterogeneity models, and based frequencies across partitions. In the Bayesian analyses the first one million generations (10%) were discarded as burn-in, based on log-likelihood tree scores. Meanwhile, the convergence diagnostic was applied and the stop probability was set to 0.01. Nodal support was evaluated individually for all trees using non-parametric bootstrapping with 1,000 ML replicates as employed in RAxML (Stamatakis, 2014), and by calculation of posterior probabilities as employed in MrBayes. Lastly, we generated a median-joining haplotype networks based on Cytb and Rag2 sequence datasets by using Popart 1.7 (http://popart.otago.ac.nz).

To date the ages of C. viridis clades, we arbitrarily subsampled 2 Cytb sequences from each clade; we also included two outgroups, including C. atrilobata (AY208524.1) and C. multilineata (AY208533.1) to calibrate date ranges as these two taxa diverged 3.1 Mya, (Quenouille, Bermingham & Planes, 2004). We created the XML BEAST input file using the software BEAUti v.1.8.2. (as implemented in BEAST), with the setting of 50 million generations under the uncorrelated relaxed clock model, and sampling trees once every 1,000 generations. To test for inter-run variation, we conducted two independent runs in BEAST 2.4.5 (Bouckaert et al., 2014), and then checked these runs for convergence with the software Tracer v1.5 (available at http://beast.bio.ed.ac.uk/Tracer/). After discarding the 20% burn-in, we pooled the remaining tree samples from the two runs into a combined file to calculate the maximum clade credibility tree. From this tree, we calculated the posterior mean divergence ages and 95% credibility intervals for all nodes using Tree Annotator v1.8.2 (implemented in the BEAST package). To compare the similarity of two gene tree topologies, Phylo.io software (Robinson, Dylus & Dessimoz, 2016) was used.

Ecological niche characterization and comparison

To characterize the ecological niches of C. viridis (C. viridis Clade A) and potential cryptic species (C. viridis Clade B), we obtained geophysical, biotic, and environmental data layers for global sea surface from Bio-ORACLE (downloaded on February 20th, 2018; Tyberghein et al., 2012; Assis et al., 2018). We then extracted the values of those factors for the locations where C. viridis samples were collected. We excluded factors such as ice thickness and sea ice concentration from our analyses because all sample sites were tropical and ice free.

To explore niche differentiation between C. viridis Clade A and C. viridis Clade B, we first compared the values of individual environmental factors between C. viridis Clade A and Clade B sampling sites using Mann–Whitney tests. Next, we plotted niche differences in a two-dimensional non-metric multidimensional scaling (NMDS) space based on the Euclidean dissimilarity of those factors. Due to the different units of the environmental factors, we standardized factor values before running the NMDS analysis. Lastly, we tested equivalency of ecological niches between the species based on environmental factors of C. viridis Clade A and Clades B sampling locations in a principle components analysis (PCA) framework. Due to potential correlations among environmental factors, we ran the PCA on all environmental factors and then used the resulting principal components to define niches of the species. Because C. viridis mainly occur in coral reefs, we restricted our PCA to regions (grid cells) where the maximum depth was equal or less than 50 m based on the bathymetry data obtained from MARSPEC (http://www.marspec.org/; Sbrocco & Barber, 2013). To make the spatial resolution of the bathymetry data layer, which is originally at 30-by-30-arc second resolution, consistent with the Bio-ORACLE data, we upscaled it to 5-by-5-arc min resolution with each pixel value in the upscaled data layer being the maximum value over 10 × 10 original pixels within that new pixel.

To measure niche similarity between C. viridis Clade A Clade B, we calculated the D and I statistics developed by Warren, Glor & Turelli (2008). These two indices range from zero and one with higher values indicating more similar niches. To test for niche equivalency between these two cryptic lineages of C. viridis, we compared the observed similarity values to null distributions of similarity obtained through a randomization procedure under a hypothesis of an identical niche (Warren, Glor & Turelli, 2008), using the “dismo” package in R (version 3.4.2).

Niche divergence tests

Results spanning 11 environmental factors (Table 2) showed no significant difference in ecological niche parameters of Clade A or Clade B of C. viridis. Moreover, there was also no significant environmental difference between the locations where two clades co-occur and those where only Clade B occurs (Table 2), suggesting that the environmental factors cannot explain why Clade B co-occurs with Clade A at some locations, but not at others.

The NMDS analysis also showed a lack of niche differentiation. The first two NMDS axes explained 86.7% of the variation in the 11 environmental factors across the 15 sampling sites. The three sampling sites where C. viridis B occurs fall at the center of the convex hull of all sampling sites in the two-dimensional NMDS space (Fig. 4), indicating complete niche overlap between Clades A and B. The first four principal components of the environmental factors, which were standardized before the analysis, had eigenvalues higher than 1. Together they explain about 83.1% of the total variation across all the grid cells with the maximum water depth less than or equal to 50 m. Based on the component scores corresponding to the locations where Clades A and B were found, niche similarity indices were high, with D = 0.841 and I = 0.973. Moreover, niche equivalency tests showed that none of the two values were significantly different from the values obtained through a randomization process under the hypothesis of an identical niche (P-values for D and I were 0.311 and 0.294, respectively).

Origin of lineage diversification

Although the divergence of Clades A and B of C. viridis is clear, the origin of this divergence is not. Within the Indo-Pacific region, the “Indo-Pacific Barrier” has been considered as a soft barrier that could be the main driving force of marine biological provinces in this region (Gaither et al., 2011). Barber & Bellwood (2005) and Cowman & Bellwood (2013) examined the importance of this barrier by evaluating the extent of temporal concordance in vicariance in three prominent families of reef fish, including Labridae, Pomacentridae, and Chaetodontidae. Both studies showed that the isolation effect of Indo-Pacific barrier on the widely distributed fishes mainly occurred mostly between end Miocene and Early Pliocene (2–6 Myr), and the majority of vicariance events occurred in a narrow time interval at approximately 2.5 Myr.

The date of divergence of Clade A and B broadly conforms to the onset of Plio-Pleisetocene sea level fluctuations (Voris, 2000). However, Clade B populations of C. viridis occur on islands surround by deep water where lowered sea levels would not result in land barriers that promote vicariance, unlike other marine taxa distributed on opposite sides of the Sunda Shelf (Barber, Erdmann & Palumbi, 2006; Crandall et al., 2008, 2014; DeBoer et al., 2014a, 2014b; Waldrop et al., 2016; Simmonds et al., 2018). Given the Indo-Pacific wide range of C. viridis and a divergence time between two cryptic clades dating to 3 Mya, isolation across the Indo-Pacific Barrier is the most likely driver of divergence between Clades A and B. Under this scenario, Clade A would be an Indian Ocean clade that has expanded into the Pacific, where it now overlaps with the Pacific Clade B, a process previously noted in Neritid snails (Crandall et al., 2008).

An alternate, but not mutually exclusive explanation for divergence of Clade B comes from recent studies of coral associated marine taxa. Similar to phytophagus insects that undergo ecological divergence associated with host switching (Berlocher & Feder, 2002; Hébert, Scheffer & Hawthorne, 2016), recent studies from the marine realm demonstrate that changes in coral host taxa can promote lineage divergence, potentially leading to speciation (Simmonds et al., 2018). Samples were not collected in a way that allows us to test this hypothesis, but future studies separating samples by coral host to determine whether individuals from Clades A and B exist in mixed schools in sympatric populations, and if so, whether those schools are associated with different coral hosts.

The geographic distribution of Clade B is curious in that it is observed in Komodo, but not in other populations in the Lesser Sunda Islands (e.g., Amed, Nusa Lembongan). Similarly, Clade B is observed on the Great Barrier Reef and Fiji, but not in New Caledonia, a population located between these too population. One explanation for this disjunct distribution is that Clade B individuals are relatively rare, and that greater sampling intensity would reveal Clade B haplotypes in adjacent ranges, as would be expected. Alternatively, it is possible that Clade A is gradually displacing Clade B populations, and the areas of sympatry represent locations where this process is incomplete. Similar arguments were made to explain sympatry of highly divergent clades of marine snails in the Pacific Ocean (Crandall et al., 2008).

In contrast to the divergence of Clade A and B, the divergence of Clade A-1, unique to the Red Sea, is more easily interpreted as the result of vicariance. The Red Sea is a semi-enclosed basin that is frequently invoked in driving population differentiation of reef organisms between the Red Sea and Indian Ocean (DiBattista et al., 2013). Additionally, the 95% HPD of Red Sea Clade is ranged from 0.013 to 0.332 Mya which aligns with the most recent closure of the Red Sea (Siddall et al., 2003). However, while mitochondrial DNA shows a well-supported Red Sea sub clade within C. viridis Clade A, this clade is not confirmed with nuclear Rag2 sequence.

Lack of concordance between the mtDNA and nuclear phylogenies is not surprising. First, there is substantially more genetic variation in the mtDNA sequences because DNA repair mechanisms are less efficient in mtDNA, resulting in much higher substitution rates (Alexeyev et al., 2013). Second, because mitochondrial genome is maternally inherited and haploid, its effective population size is one‐fourth that of a nuclear gene, meaning that lineage sorting occurs more rapidly in mtDNA (Hare, 2001). Thus, mitochondrial gene trees have a substantially higher probability of accurately recovering recently divergence lineages with short internode distances than do nuclear genes (Moore, 1995). Given that the Red Sea mtDNA clade only dates to 0.138 Mya, it is unsurprising that divergence in this region is not recovered in the Rag2 sequences. Thus, while the Red Sea population is clearly a sub-population of C. viridis Clade A, it is unlikely a cryptic species, like Clades A and B.

Cryptic species

Many phylogeographic studies in the Indo-Pacific find highly divergent clades, often with regions of sympatry (Barber et al., 2002; Barber, Erdmann & Palumbi, 2006; Crandall et al., 2008; Gaither et al., 2011; Liu et al., 2012; DeBoer et al., 2014a; Bowen, Karl & Pfeiler, 2007). In many cases, these cryptic clades may represent cryptic species, similar to Liu et al. (2012) where a divergent clade of Pomacentrus coelestis in Micronesia was subsequently described as a new species (Liu, Ho & Dai, 2013). The concordant phylogeographic patterns in mtDNA and nuclear DNA, and the depth of this divergence suggests that Clades A and B of C. viridis may represent two cryptic species with different geographic ranges; Clade A is a widely distributed Indo-Pacific, while Clade B is only found in Lizard Island, Komodo, and Fiji.

Comparison of Cytb and Rag2 trees revealed a limited amount of discordance. Given the relatively recent divergence, one potential explanation for this pattern is incomplete lineage sorting (Tang et al., 2012). Incomplete lineage sorting is the simplest explanation for individuals with Clade B C. viridis haplotypes in the mtDNA tree having Clade A C. viridis Rag2 sequences (Table S1). However, incomplete lineage sorting could not explain the sequence of komo4 clustered with C. viridis clade B in the Rag2 tree but clustered with C. viridis Clade A in the mtDNA tree. The later finding could likely explain by the consequence of hybridization.

Hybridization is in C. viridis is unsurprising as closely related Pomacentrids often have overlapping geographic distributions, co-occur in the same microhabitats (e.g., colonies of branching corals) (Randall & Allen, 1977). Hybridization has been observed in several sibling species of damselfish, including Abudefduf abdominalis × Abudefduf vaigiensis (Coleman et al., 2014), Amphiprion chrysopterus × Amphiprion sandaracinos (Gainsford, Van Herwerden & Jones, 2015), Amphiprion mccullochi × Amphiprion akindynos (Van Der Meer et al., 2012), Dascyllus carneus × D. marginatus (DiBattista et al., 2015), and D. reticulatus × D. aruanus (He et al., 2017). As such, both hybridization and incomplete lineage sorting are likely responsible for the discordant patterns in the two markers.

Incongruence between genetic and ecological niche divergence

Hutchinson (1978) proposed that two species cannot occupy the same ecological niche, yet niche conservatism predicts that closely related taxa retain ancestral ecological affiliations and persists in similar environments (Lord, Westoby & Leishman, 1995; Webb et al., 2002; Wiens & Graham, 2005). Allopatric sister taxa are often characterized by niche conservatism, but because geographic isolation drives speciation (Peterson, Soberón & Sánchez-Cordero, 1999; Peterson et al., 2001; Kozak & Wiens, 2006) sibling species can’t compete for the same niche space, because, by definition, they have non-overlapping geographic ranges.

The alternative to niche conservatism is niche divergence, in which sister taxa occupy different niches (Dayan & Simberloff, 2005). Niche divergence is typically associated with sympatric speciation as diversification results from reduction in gene flow associated with divergence in traits with ecological function (e.g., habitat segregation, pollinator divergence, behavioral changes, phenological shifts, and mating system shifts). Under niche divergence, young sister species with high degrees of range overlap should ecologically diverge (Dayan & Simberloff, 2005; Davies et al., 2007). While ecological modeling and phylogenetics are commonly integrated to understand the relationship between evolutionary and ecological divergence of sibling species (Kalkvik et al., 2012; Schorr et al., 2012), few studies examine the intra-specific (cryptic lineage) level (Gutiérrez‐Tapia & Palma, 2016).

Our results show that Clade A and B of C. viridis are sympatric, with Clade B haplotypes occurring only at three localities. There are two possible explanations for this pattern. First, lineage diversification could have occurred in sympatry, as seen in Indo-Pacific coralivorous snails (Simmonds et al., 2018). Alternatively, diversification could have occurred in allopatry with secondary overlap in geographic ranges as proposed for neritid snails (Crandall et al., 2008). Niche comparisons based on 11 environmental factors showed that Clade A has a broader ecological niche than Clade B does (Fig. 4). Despite this difference, there is no significant niche differentiation between these two lineages. As such, diversification is more likely to have arisen with niche conservatism rather than niche differentiation.

Given the absence of niche differentiation, our results suggest that the diversification of Clades A and B most likely occurred in allopatry, and that secondary contact is the most likely explanation for their current distribution. If true, niche conservatism could help explain the disjunct distribution of Clade B. Under this scenario, Clade A and Clade B would compete for the same niche space when they occur in sympatry. If Clade A is a superior competitor, it would gradually displace Clade B populations, potentially explaining the relative rarity of Clade B haptotypes and their disjunct distributions.

Alternatively, Pyron et al. (2015) suggested that niche conservatism could result from soft allopatry where there is low environmental heterogeneity. As noted above, the shallow Sunda and Sahul continental shelves exposed during low sea level stands (Voris, 2000), forming long land bridges that restricted larval exchange between the tropical Indian Ocean and the western Pacific (reviewed in Randall, 1998). However, because deep water passages in what is modern day Indonesia remained open, and because low sea level stands were intermittent, the Indo-Pacific Barrier has been considered a soft dispersal barrier for marine taxa (Cowman & Bellwood, 2013). The 11 environmental factors we used to define the niches of the two lineages of C. viridis only reflect broad-scale environmental conditions. However, there could be subtle variation in environmental factors, or variation could occur at finer scales (e.g., microhabitat), resulting in soft-allopatric divergence.

It is important to note that the broad environmental variation we examine does not capture the potential diet shifts observed in sibling species (Goodheart et al., 2017), novel traits (Liu et al., 2018), and/or micro habitat preference (Whitney, Donahue & Karl, 2018) that could act to drive or reinforce lineage diversification. The latter is particularly intriguing given recent studies demonstrating ecological speciation resulting from shifts in coral hosts (Simmonds et al., 2018). Further studies such as stable isotope analysis and detailed morphological examinations in a phylogenetic framework are needed to better understand the ecological and morphological divergence between these two lineages.