Distinct patterns of hybridization across a suture zone in a coral reef fish (Dascyllus trimaculatus)

Abstract Hybrid zones are natural laboratories for investigating the dynamics of gene flow, reproductive isolation, and speciation. A predominant marine hybrid (or suture) zone encompasses Christmas Island (CHR) and Cocos (Keeling) Islands (CKE), where 15 different instances of interbreeding between closely related species from Indian and Pacific Oceans have been documented. Here, we report a case of hybridization between genetically differentiated Pacific and Indian Ocean lineages of the three‐spot dascyllus, Dascyllus trimaculatus (Rüppell, 1829). Field observations indicate there are subtle color differences between Pacific and Indian Ocean lineages. Most importantly, population densities of color morphs and genetic analyses (mitochondrial DNA and SNPs obtained via RADSeq) suggest that the pattern of hybridization within the suture zone is not homogeneous. At CHR, both color morphs were present, mitochondrial haplotypes of both lineages were observed, and SNP analyses revealed both pure and hybrid genotypes. Meanwhile, in CKE, the Indian Ocean color morphs were prevalent, only Indian Ocean mitochondrial haplotypes were observed, and SNP analysis showed hybrid individuals with a large proportion (~80%) of their genotypes assigning to the Indian Ocean lineage. We conclude that CHR populations are currently receiving an influx of individuals from both ocean basins, with a greater influence from the Pacific Ocean. In contrast, geographically isolated CKE populations appear to be self‐recruiting and with more influx of individuals from the Indian Ocean. Our research highlights how patterns of hybridization can be different at scales of hundreds of kilometers, due to geographic isolation and the history of interbreeding between lineages.


