Geographic isolation facilitates the evolution of reproductive isolation and morphological divergence

Abstract Geographic isolation is known to contribute to divergent evolution, resulting in unique phenotypes. Oftentimes morphologically distinct populations are found to be interfertile while reproductive isolation is found to exist within nominal morphological species revealing the existence of cryptic species. These disparities can be difficult to predict or explain especially when they do not reflect an inferred history of common ancestry which suggests that environmental factors affect the nature of ecological divergence. A series of laboratory experiments and observational studies were used to address what role biogeographic factors may play in the ecological divergence of Hyalella amphipods. It was found that geographic isolation plays a key role in the evolution of reproductive isolation and divergent morphology and that divergence cannot be explained by molecular genetic variation.


| INTRODUCTION
Geographic isolation in novel environments often results in rapid ) parallel and convergent evolution Muschick, Indermaur, & Salzburger, 2012). Reproductive isolation has been shown to evolve rapidly in populations adapting to novel environments (Hendry, Wenburg, Bentzen, Volk, & Quinn, 2000), presumably resulting in ecological speciation. However, identifying and quantifying the potentially multifarious processes that contribute to the evolution of reproductive isolation remains a challenge (Garant, Forde, & Hendry, 2007;Nosil, Harmon, & Seehausen, 2009;Nosil et al., 2012;Rundell & Price, 2009). These processes might include ecological, physiological, or morphological adaptation to novel environments, along with biogeographic processes that promote differentiation or limit gene flow. Identifying the contributors to reproductive isolation can be especially difficult in recently diverged, or rapidly diverging lineages, or in lineages that contain cryptic diversity. For example, morphological similarity is not a reliable predictor of interfertility among cryptic lineages. Reproductive isolation has been found to exist within nominal morphological species revealing the existence of cryptic species complexes (Dincă et al., 2013;Gebiola, Kelly, Hammerstein, Giorgini, & Hunter, 2016;Ishikawa et al., 2013;Paterson et al., 2016). Cryptic species complexes may also be paraphyletic in many instances due to selection driving morphological conformity across several unrelated populations occurring in similar habitats Westram, Panova, Galindo, & Butlin, 2016). Despite these challenges, cases of cryptic divergence provide opportunities for study of the evolution of reproductive isolation (Rosenblum & Harmon, 2011) and the development of approaches that can be used to test hypotheses about the factors contributing to reproductive isolation.
Herein, we combine data on morphological and molecular genetic variation with experimental quantification of the strength of reproductive isolation among populations with varying degrees of geographic isolation.
We focused on the freshwater genus of talitrid amphipods, Hyalella (Smith 1874). The nominal species H. azteca (Gonzalez & Watling, 2002) has been found to contain extensive cryptic diversity (Dionne, Dufresne, & Nozais, 2017;Dionne, Vergilino, Dufresne, Charles, & Nozais, 2011;Vergilino, Dionne, Nozais, Dufresne, & Belzile, 2012;Witt & Hebert, 2000;Witt, Threloff, & Hebert, 2006). Also belonging to this genus are numerous morphologically distinct nominal species, each endemic to just a single locality (Baldinger, 2004;Baldinger, Shepard, & Threloff, 2000;Cole & Watkins, 1977;Stevenson & Peden, 1973;Witt et al., 2006). Some of these populations have been found to occur sympatrically with populations of H. cf. azteca (Cole & Watkins, 1977;Stevenson & Peden, 1973;Witt et al., 2006), suggesting that reproductive isolation has allowed the two forms to coexist sympatrically without the endemic form going extinct due to introgression. This assertion is supported by al lack of evidence for niche partitioning among sympatrically occurring populations of Hyalella (Dionne et al., 2017). The presence of cryptic lineages, variation in the degree of geographic isolation among lineages, and the evidence of local adaptation in the narrowly distributed lineages makes Hyalella an ideal system for quantifying the factors that contribute to the evolution of reproductive isolation.
The objectives of this study were to address three questions regarding the patterns of reproductive isolation among Hyalella lineages: (i) "What role does geographic isolation play in determining the patterns of molecular and morphological differentiation?", (ii) "How are levels of reproductive isolation related to morphological and molecular genetic differentiation?", and (iii) "Can variation in reproductive isolation be explained by biogeography?" In order to answer these questions, we surveyed morphological and molecular genetic variation, and experimentally quantified the strength of reproductive isolation among lineages, and the answers to these questions contribute to a foundation for understanding the multifarious processes that shape the evolution of reproductive isolation.

