Species integrity, introgression, and genetic variation across a coral reef fish hybrid zone

Abstract Hybridization and introgression are evolutionarily significant phenomena breaking down species boundaries. “Hybrid zones” (regions of species overlap and hybridization) enable quantification of hybridization frequency and examination of mechanisms driving and maintaining gene flow. The hybrid anemonefish Amphiprion leucokranos is found where parent species (A. chrysopterus; A. sandaracinos) distributions overlap. Here, we examine geographic variation in hybridization and introgression, and potential impacts on parent species integrity through assessing relative abundance, social group composition, and genetic structure (mtDNA cytochrome b, 21 microsatellite loci) of taxa at three hybrid zone locations: Kimbe Bay (KB) and Kavieng (KA), Papua New Guinea; the Solomon Islands (SO). Relative abundances of and size disparities between parent species apparently drive hybridization frequency, introgression patterns, and genetic composition of taxa. Conspecific groups are most common in KB (65%) where parent species are similarly abundant. Conversely, mixed species groups dominate SO (82%), where A. chrysopterus is more abundant. Hybrids most commonly cohabit with A. sandaracinos in KB (17%), but with A. chrysopterus in KA (22%) and SO (50%). Genetic differentiation (nDNA) analyses indicate that parent species remain distinct, despite ongoing hybridization and hybrids are genetically similar to A. sandaracinos—resulting from persistent backcrossing with this smallest species. This study shows that hybridization outcomes may depend on the social and ecological context in which taxa hybridize, where relative abundance and disparate size of parent species explain the frequency and patterns of hybridization and introgression in the A. leucokranos hybrid zone, reflecting size‐based dominance behaviors of anemonefish social groups.