| INTRODUC TI ON
Hybridization, or the interbreeding between species or subpopulations, is common in areas where allopatric lineages overlap (Hewitt, 1988). Hybrid zones located at the edge of biogeographic provinces (i.e., suture zones) are areas where regional biotas interact, promoting reproduction between closely related groups (Remington, 1968).
These areas are widely recognized as natural laboratories for investigating the dynamics of gene flow, reproductive isolation, and speciation (Hewitt, 1988;Remington, 1968). Traditionally, hybridization has been associated with hampering differentiation, by allowing the exchange of genetic material (Abbott et al., 2013), which can ultimately lead to the fusion of divergent lineages (Rudman & Schluter, 2016). Hybridization has also been associated with the production of maladapted individuals, leading to evolutionary dead-ends (Barton, 2001). However, more recent studies suggest this process can also have beneficial outcomes, as advantageous traits can be shared between multiple groups (Pardo-Diaz et al., 2012;Runemark, Vallejo-Marin, & Meier, 2019), hybrids have higher fitness than their parental lineages (i.e., hybrid vigor; Chen, 2013), and hybridization can also lead to the formation of new species and generate adaptive radiations (Meier, Marques, Wagner, Excoffier, & Seehausen, 2018;Meier et al., 2019;Seehausen, 2004).
Until recently, hybridization was considered relatively rare in marine environments (Arnold, 1997). However, ongoing ecological and molecular research in marine suture zones has led to a stark increase in the number of confirmed cases of hybridization (DiBattista et al., 2015;He et al., 2019;Hobbs, van Herwerden, Pratchett, & Allen, 2014;Montanari, Hobbs, Pratchett, & van Herwerden, 2016). To date, two major marine suture zones have been recognized for coral reef fishes: the well-studied zone that encompasses Christmas Island (CHR) and Cocos (Keeling) Islands (CKE) in the eastern Indian Ocean Hobbs, Frisch, Allen, & Van Herwerden, 2008), and the less explored zone in the Socotra Archipelago (DiBattista et al., 2015). In the former, glaciations led to sea level drops that exposed the Sunda Shelf, restricting water (and consequently gene flow) between the tropical Indian and Pacific Oceans, which resulted in divergence of many marine populations and species (Barber, Palumbi, Erdmann, & Moosa, 2000;Gaither & Rocha, 2013;Ludt & Rocha, 2015). Since the end of the last ice age, sea levels rose and the previously exposed land masses are now shallow reef habitat, allowing the admixture of Indian and Pacific Ocean lineages.
To date, 15 cases of hybridization have been observed in the Cocos-Christmas suture zone, involving 27 species across eight families ). These cases commonly follow Hubbs' principle, which takes place when one or both parental species exhibit low frequencies (DiBattista et al., 2015;Hubbs, 1955). Previous studies suggest a potential bias when reporting these instances, as most cases belong to sister species that produce hybrids with a clear intermediate coloration, such as butterflyfishes (Chaetodontidae; . Further, 14 out of those 15 hybrid cases are from broadcast spawning species, the other case is for substrate spawners Chromis fieldi and C. margaritifer (He et al., 2019). Damselfishes (family Pomacentridae) are one of the most diverse families of coral reef fishes, and hybridization events are not uncommon (i.e., anemonefishes, Litsios & Salamin, 2014). However, relatively few hybrids have been reported in the Cocos-Christmas suture zone, raising the issue of potential cryptic hybridization. Cryptic hybrids are common in the marine environment (Richards & Hobbs, 2015), and they could represent a significant portion of the cases of hybridization in the Cocos-Christmas suture zone ).
The three-spot dascyllus (Dascyllus trimaculatus, (Rüppell 1829) species complex offers interesting opportunities for understanding the dynamics of speciation and introgression. They can be split  (Bernardi, Holbrook, & Schmitt, 2001;Bernardi, Holbrook, Schmitt, & Crane, 2003;Bernardi, Holbrook, Schmitt, Crane, & DeMartini, 2002;Leray et al., 2010;McCafferty et al., 2002). A widespread population genetic study of this group suggested a combination of historical and ecological factors promoted the divergence of this complex (Leray et al., 2010). The present study aims to elucidate the genomic divergence of the Pacific and Indian Ocean lineages of D. trimaculatus (isolated by limited water circulation in the Sunda Shelf during the Pleistocene; Leray et al., 2010), and their interaction in the Cocos-Christmas suture zone using a combination of genomic and field approaches. We hypothesize that hybridization is occurring between these lineages in the suture zone and that the genomic structure of hybrids will be different for the islands due to their geographic positioning. As CHR is closer to the Pacific Ocean and CKE is closer to the Indian Ocean, our expectation is that hybrid individuals from CHR will have a higher proportion of Pacific ancestry, and conversely, CKE hybrid individuals will have a higher proportion of Indian ancestry. Our study provides an example on the importance of suture zones for genomic exchange between closely related lineages, while highlighting how dynamics of hybridization can be different at small geographic scales (hundreds of kilometers).

| MATERIAL S AND ME THODS
The species D. trimaculatus is a planktivore commonly found in lagoons and reef walls of the Indo-Pacific (Bernardi et al., 2001;Randall & Allen, 1977). They are sexually monomorphic and display complex reproductive behavior where the males attract females with acoustic signals and movements (Fishelson, 1998). Eggs are demersal and the pelagic larval duration is estimated to be 22-26 days (Wellington & Victor, 1989). Juveniles are associated with anemones, urchins, or branching corals. As adults, they shift away from these associations, but typically remain nearby (Randall & Allen, 1977).

| Sampling
Previous studies have indicated significant genetic differences between the Indian and Pacific lineages of D. trimaculatus (Leray et al., 2010). Hence, the sampling design aimed to gather individuals from the suture zone and locations representing pure Pacific or pure Indian Ocean lineages (Figure 1). The Indian Ocean lineage was represented by samples of the Arabian Peninsula (APE), the western Indian Ocean (Europa; EUR), and the central Indian Ocean (Diego de Garcia; DGA). The Pacific Ocean lineage was represented by samples from Indonesia (IND), the Philippines (PHI), and Okinawa, Japan (OKI). Fish were collected between 2010 and 2014, and the samples from IND and three from CKE were part of an earlier study (Leray et al., 2010). All experiments were performed in accordance with UCSC Institutional Animal Care and Use Committee (IACUC/ BERNG-1601).

| Color morphs and visual surveys
We explored whether there were consistent differences in color between Pacific and Indian Ocean lineages across the species' distribution range using personal field observations and photographs. The abundance and color morph frequencies of D. trimaculatus at CHR and CKE were estimated using standard underwater visual census methods. Three replicate 50 × 5 m belt transects were conducted at two depths (5 and 20 m) across eight sites at CHR and at seven sites at CKE.

| Mitochondrial DNA sequencing and data analysis
Genomic DNA was extracted using the Qiagen DNeasy animal blood and tissue kit (Qiagen). A 495 base pair fragment of the mitochondrial control region (D-loop) was amplified in 121 individuals using CR-A and CR-E primers (Lee, Conroy, Howell, & Kocher, 1995). PCRs (10 μl) were performed using 5 μl of multiplex PCR mix (Qiagen), 0.5 μl of CRA 10 μM, 0.5 μl of CRE 10 μM, 3 μl of water, 1 μl of DNA (30-100 ng/μl). Touchdown PCRs were set up as follows: 15 min at 95°C, followed by 20 cycles of 30 s at 95°C, 60 s at 58°C, 90 s at 72°C. During these cycles, the temperature was decreased −0.4°C every minute. This was followed by 15 cycles of 30 s at 95°C, 60 s at 50°C, 90 s at 72°C; and a final extension of 72°C for 10 min. DNA was purified using exonuclease I and FastAP™ thermosensitive alkaline phosphatase (ExoFAP; USB), running it for 60 min at 37°C, followed by 15 min at 85°C. DNA was sequenced in the forward and reverse direction using fluorescent-labeled dye (BigDye 3.1, Applied Biosystems Inc.) using an ABI 3730xl Genetic Analyzer. Sequences were aligned and trimmed using Geneious 6.06 (Kearse et al., 2012), and final edits were performed manually.
In total, 121 individuals from eight populations ( Figure 1) were used in the mitochondrial DNA data analysis. The haplotype and nucleotide diversity of each lineage and of the populations in the suture zone was estimated with Arlequin 3.5.1.2 (Excoffier & Lischer, 2010).
In order to assess genetic divergence between pure lineages and hybrid sites, estimates of pairwise Phi-st (Φ ST ) were calculated with the same software. For the above analysis, we pooled the population samples of each pure lineage in the Indian Ocean (APE, EUR, DGA) and Pacific (IND, PHI, OKI; see Table 1). Significance of Φ ST values was tested with 10,000 permutations. The Kimura 2 Parameter was selected as the best model of sequence evolution in JModeltest2 (Darriba, Taboada, Doallo, & Posada, 2012). Sequential Bonferroni corrections were applied for multiple comparisons (Rice, 1989).
In order to assign each individual haplotype from the hybrid zone to the Indian and Pacific lineages, a minimum spanning haplotype network was constructed using the software PoPArt (Leigh & Bryant, 2015).

| RADSeq library preparation, sequencing, and assembly
Genomic analyses were done on single nucleotide polymorphisms (SNPs), obtained via double-digest restriction site-associated DNA sequencing (ddRAD). The library preparation was performed for 224 samples for two different projects including this study. Library preparation followed the protocol described by Peterson, Weber, Kay, Fisher, and Hoekstra (2012), starting with 500 ng of total DNA per sample. Samples were digested for 3 hr at 37°C, using Hohenlohe, Bassham, Amores, & Cresko, 2013). Raw data were demultiplexed and filtered using the "process_rad_tags" script. Average quality scores were determined within a sliding window 15% the length of the sequence, and reads with any score below 90% of being correct were discarded. Sequences were trimmed to 95 bp, and loci were assembled using the stAcks "de_novo_map.pl" pipeline, using a minimum of three identical reads to create a stack (m = 3), three mismatches allowed between loci within an individual (M = 3), five mismatches when aligning reads (N = 5), and two mismatches when building the catalog (n = 2). We used the "populations" script to obtain loci shared between the eight populations (p = 8), in at least 65% of individuals within a population (r = 0.65) and with coverage of 8× (m = 8). We used only the first SNP of each sequence and removed loci with minor allele frequencies lower than 5%. The filtering resulted in a total of 128 individuals with 2,818 loci, with 83% completeness (see Figure 1 for samples per site). For all downstream analyses, we used the structure output file produced by stAcks which was converted to other file formats using PgdsPider 2.0 (Lischer & Excoffier, 2012).