| METHODS
This study used five populations (Table 1) belonging to Hyalella-a widespread and abundant taxon of freshwater amphipods distributed across North America. Two of the populations studied herein fit the expectations of the H. cf. azteca morphotype (as defined by Gonzalez & Watling, 2002) while the other three are noticeably morphologically distinct. However, only one of these populations has been formally described (Hyalella texana Stevenson & Peden, 1973). One of our sampling locations (San Marcos Springs; referred to herein as SMS) was found to have a population of H. cf. azteca co-occurring with a morphologically distinct undescribed spring-endemic species (referred to herein as SMS Hyalella sp.). SMS Hyalella sp. and H. texana are both documented to be endemic to physicochemically stable springs separated by hundreds of kilometers of ambient surface water, suggesting that physiological limitations are responsible for the geographic isolation of these populations. We compared morphological and genetic variation, in combination with attempted mating experiments and study of biogeographic distributions, in an attempt to explain factors contributing to reproductive isolation.

| Establishment of stock cultures
Stock cultures of Hyalella were established to provide a continuous source of live animals for experimentation and to control for the possibility of morphological differentiation due to phenotypic plasticity in situ. Cultures of Hyalella spp. were collected from the four localities listed in Table 1. All amphipods were collected from source localities using dip nets, turkey basters, or a Ponar grab sampler. Cultures were established with at least three separate sampling events for each population between January and August of 2014 and maintained under essentially identical conditions in separate 20-L buckets for each population. Each bucket was given a sand substrate and filled with artesian water with water changes twice monthly. Buckets were maintained at a constant 22°C and kept on a 12 hr/12 hr-light/dark cycle.
All cultures were fed the same diet of Amblystegium sp. and organic detritus ad libitum on a daily basis. Great care was taken when handling cultures to ensure that organisms did not get moved between cultures.

| Quantifying morphological variation
To quantify morphological variation while controlling for potential effects of phenotypic plasticity, cultures from all five populations were T A B L E 1 Collection localities and count of dorsal mucronation for each Hyalella population reared in common-garden replicates. As conditions in all stock cultures (described above) were maintained under the same conditions, we anticipated that five generations would be sufficient to control for environmental or maternal effects; a lack of variation in neonate size across all the populations and generations in captivity suggests that maternal effect was not a factor (Glazier, 2000; Table S3). Therefore, after at least five generations of raising amphipods in stock cultures, morphology was compared between cultures. Twenty individuals were gently wet-mounted (taking care to avoid harming experimental individuals) and photographed at 10× magnification with a calibrated scale bar superimposed on each photograph using an Olympus cellSens camera system and software. Morphometric characters (total length, longest mucronation length, and head length; see Figure S1 for explanation) were estimated from these photographs using Digimizer software (www.digimizer.com). We also counted the number of dorsal mucronations and calculated the ratio of the length of the longest mucronation to total length.
Principal components analysis (PCA) was used to examine morphological characteristics of wild-caught populations. The PCA analysis included three independent variables: head length/total length ratio and longest spine as continuous variables, and spine count as a meristic variable. PCA was conducted in R using the "princomp" function.
Using the same morphological variables, the degree that commongarden populations morphologically differentiated from wild-caught populations, if at all, was assessed. Pairwise permutational multivariate analysis of variance (perMANOVA) was used to test for differences across wild-caught populations, across cultured populations, and between wild-caught and cultured populations using Bray-Curtis distance and the "Adonis" function using the statistical package "vegan" in R. For each perMANOVA analysis, a sequential Bonferroni was applied to the results. The statistical package vegan in R was used for the perMANOVA analyses (Oksanen et al., 2016).
To test for a relationship between morphological and molecular variation among Hyalella, the Euclidean distance between the centroids of each population in PCA space was compared to pairwise Bayesian model-corrected genetic distances using the "ade4" package (Dray & Dufour, 2007) in R. The relationship between morphology and phylogeny was visualized using the R package "phytools" (Revell, 2012).