| INTRODUC TI ON
Hybridization among closely related species is common and can play a significant role in evolution and speciation (Mallet, 2005). Originally thought to be an evolutionary dead end (Coyne & Orr, 2004;Dowling & Secor, 1997;Mayr, 1942), it is now clear that hybridization can contribute to evolutionary change in various ways, most notably through introgression where genetic information from one species transfers to another through repeated backcrossing (Abbott et al., 2013).
Hybridization can promote evolutionary novelty within a system faster than through mutation alone (Grant & Grant, 1994;Kunte et al., 2011). The outcomes of hybridization events are diverse, include fusion of species, reinforcement of reproductive barriers, and generation of new distinct populations of mixed ancestry, and may provide the foundation for speciation and diversification to occur (Abbott, Hegarty, Hiscock, & Brennan, 2010;Mallet, 2007,;Meier et al., 2017;Servedio & Noor, 2003;Taylor et al., 2006;Via, 2009;Wu, 2001). Studying young hybrid taxa therefore allows contemporary insights into potential speciation events or species coalescence in progress, occurring at secondary contact zones between closely related taxa that may be undergoing rapid adaptive radiations (Gourbiere & Mallet, 2010;Meier et al., 2017;Price & Bouvier, 2002;Seehausen, 2004). However, in nature the mechanisms driving and maintaining hybridization, determining patterns of introgression, and maintaining species integrity remain poorly understood.
Hybrid zones provide natural laboratories for studying hybridization and investigating patterns of variation among hybridizing species. Hybrid zones may vary spatially and temporally, with taxa subjected to demographic processes in which novel ecological opportunities may arise (Abbott et al., 2013). Ecological factors often associated with hybridization include abundance disparities between closely related taxa and the shared use of a limited resource (i.e., host, food source, and habitat). As such, the causes and consequences of hybrid zones are complex and varied, and patterns of gene flow represent single observations in time of a dynamic interaction between species (Abbott et al., 2013).
The line of convergence between Indo-Australian and Pacific plates from north-western Papua New Guinea (PNG) to the Solomon Islands (SO) represents a third, but lesser known suture zone ("PNG-Solomon Islands suture zone") where the ranges of many sister species also overlap and taxa hybridize (Gainsford, van Herwerden, & Jones, 2015;Hobbs, van Herwerden, Pratchett, & Allen, 2013;McMillan, Weigt, & Palumbi, 1999). For example, the two butterflyfish species Chaetodon punctatofasciatus and C. pelewensis commonly hybridize here, where McMillan et al. (1999) found a greater frequency of hybrid phenotypes in comparison with parental phenotypes, suggesting greater fitness of hybrids to parental species within the hybrid zone. The "PNG-Solomon Islands suture zone" has a dynamic history of disturbance associated with climatic changes and sea level fluctuations. It falls within the eastern part of the Coral Triangle-the global center of marine biodiversity (Hughes, Bellwood, & Connolly, 2002)-where many closely related species share habitats, increasing hybridization opportunities. Thus, the 'PNG-Solomon Islands suture zone' can provide unique insights into processes promoting hybridization between cohabiting species and how hybridization affects biodiversity on coral reefs.
Anemonefishes are an evolutionarily young, rapidly diversifying group that is prone to hybridization (Santini & Polacco, 2006;Timm, Figiel, & Kochzius, 2008), providing an ideal system to test evolutionary questions on hybridization (Abbott et al., 2013). Anemonefish groups are structured based on size, where individuals queue to breed (Buston, 2004;Buston & Cant, 2006). Females are largest and dominant, followed in size by subdominant males, and progressively smaller nonbreeding subordinates (Fricke, 1979;Hattori, 1991). In the 'PNG-Solomon Islands suture zone' hybridization occurs between anemonefish species Amphiprion chrysopterus and Amphiprion sandaracinos which cohabit the same anemone host species (Gainsford et al., 2015). Due to distinctive colouration, the hybrid of these species was initially described as a nominal species, A. leucokranos, but later confirmed to be a hybrid based on intermediate morphology along with genetic and ecological similarities (Gainsford et al., 2015). The parent species have predominantly allopatric distributions but cohabit and hybridize within a narrow area of overlap, hereafter termed the A. leucokranos hybrid zone. Size differences between hybridizing species in the context of anemonefish hierarchical behavior and protandrous hermaphroditism were most important in driving ecological and evolutionary patterns observed in this hybridization. The larger A. chrysopterus always mates as the female when reproducing with the smaller A. sandaracinos (Gainsford et al., 2015). Factors such as abundance disparities, overlapping patterns of resource use, and breakdown in assortative mating all promote hybridization in marine fishes (Montanari et al., 2016); however, these factors may vary across the hybrid zone. To understand which underlying mechanisms maintain the hybrid zone and integrity of hybridizing species, knowledge of the geographic abundance patterns of parent species; levels of cohabitation of parent species and hybrids within anemone hosts, and patterns of backcrossing and introgression are required. If relative abundances of hybridizing species differ or cohabitation levels change, the prevalence of hybrids and introgression levels between species is likely to differ strikingly.
Furthermore, geographic variation in genetic differentiation levels between parent species may impact on interbreeding propensity and/or the magnitude of introgression between species.
The overall aim of this study was to investigate geographic variation in abundance, cohabitation, phenotypic characteristics, and genetic composition of parent species and hybrids across the A. leucokranos hybrid zone. A combination of ecological observations, phenotypic measurements, and genetic analyses (using mitochondrial and nuclear DNA markers) was applied to answer five specific questions: (1) How does relative abundance of parent species, hybrid frequency, and cohabitation among taxa vary across the hybrid zone? (2) Are hybrid phenotypic characteristics and likely patterns of backcrossing based on phenotypic characters consistent across the hybrid zone? (3) Does host anemone use by these taxa vary across the hybrid zone? (4) Are patterns of historic (mtDNA) and contemporary (nDNA) genetic structure among parent species and hybrids consistent across the hybrid zone or do patterns differ? and (5) How does genetic structure of taxa across the hybrid zone relate to parent species regional abundances? Answering these questions will inform of mechanisms promoting and maintaining hybridization, pathways of introgression, and likely consequences thereof to the resilience of hybridizing taxa into the future.

| Study taxa and locations
This study was conducted at sites within the A. leucokranos zone Reef in Australia, eastward to French Polynesia (Fautin and Allen, 1997). The A. leucokranos hybrid zone is found where these parent Fish were captured using hand nets, anaesthetized in situ using clove oil, and released into their home anemone postsampling.