| Analysis of RADSeq data
A total of individuals from the eight populations mentioned above

Indian Ocean
No significant genetic differences were observed between these localities, so they were combined into regional pools.
To quantify the extent of genetic differentiation between lineages and the hybrid zone, pairwise F ST values were estimated between CHR, CKE, and the Pacific/Indian Ocean lineages using Arlequin with 10,000 permutations. For this analysis, we pooled the population samples of each pure lineage in the Indian Ocean (APE, EUR, DGA) and Pacific (IND, PHI, OKI; see Table 1). In addition, pairwise F ST was calculated per population separately. Sequential Bonferroni corrections were applied for multiple comparisons (Rice, 1989). To further investigate population structure, genetic assignment of individuals of each population was calculated with structure (Pritchard, Stephens, & Donnelly, 2000) using correlated allele frequencies in an admixture model, for 1 million Markov chain Monte Carlo (MCMC) iterations and 10% burn-in. We ran 10 simulations for each group (K = 1 to K = 9). The most likely number of clusters (K) was determined with the Evanno method (Evanno, Regnaut, & Goudet, 2005) using structure HArvester (Earl & vonHoldt, 2012).
Results of the 10 simulations of the most likely K were aligned using cluMPP (Jakobsson & Rosenberg, 2007), and later graphed with distruct (Rosenberg, 2004).
In addition, to visualize the relationship between populations in the suture zone and the Indian and Pacific lineages, we ran a discriminant analysis of principal components (DAPC; Jombart, Devillard, & Balloux, 2010). DAPC describes genetic clusters using a multivariate method that combines discriminant analysis (DA) with principal component analysis (PCA) to best summarize differences between groups (populations) while minimizing variation within groups. DAPC also provides the probability of assignment of individuals to specific groups based on the retained discriminant functions, which can be interpreted as the proximity of individuals to the different clusters (Jombart et al., 2010). The analysis was executed using Adegenet (Jombart, 2008) for R (R Development Core Team, 2015), using the best number of principal components (PCs) identified with the cross-validation method (xValDapc function). To investigate the putative number of clusters (groups) in the data, we applied the "find.
clusters" algorithm, which finds the number of groups maximizing the variation between groups, and compares different clustering solutions using Bayesian information criterion (BIC, Jombart et al., 2010). The DAPC was generated using the number of PCs and clusters identified.
To investigate if the structure signal is driven by SNPs highly differentiated between Pacific and Indian Ocean populations, we ran additional structure analysis using a subset of neutral loci and another subset of highly differentiated loci. To classify loci, we ran the F ST outlier approach implemented in Arlequin. The method consists of a modified fdist procedure under a hierarchical island model that simulates a null distribution across loci as a function of heterozygosity and determines outliers as those that show either significantly  (Excoffier, Hoffer, & Foll, 2009). We ran 50,000 simulations with 100 demes per group, with minimum and maximum expected heterozygosities of 0 and 0.5, respectively. Based on the results of this analysis, we classified each locus into one of three categories: neutral loci (between the 95% quantiles), balancing loci outliers (below the 1% quantile), or divergent loci outliers (above the 99% quantile).
In addition, we used this classification to compare allele frequencies across the sampled populations. We calculated the allelic frequencies of each locus with Arlequin, and results were plotted with ggplot2 in R.
In order to identify fixed loci between Pacific and Indian Ocean populations, we grouped all individuals from the Indian Ocean into a single population (APE, EUR, DGA) and did the same for all samples from the Pacific Ocean (IND, PHI, OKI). Samples from the suture zone of CHR and CKE were excluded from this analysis. We used the "populations" script to identify loci shared between the Pacific and Indian ocean populations (p = 2), in at least 65% of individuals within a population (r = 0.65). We used only the first SNP for each locus, estimated the pairwise F ST between the groups and only kept fixed loci, resulting in 50 loci with F ST = 1. We generated a whitelist of all 50 fixed loci to run a population script on four populations (Pacific, Indian Ocean, CHR, and CKE). We then generated a genepop file, to obtain the allelic frequencies of the 50 fixed loci among the Pacific and Indian ocean lineages, as well as for the two populations in the suture zone (CHR and CKE). Of those 50 loci, 5 were not found in the suture zone populations. We used the remaining 45 loci to estimate the average allelic frequencies between the Indian Ocean, Pacific Ocean, CHR, and CKE.