| Molecular methods
A molecular phylogeny based on the mitochondrial cytochrome C oxidase subunit I (COI) locus was constructed in order to analyze the relationship between morphological similarity, geographic factors, and a history of shared common ancestry. During collection of organisms for the establishment of stock cultures, some specimens were preserved in 95% EtOH and stored at room temperature. DNA was extracted from these individuals (n = 3 per population) by placing all or part of individuals in microtubes containing a chelating resin (Chelex 100, Sigma Aldrich), heated to 60°C for 20 min, then 100°C for 20 min. A fragment of the COI gene was PCR-amplified using the primer TrpPar1 (5′-GTTATATAAACTATTAGCCTTCCAA-3′) paired with either COIaV9 (5′-ACTGCCACAACAGAYAARTAMGACCC-3′) or COIaV10 (5′-ACAGCAACAACAGATAARTARGACC-3′). PCR was carried out with TopTaq (Qiagen) kits in 50μl reactions containing 2.0 μl of template DNA. Cycling conditions were 5 min at 94°C, followed by 35 cycles of 45 s at 94°C, 45s at 51°C, and 60 s at 72°C, followed by 5 min at 72°C. PCR products were gel purified and sequenced with TrpPar1 (COI) using an ABI 3730 automated sequencer.
All sequences generated by this study were queried in GenBank using a BLAST search (blast.ncbi.nlm.nih.gov) which returned 332 Hyalella sequences (including the 15 sequences generated herein) as well as 242 sequences belonging to ten other families of amphipods (a subset of one sequence per amphipod family was randomly selected to serve as outgroups

| Quantifying reproductive isolation
To quantify reproductive isolation, we performed a series of nochoice within (control groups) and between-population (experimental groups) mating experiment where we attempted to achieve all possible combinations with respect to both population source and sex (Table 2). Pairings were established using stock cultures by selecting one female from one population source and selecting one male from the same (control groups) or one male from a different population (experimental groups). Only females that were not brooding eggs or young in their marsupia were selected for the experiments.
Body lengths of all individuals were measured prior to pairing by gently wet-mounting and estimating length with a calibrated reticle.
Females were paired with males that were at least equal in length but not greater than twice as long in order to control for size-assortative effects on mating success (Bollache & Cézilly, 2004). This is a conservative approach to estimating reproductive isolation because reproduction is rarely successful between pairs where males are smaller than females. Some combinations could not be achieved because it was difficult to find suitable males (i.e., males that were larger than the respective female) of the various Hyalella types for female H. texana because H. texana is appreciably larger in size than most other Hyalella species.
An experimental replicate consisted of one male and female pair.
Replicates were placed individually in sealed 150-ml containers. Each container was given the same sand substrate and fed a diet consisting of Amblystegium sp. and organic detritus ad libitum. Cultures were maintained at a constant 22°C, and water in containers was refreshed weekly.
Mating trials were run for 8 weeks and were checked once weekly for the production of offspring. After 8 weeks had elapsed, any pairs that had not reproduced were considered to represent unsuccessful crosses. If free swimming neonates were observed, the adults were removed. The length was measured for each of the neonates and the mean length was used to estimate age of the brood following the equation: A = L-1, where A = age in weeks and L = length in mm (see Tables S3 and S4).
Estimated age was used to estimate the date that hybrid offspring had hatched and the date at which they would become 8 weeks old (since 8 weeks is the age at which most Hyalella species are thought to have finished most of their ontogenetic growth; Strong, 1972). At 8 weeks of age, surviving hybrid offspring were either paired with siblings or individuals other than their parents of one or both of the parental populations to test F 1 ′s for interfertility and backcross fertility. These pairings were allowed to run for 8 weeks and were checked once weekly for the production of offspring.

| Correlates of reproductive isolation
To evaluate potential factors that might explain the occurrence of reproductive isolation, the results from the reproductive isolation experiment were arranged into a matrix. This matrix was compared to a matrix of pairwise Bayesian model-corrected genetic distances.
Matrices were compared using the "ade4" package (Dray & Dufour, 2007) to run Mantel tests in R. To assess the possibility that the relative degree of geographic isolation may potentially lead to reproductive isolation, each population was scored as either reproductively isolated (1) or not (0). Two one-way ANOVAs were used to test for a relationship between the reproductive isolation score and (i) the number of populations of each clade (as determined by molecular analysis), and (ii) the length of reach occupied by each population (Table 3). was computed from the resulting trees (Figure 1). Appreciable molecular divergence was detected within Hyalella with evidence of saturation ( Figure 2). To facilitate discussion of the phylogeny, haplotypes are grouped into clades (Figure 1, Table S5).

| Evaluation of morphological variation
Principal component axes I and II cumulatively explain 89% of the morphological variation in the characters measured. PC I explained 61% of the variance while PC II explained another 28%. The morphological gradient along PC I shifted from negative loadings for longest spine and spine count to positive loadings for head/total length; PC axis II had a gradient of longest spine to head/total length ( Figure 3, Table 4). PerMANOVA (as well as linear discriminant analysis and discriminant function analysis) on wild-caught amphipods showed that all five wild-caught amphipod populations were morphologically distinct (Table 5; Tables S1 and S2). Wildcaught versus the common-garden amphipods showed that all F I G U R E 1 Bayesian phylogeny based on 495-bp region of the COI gene. Terminal nodes represent unique haplotypes. Haplotypes were grouped into clades where applicable. Bayesian posterior probabilities are given at all major nodes common-garden amphipods differed significantly from their wildcaught ancestors except for the SMS population (Table 5). Despite an apparent shift in morphology after five generations in captivity under essentially identical conditions, all common-garden populations remained morphologically distinct from each other at p < .05 (Table 5). Distribution of centroids in PCA space was not found to be significantly correlated to genetic distance (r = 0.75, p = .19; Figure 4).

| Reproductive isolation
After 8 weeks, all conspecific controls had successfully produced offspring while only three of the potential crosses successfully produced offspring (Table 6). In some cases, certain between-population pairings resulted in predation by one individual on their potential mate; this was never observed with any of the within-population pairings. Among the replicates that successfully produced offspring, there was noticeable resistance by the heterospecific pairs to mate. Conspecific control pairs produced offspring as early as 2 weeks into mating trials while none of the successful heterospecific pairs produced offspring until after at least 4 weeks ( Figure 5). This result is consistent with interfertile heterospecific populations having some degree of prezygotic reproductive isolation. After rearing hybrid offspring to adulthood, all hybrid offspring successfully produced offspring suggesting that hybrids are fertile.

| Evaluation of factors contributing to evolution of reproductive isolation
Reproductive isolation was not found to be significantly explainable by genetic distance (r = 0.54, p = .16). However, geography was found to be an important factor (Figure 6) as the number of populations of each clade and the length of reach occupied by each population was both found to significantly explain the occurrence of reproductive isolation ( Figure 6, Table 7).

| DISCUSSION
Appreciable morphological and molecular differentiation were observed in the five populations of Hyalella in this study (Figure 1,   Figure 3). The molecular analysis suggests that ( paraphyletic complex of populations with appreciable molecular divergence between several distinct lineages despite fitting the usual expectations of a morphological species (Figure 1, Witt & Hebert, 2000;Witt et al., 2006;Wellborn & Broughton, 2008;Dionne et al., 2011). Based on the depth of molecular divergence between populations, and the paraphyletic distribution of populations conforming to the H. cf. azteca morphotype, it is likely that this morphology is due to selection more so than common ancestry; although it is unclear if selection is stabilizing, causing the retention of ancestral morphology, or directional causing convergence. The observation of morphological diversity not conforming to an inferred history of shared common ancestry is not unique to Hyalella (Faria et al., 2014;McGee, Neches, & Seehausen, 2016).
The common-garden populations experienced some degree of morphological differentiation from wild-caught ancestors after just five generations in captivity (Table 5). It is likely that differentiation of lab stock occurred via drift or plasticity due to inevitable bottlenecks when establishing populations in captivity, as morphological divergence and local adaptation have been shown to occur rapidly in captivity (Fragata et al., 2014). However, differentiation was convergent toward the H. cf. azteca morphotype which could explain the pervasiveness of this form.
Only some of the populations were found to be interfertile and this did not strongly correlate with history of common ancestry or morphological similarity (Figure 4).  2001;Ellwood & Gose, 2006;Hall & Penner, 2013;Nordt, Boutton, Hallmark, & Waters, 1994;Russ, Loyd, & Boutton, 2000). Therefore, divergence between populations likely occurred in the absence of gene flow; thus, sympatrically occurring populations likely represent secondary contact. It is unclear if divergence occurred directionally due to selection or drift during periods of geographic isolation.
However, molecular distance did not predict morphology or reproductive isolation, but geographic range size was found to be negatively correlated with interfertility ( Figure 6, Table 7). This result is consistent with the hypothesis that geographic isolation drives accelerated divergence (Woolfit & Bromham, 2005) as spring-endemic populations were both found to be completely reproductively isolated and morphologically distinct. Shared morphology among spring All numbers are p-values. The diagonal (bolded) represents the comparison of laboratory-reared common-garden populations with respective wild-caught source populations. Below the diagonal are comparisons across wild-caught populations, and above the diagonal are comparisons across commongarden populations laboratory-reared for five generations.
T A B L E 5 Results of perMANOVA tests of morphometrics from wild-caught and common-garden stock F I G U R E 4 Phylomorpho plot of population centroids with phylogenetic relationship. Genetic similarity is not related to distribution of centroids in principal components analysis space. Decimals along branches represent the Bayesian model inferred number of substitutions  (Cole & Watkins, 1977;Stevenson & Peden, 1973), and repeated parallel adaptations have been shown to occur when selection is strong and similar across different populations Westram et al., 2016). It is possible that the ecological gradient that is responsible for morphological divergence in spring Hyalella is also associated with the evolution of reproductive isolation (Rundle & Nosil, 2005). If so, Hyalella, particularly spring-endemic Hyalella, may represent a good model for the study of ecological speciation.
It is remarkable that Hyalella was recovered as a monophyletic taxon as the depth of divergence between different Hyalella lineages is comparable to the depth of divergence observed among the other amphipod families included in our analysis (depth of divergence between outgroups in Figure 1 is comparable to the divergence found within Hyalella). However, this could be due to a generation time effect leading to different rates of molecular divergence (Chao & Carr, 1993;Ohta, 1993;Thomas, Welch, Lanfear, & Bromham, 2010) as Hyalella has roughly four generations per year-a much faster rate of reproduction than observed in the other amphipod families discussed herein (Crawford & Tarter, 1979;Welton, 1979).
We also only used a single mitochondrial locus because of the abundance of archived COI sequences for amphipods; however, the rate of divergence may be too rapid at the COI locus to properly estimate relationships with such deep divergence (Figure 2). It is possible that divergence between lineages is approaching saturation which ap- It is important to point out that the present study recovered fewer haplotypes than previous authors despite sequencing the same locus and using the same sequences published by other authors on GenBank. This is likely due to the trimming of sequences to much fewer base pairs in order to have a complete alignment as different authors amplified different regions of the COI locus. Therefore, it is likely that variable sites were eliminated that other authors used to identify haplotypes.
Divergence in isolation may lead to scenarios that allow for trait deterioration due to genetic drift (Bromham, 2009;Woolfit & Bromham, 2005), including attributes that affect reproductive isolation. Presumably, there is strong stabilizing selection within a population to maintain interfertility with other members of the same population. In smaller populations, individuals that are divergent in reproductive compatibility have a greater proportional effect on the gene pool of the population (Gillespie, 2001;Woolfit & Bromham, 2005). Therefore, it is more likely that genomic changes that lead to barriers to interfertility will be retained in smaller pop-  (4) Diagonal represents same-population controls.
T A B L E 6 Cumulative percentage of crosses that successfully produced offspring after 8 weeks, with number of replications in parentheses (which was sometimes limited by the availability of individuals size-matched for compatible pairings) F I G U R E 5 Cumulative proportion of successfully reproducing pairs across time. By the second week, conspecific pairs had produced offspring; however, none of the heterospecific crosses produced offspring until at least 4 weeks had elapsed. Only the heterospecific crosses that successfully produced offspring are depicted. None of the heterospecific pairings including H. texana or SMS Hyalella sp. successfully produced offspring light on the factors that contribute to the evolution of reproductive isolation.

ACKNOWLEDGMENTS
We are deeply grateful to Gary Wellborn for his invaluable contributions.
Without the numerous volunteers that assisted with field collections, es-

DATA ACCESSIBILITY
Genbank accessions: see Table S5. Geography in both size of distribution and number of known localities for each haplotype was found to significantly explain the occurrence of reproductive isolation.