| Abundance and cohabitation
All individuals of the parent species A. chrysopterus and A. sandaracinos, and hybrid A. leucokranos observed at the three locations were recorded. Levels of cohabitation were compared by recording: (1) The number of individuals of conspecific A. chrysopterus groups, (2) conspecific A. sandaracinos groups, (3) hybrid only groups, (4) heterospecific A. chrysopterus and A. sandaracinos groups, (5) heterospecific A. chrysopterus and hybrid groups, (6) heterospecific A. sandaracinos and hybrid groups, and (7) groups containing both parents and hybrids.
Additionally, habitat use by parent species and hybrids were characterized among regions within the hybrid zone. Most individuals were encountered between 1 and 20 m depth, where depth, host anemone species, immediate surrounding habitat, and reef zone (reef flat, crest, and slope) were recorded for all groups examined in this depth zone. For each individual captured, the following data were recorded: phenotype (photographed), total length (measured to the nearest mm), and sex (assigned based on relative social position). The presence of egg clutches was recorded when observed, and putative parent species identified.

| Phenotypic characteristics
The relative frequency of hybrid phenotypes was calculated from photographs using seven qualitative traits including: tail shape and color, dominant body color, presence and completeness of dorsal stripe and side bars, as well as latitudinal body shape (see Table S1 for phenotypic categories ranging from pure A. chrysopterus to hybrids, and pure A. sandaracinos, and Table S2 for relative frequency of qualitative phenotypic traits).

| Genetic structure
The population genetic and phylogenetic structure within and outside the hybrid zone was compared to assess regional variation in  individual was fin-clipped for genetic analyses. Small (4 mm 2 ) caudal fin clips were taken from all captured fish and preserved in 80% ethanol. Samples of "pure" parental species from outside the hybrid zone were included to allow for comparison of species-specific genetic signals. Both mitochondrial cytochrome b and nuclear microsatellite markers were employed to estimate historical and contemporary gene flow, respectively.

| Laboratory protocols for genetic analyses
For all laboratory methods described herein, genomic DNA was isolated from fin clips using a standard salting out protocol (Sunnucks & Hales, 1996), and polymerase chain reaction (PCR) products were column purified with GE Illustra Sephadex G-50 for sequencing.
Forty-two Amphiprion spp. microsatellite markers, including 8 novel loci, were tested on seven to eight individuals each of A. sandaracinos, A. chrysopterus and hybrid taxa. Novel primer development and cross-amplification success of markers is detailed elsewhere (Gainsford, Jones, Gardner, & van Herwerden, 2020). Of all markers tested, 23 highly polymorphic loci that consistently cross-amplified in all study taxa across regions were used in optimized multiplex reactions, based on locus sizes using Multiplex Manager 1.0 software (Holleley & Geerts, 2009). PCRs of seven multiplex sets of two to six markers were carried out in 10 µl reactions with 50 ng template, 2X Type-it Multiplex PCR Master Mix (QIAGEN), and 2 µM each primer (forward and reverse). PCR products were visualized by gel electrophoresis using 2.0% agarose, purified as above, and genotyped on an ABI 3730XL Genetic Analyser (Applied Biosystems) with GeneScan LIZ-600bp size standard.