| Color morph and visual survey results
The photographs obtained in the field, surveys, and georeferenced images/videos from Flickr used in this study revealed consistent pat-  and/or white were observed in the suture zone (Figure 3 Figure A2). The suture zone appears to be more differentiated from the Pacific Ocean than the Indian Ocean lineage ( Figure A3).

| RADSeq analyses
Assignment tests with structure (all 2,818 loci) showed that the most likely number of clusters is K = 2 (∆K = 14,566), and the second most likely is K = 3 (∆K = 489; Figure 6; Figure A4). The individuals from CHR consisted of both pure and hybrid genotypes, and the analysis identified F1 hybrids and admixed individuals.
In CKE, no pure genotypes were found, and all the individuals were backcrosses with >80% of their genotype composition assigning to the Indian Ocean lineage ( Figure 6). The DAPC analysis The eigenvalues of the analysis showed that 70% of the genetic structure was captured by the first two principal components ( Figure 7).
Of the 2,818 loci, we identified 1,949 neutral loci (between the 95% quantiles), 118 divergent loci (high F ST outliers, above the 99% quantile), and 289 balancing loci (low F ST outliers, below the 1% quantile). A total of 462 loci were unclassified (between the quantiles).
structure results with the neutral loci ( Figure A5) and the divergent loci ( Figure A6) showed similar results than the one with all 2,818 loci. Further, the most likely number of clusters remained as K = 2 ( Figures A5 and A6), while analyses with the balancing loci showed no population structure.
The graphs of the allele frequencies showed the expected patterns, as the highly divergent loci revealed a steep cline on the suture zone, neutral loci exhibit a less steep transition, and the balancing ones, showed little variation (Figure 8 and Figures A7-A9).
Of the 50 loci that were fixed between Pacific Ocean and Indian Ocean individuals (F ST = 1), 17 correspond to protein coding regions (Table A3). Allelic frequencies at these loci varied among populations. By definition, allelic frequency was 100% for alternative alleles in the Indian Ocean and the Pacific. For the suture zone, allelic frequencies were very different. Indeed, the CHR population had a higher allelic frequency for the Pacific alleles (frequency of allele q = 0.685, with a narrow range, q = 0.750-0.633, Figure 5). In contrast, for the CKE population, allelic frequencies were highest for the Indian Ocean alleles (frequency of allele p = 0.925, with also a relatively narrow range of frequencies, p = 0.960-0.813, Figure 5). On average, allelic frequency of the most common allele was higher in CKE populations than in CHR populations (p = 0.925 vs. q = 0.685, Figure 5).

| Contrasting patterns of hybridization across taxa
In general, decades of studies in fishes suggest hybridization is facilitated by two components: differential abundance and spawning mode.
Indeed, for hybridization to be facilitated, the abundance of hybridizing species is expected to be very different, with one species being abundant and the other being rare (Hubbs, 1955), and when reproduction occurs, broadcast spawning is expected to increase the likelihood of successful hybridization (Nydam & Harrison, 2011). This pattern is observed in the Cocos-Christmas suture zone, as 14 out of the 15 hybrids are broadcast spawners and usually one of the two species involved in the hybridization is rare  It is possible that the three-spot dascyllus population of CKE is largely relying on self-recruitment at present, as it is located ~900 km from CHR and more than 1,000 km from Indonesia or any other land-mass (Hobbs, Newman, et al., 2014). Based on biogeographic observations, it has been suggested that populations of reef fish at CKE are likely maintained by self-recruitment (Hobbs, Newman, et al., 2014). This idea is also supported by the Dascyllus trimaculatus molecular indices, as the gene diversity was the lowest compared to other Indian Ocean locations, and by the DAPC results, as CKE was characterized as a separate cluster. The D. trimaculatus depth limit is 55 meters (Allen, 1991) travel large distances (~1,000 km) and settled individuals need to use the scarce available habitats (Hobbs, Newman, et al., 2014). Larvae from elsewhere may not reach the island frequently or die in transit, due to lack of food or complex current reversals surrounding this area (Nieblas et al., 2014). Alternatively, recently settled larvae arrive but only hybrids survive, due to local environmental or ecological conditions. Overall, CKE may be colonized by infrequent vagrants from pure populations, but not enough to maintain a genomic signature of ongoing hybridization or to show genetically pure lineages in our population sample. Self-recruitment at CKE is likely leading to the formation of a hybrid lineage (as shown in Figure 7) that may continue to evolve into a new species. The evolution of new species from hybridization (hybrid speciation) has been shown for freshwater fishes (Abbott et al., 2013;Meier et al., 2017) and it has been suspected for some anemonefishes (Litsios & Salamin, 2014).
Interestingly, two mitochondrial haplotypes were found in CHR, whereas only one mitochondrial haplotype was found in CKE. The CKE haplotype corresponds to the one in the Indian Ocean lineage.
Further, there appears to be mitonuclear discordance as there are, at CHR, individuals with predominant SNPs from the Indian Ocean that have mitochondrial haplotypes from the Pacific (Figure 4). The opposite is also true (i.e., predominant Pacific SNPs with Indian Ocean mitochondria; Figure 4), which further suggests there is rampant hybridization in the suture zone at CHR, but not at CKE. This has been previously reported in marine fishes, where hybridization can lead to different mitonuclear combinations across closely related groups (Bernal, Gaither, Simison, & Rocha, 2017). These events will continue until there is some selective advantage maintaining certain combinations of mitochondrial haplotypes with specific nuclear backgrounds (Hill, 2016), which could already be happening at CKE.