| Data compilation and analyses
Cytochrome b sequences were MUSCLE aligned (Edgar, 2004a(Edgar, , 2004b and manually edited in Geneious v9.0.4. An alignment including sequences from all regions sampled, including those outside the hybrid zone, was used to estimate phylogenetic evolutionary history of taxa and relationships among haplotypes. The best substitution model for the alignment was the HKY + G model chosen from 21 models using a likelihood approach under default settings in MEGA6 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013). Phylogenetic relationships were inferred using standard approaches including maximum parsimony (MP) and maximum likelihood methods in MEGA6 (Tamura et al., 2013), and Bayesian inference (BI) using the MrBayes 3.2 plug-in (Ronquist et al., 2012) through Geneious ( Figures S2-S4). For all analyses, the HKY + G substitution model was implemented, and trees were outgroup rooted using individuals from Amphiprion ocellaris (DQ343956-7, KF264293-4). MP analyses included 10 independent runs using 1,000 bootstrap replicates, with all ten best MP trees recovered having identical length and topology. ML analyses were performed using 1,000 bootstrap replicates under a likelihood approach, and BI analyses were conducted with 1,100,000 iterations and 100,000 tree burn-in.
Structure was run using the admixture ancestry model informed by location, with correlated allele frequencies for each K value for 10 individual repetitions, at 1,000,000 MCMC iterations following a 100,000 burn-in. Structure Harvester (Earl & von Holdt, 2012) was used to assess the best K following Evanno's method (Evanno, Regnaut, & Goudet, 2005). To visually assess relationships between predefined population clusters, a discriminant analysis of principle components (DAPC) was executed using the adegenet package (Jombart, 2008)

| Relative abundance and patterns of cohabitation
Across the hybrid zone, the relative frequency of hybrid individuals was comparable among the three surveyed locations (22%-30%;  at Solomon Islands (12% and 5%, respectively; Figure 4). Therefore, the proportion of heterospecific groups, and thus the incidence of hybridization, varied across the hybrid zone.
The formation of hybrid only groups was low and varied across the hybrid zone: between 8% (Kavieng) and 1% (SO

| Phenotype variation and relative frequency
Relative frequency of phenotypic traits revealed considerable regional variation in both A. chrysopterus and hybrids, but minimal variation across A. sandaracinos populations (  Figure S1).

| Genetic variation across the hybrid zone
Four hundred and thirty mtDNA cytochrome b sites were resolved for 388 individual anemonefishes (Table 1)

| Historical population genetic structure (mtDNA)
The level of population differentiation was high for pairwise population comparisons of cytochrome b (Tables 2 and 3

), revealing
A. sandaracinos to be most differentiated from other taxa examined.
All "pure" A. sandaracinos populations appeared highly differentiated from all other populations, with the exception of Kavieng and Kimbe Bay A. sandaracinos, which were not significantly different (Table 2).  A. chrysopterus, estimated under a median-joining algorithm. Each "pie" represents an individual haplotype, the size of which is proportional to the total number of individuals sharing that haplotype, where individual population identity is indicated by color according to the key (see inset). Substitutions separating haplotypes are indicated in the legend for one, five, ten, and twenty substitutions, respectively. Phylogenetic relationship structure is inferred from MP and ML bootstrap support values (%), and BI posterior probabilities (HPD, maximum value of 1). See Table 1 (Fu, 1997). High haplotype diversity (h > 0.5) and low nucleotide diversity (π < 0.5) were recorded among parental and hybrid populations across the hybrid zone (Table 1), consistent with recent population expansions. Together these results provide evidence for a historical bottleneck followed by population expansion in the hybrid zone.

| Contemporary genetic structure (nDNA)
Summary statistics for 21 microsatellite loci are presented in Table S3. Sample codes used are as above in Table 1 caption.  (Table 3). Although significant, variation among populations within species groups and among individuals within populations contributed only 5.01% and 2.15%, respectively, to overall variation (Table 3).

| Genetic structure relative to parent species abundance
The degree to which parent species nDNA contributed to hybrid populations varied regionally regardless of the relative abundance of species. In Kimbe Bay, where the abundance of each pure parent species and hybrids were near equal (Figure 3), there is an asymmetric 25:75 contribution by A. chrysopterus and A. sandaracinos to hybrid populations (Figure 5c). In Kavieng and Solomon Islands, there is a near 50:50 input by A. chrysopterus and A. sandaracinos (Figure 5c), despite relatively high and low abundance, respectively, of A. sandaracinos compared to A. chrysopterus at these two locations (Figure 3).

| D ISCUSS I ON
Regional disparities in parent species frequency and inherent size disparities between hybridizing species drive variation in the genetic structure among taxa across the A. leucokranos hybrid zone.
The relative abundance of parent species and hybrids varied across the hybrid zone regionally and observed levels of cohabitation did not reflect a scenario whereby rare species 'seek out' heterospecific mates in the absence of conspecifics. Subsequently, hybrid phenotypes were highly variable across the hybrid zone, reflecting the de- Collectively, results suggest the hybrid (originally described as A. leucokranos) is no less fit than the parent species are and may persist in the hybrid zone to differentiate completely from parent species over time. This study shows that the outcome of hybridization is dependent on the social and ecological context in which taxa hybridize.

| Regionally disparate abundance and cohabitation
In Kavieng and Solomon Island regions, where abundance disparities between parent species are evident, significantly more mixed species group assemblages occur than in Kimbe Bay, where conspecific groups are twice as common. In contrast, the frequency of each parent species in Kimbe Bay is relatively equal and overall conspecific assemblages predominate. Abundance disparities between species are considered a key factor facilitating hybridization between sister taxa in regions of range overlap (Hubbs, 1955). In a recent review of fish hybridization, rarity of one or both parent species was reported as the primary ecological factor implicated in promoting hybridization among marine fishes (Montanari et al., 2016). This was followed by shared resource use, specifically the degree of habitat and dietary overlap; however, these factors are not often empirically tested and rather proposed to explain this phenomenon. Mate choice experiments on hybridizing marine fishes are not currently available; however, experimentally altering the relative abundance of two largely sympatric grasshopper species increased hybridization propensity when relative frequencies of sister taxa were increasingly disparate, due to additional inter-species encounters (Rohde, Hau, Weyer, & Hochkirch, 2015). Authors concluded that abundance disparities are a major driver of hybridization and experimentally found for the first time that hybridization probability increased with decreasing relative frequency of conspecific taxa (Rohde et al., 2015). Hybrid systems in which one species is rare and the other abundant are widely reported, where rare species are generally purported to choose mates from an abundant sister species in the absence of conspecifics (Allen, 1979;Hobbs & Allen, 2014;Marie et al., 2007;Montanari, Hobbs, Pratchett, Bay, & van Herwerden, 2014;Moyer, 1981;Randall, Allen, & Steene, 1977;van Herwerden et al., 2006). The data show that common species mate with less common species. Pyle and Randall (1994) asserted that the general assumption of rare species seeking out heterospecific mates does not consider why individuals from a common species might choose to mate with individuals from a rare species when conspecifics are abundant.
It was suggested that particular social systems may provide alternative opportunities for reproduction at more favorable times for dominant individuals of a particular sex (Pyle & Randall, 1994), such as in the harem forming Centropyge species that hybridize (Kosaki, Pyle, Randall, & Irons, 1991;Lutnesky, 1992aLutnesky, , 1992bMoyer, 1981Moyer, , 1990.
However, in Centropyge spp., gender frequency disparities are apparently more important drivers than species abundance disparities.
Here, we propose that in the A. leucokranos hybrid zone, where abundance disparities clearly appear to be associated with hybridization propensity, the underlying reason that abundant A. chrysopterus choose to hybridize may be associated with demand for a limited resource-the host anemone on which groups live and reproduce. In this hybrid zone, intraspecific competition for limited host anemones is great, and the larger species in a given scenario holds a significant size advantage when joining and living in mixed groups.

| Drivers of population structure across this hybrid zone
What is driving the structure found across the A. leucokranos hybrid zone, where abundance disparities appear to promote hybridization? In considering preferences for conspecific or interspecific group formation, it is widely assumed that all individuals have equal choice in determining breeding partners. However, the assumption that mate choice is a level playing field in hybridization between species in hierarchical groups is fundamentally false, as the factor on which dominance depends may not be equally distributed among taxa (Bronson, Grubb, Sattler, & Braun, 2003;Reudink, Mech, & Curry, 2006). In the case of anemonefish, dominance is dependent on size, and anemonefish are well known for living in hierarchical groups in which size dominance determines an individual's right to reproduce as either a female or male (Buston, 2003(Buston, , 2004Buston & Cant, 2006;Fricke & Fricke, 1977;Moyer & Nakazono, 1978 (Ang & Manica, 2010;Keller & Reeve, 1994;Reeve & Keller, 2001;Vehrencamp, 1983;Wong, 2011). In this way, moderately sized hybrids and small A. sandaracinos are disadvantaged against the more dominant species and must continue to queue in the hope of reproducing, rather than facing eviction and becoming vulnerable to mortality outside the group (Buston, 2004;Wong, Munday, Buston, & Jones, 2008). Accordingly, these individuals may benefit from remaining in queues to maintain their reproductive position, through reduced eviction risk and shorter queue times, rather than recruiting to other groups (Mitchell, 2005;Wong, Buston, Munday, & Jones, 2007), and that may compensate for and potentially overcome size disadvantages in fish social hierarchies (Alcazar, Hilliard, Becker, Bernaba, & Fernald, 2014

| Persistence of hybrid A. leucokranos
High haplotype and low nucleotide diversities throughout the hybrid zone suggest that hybridizing species have historically experienced a population bottleneck followed by rapid population growth, which has led to an accumulation of mutations (Avise, Neigel, & Arnold, 1984;Grant & Bowen, 1998;Rogers & Harpending, 1992 & Palumbi, 1995;Palumbi, 1994;. Anemonefish species studied in the Indo-Pacific are highly related based on morphometrics, phylogenetic, and population genetic data (Gainsford et al., 2015;. Specifically, A. leucokranos hybrids regardless of region have greater genotypic diversity than parent species. Hybrid zones, where the process of divergence is underway, offer insights into the importance of biogeography and ecology in shaping population histories and future evolutionary patterns. Despite drawing the conclusion that A. leucokranos may be a true species, Santini and Polacco (2006)  are also evident, albeit rare, due to the size-dominant behavior structuring anemonefish groups.
There is a strong case for recognizing the status of this hybrid as more than an evolutionary dead end in light of overwhelming evidence for the importance of hybridization to the two parent taxa, as well as the indication that A. leucokranos appears to be differentiating from parent taxa. Recently, authors have highlighted the importance of acknowledging hybrid species in light of legislation which is inherently vague and generally does not consider protection or conservation policy measures (Allendorf, Leary, Spruell, & Wenburg, 2001;Chan, Hoffmann, & van Oppen, 2019;Richards & Hobbs, 2015). Losses of taxonomic evolutionary novelty and phylogenetic diversity, as well as increased species extinctions are predicted from inadequate management of hybridization and hybridizing lineages (Chunco, 2014;Dowling & Secor, 1997;Forest et al., 2007;Van Dyke, 2008). Pertinently, A. leucokranos, a highly prized aquarium trade species that is iconic, rare and easily caught due to reliance on sessile anemone hosts, is likely to be detrimentally impacted by removal of its current species status.
Taxonomic delisting of this species, already prized by aquarium traders, may lead to increased harvest and rarity of this already locally rare and endemic taxon, simultaneously driving an increase in market value of individual fish and hence greater motivation for trade. Richards and Hobbs (2015) concluded that in order to conserve coral reef biodiversity, and the processes that are implicit in initiating and maintaining biodiversity, such as hybridization, policies regarding conservation and management must be addressed on an individual case basis, as removal of species status or lack of protection may indirectly impact evolution and biodiversity of species overall.
Extensive investigation of the Amphiprion leucokranos hybrid zone revealed that parent species abundance and size disparities drive regional ecological patterns and genetic structure among taxa. The size of parent species, rather than the species itself, better explains the historical and existing genetic structure, reflecting the characteristic size-based dominance behavior of anemonefish.
This study demonstrated that rare species may not always choose to hybridize with abundant species when abundance disparities arise, such as along the edges of their biogeographical distributions. High haplotypic diversity and low nucleotide diversity in all populations examined suggest a bottleneck followed by recent population expansion that has led to initiation and persistence of this hybrid zone, where the hybrid A. leucokranos appears to be differentiating from the parent taxa. This study emphasizes the need and importance of protection for hybrid species. Not only are A. leucokranos vulnerable to over-harvesting by aquarium traders, but they are also important contributors to both the evolutionary resilience of hybridizing parent species and the biodiversity of coral reef systems.

ACK N OWLED G EM ENTS
The authors are grateful for the logistical support of Mahonia Centre at James Cook University. James Cook University animal ethics approval numbers A1633 and A1791.

CO N FLI C T O F I NTE R E S T
The authors of this manuscript declare no conflict of interest, financial, or otherwise.