| Evolutionary implications
The Cocos-Christmas hybrid zone is a product of secondary contact; the seaway connecting Indian and Pacific populations has opened and closed several times during glacial cycles. Many other sister species with Pleistocene origins overlap at this location (Gaither & Rocha, 2013). Our data do not show any evidence of introgression beyond CHR and CKE due to restricted gene flow.
The Sunda Shelf barrier might be strong enough to restrict gene flow between Pacific and Indian Ocean populations, even during high sea level stands, and maintain their differences. Additionally, the Pacific Ocean lineage may be maladapted in the Indian Ocean, and vice versa (Hobbs & Salmond, 2008), making their expansion improbable. Finally, since multiple hybrids and second-generation backcrossed individuals were found in the suture zone, it is possible that geographic isolation of the two lineages is what maintains the difference.
In addition, sexual selection may play a role in the maintenance of the hybrid zone. Some Dascyllus species have UV vision capabilities (Losey, 2003) and there are pomacentrids that have dis-

| Cryptic hybridization
The D. trimaculatus complex is a group of species and divergent lineages that lends itself well to studying cryptic hybridization (Leray et al., 2010 (Hobbs & Salmond, 2008) and in areas of Indonesia (Gaither & Rocha, 2013

| CON CLUS IONS
Our research shows that in two locations of a single suture zone, frequency of hybridization is variable, depending on the geographic iso-

ACK N OWLED G EM ENTS
We want to thank multiple funding sources for making this project fi- We also thank D Wagner and SA Jones for assistance with field collections. We thank T. Sinclair-Taylor and K. Willshaw for pictures at Cocos and Christmas; E. Tandjaja, T. Sinclair-Taylor, J. DiBattista for samples and field observations; P. Saenz-Agudelo for assistance in the laboratory; J. Chadha for graphic design, and I. Fernandez, C.
Prada, and H.T. Pinheiro for insightful discussions. We thank D.
Hogan and anonymous reviewers for the comments that greatly improved the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

Population Frequency
Allele_Cat a b F I G U R E A 7 Allele frequencies of a random subsample of 25 loci with the highest Fst (i.e. outliers potentially under divergent selection). The numbers at the top of each graph are the names of the loci. In the Y-axis, the populations are represented by numbers: Arabian Peninsula (APE, 1); = Europa, Scattered Islands (EUR, 2); Diego Garcia, Chagos Archipelago (DGA, 3); Cocos (Keeling) Islands (CKE, 4); Christmas Island (CHR, 5); Indonesia (IND, 6); Philippines (PHI, 7); Okinawa (OKI,8). Dashed lines indicate the location of the suture zones

Population Frequency
Alelle_cat a b F I G U R E A 9 Allele frequencies of a random subsample of 25 loci with the lowest Fst (i.e., outliers potentially under balancing selection). The numbers at the top of each graph are the names of the loci. In the Y-axis, the populations are represented by numbers: Arabian Peninsula (APE, 1); = Europa, Scattered Islands (EUR, 2); Diego Garcia, Chagos Archipelago (DGA, 3); Cocos (Keeling) Islands (CKE, 4); Christmas Island (CHR, 5); Indonesia (IND, 6); Philippines (PHI, 7); Okinawa (OKI,8