Genetic introgression, allele exchange across species boundaries, is a commonly recognized feature of animal evolution. Under such a paradigm contemporary contact zones provide first-hand and complementary insight into the geographic, phenotypic, and genetic details of introgression. Also, when mate choice phenotypes are conspicuous and variable in hybrids, contact zones provide potential insight into how sexual selection interacts with species boundary maintenance, particularly when postzygotic reproductive isolation is weak. The Habronattus americanus subgroup includes several recently evolved jumping spider species, with an estimated age of about 200,000 years, and substantial evidence for hybridization and introgression. We explored a contact zone involving H. americanus (Keyserling, 1885) and H. kubai (Griswold, 1979) on Mount Shasta, California, in alpine habitats that would have been unavailable (under ice) at the Last Glacial Maximum. We characterized morphological diversity within the contact zone, including the fine-scale geographic distribution of hybrid and parental individuals, and assessed genetic variation using ddRADseq data. Combined results indicate a lack of measured genomic differentiation between specimens with distinct morphologies, including individuals with phenotypes of the parental species. We identified a diverse array of hybrid morphologies, with phenotypic evidence for backcrossing, essentially forming a phenotypic bridge between parental taxa. The study area is characterized by more hybrid than parental individuals, with a significantly larger number of red-palped morphologies than white and/or yellow-palped morphologies; the novel, white-palped phenotype is perhaps transgressive. Overall, these results contribute to a better understanding of the expected ebb and flow of lineage interactions during the early stages of speciation.

Specimen Collection. The Mount Shasta study site is located just west of Panther Meadows in the Shasta-Trinity National Forest, at an elevation of approximately 2300 meters, on the south flank of Mount Shasta. We surveyed for spiders in an approximately 270,000 square meter area. Habitat included a southward extending open area of mostly treeless high-elevation pumice field, with herbs and low shrubs, and adjacent high-elevation coniferous forest. A small, paved road winds through the study site. Habronattus americanus subgroup members prefer open habitats, so we spent most survey time out of the coniferous forest. We did not sample the more easterly Panther Meadows proper, with meadow edge microhabitats that we typically associate with H. kubai. Special permits were not required to collect spiders from this location.

Evidence of hybridization at this location was previously identified by the collection of male specimens with combined traits of both parental species. In July 2018 a team of three persons collected 89 specimens from across the study site, including 78 adult males, 9 females, and two penultimate males. We visually scanned ground habitats at random and attempted to collect all adult males seen (hand-collected into glass vials); because females of the two different species cannot be morphologically separated, we limited our collection of adult females. All specimens were georeferenced where collected using a Garmin Etrex 20x handheld GPS device.

Morphological Data Collection and Analysis. Eleven discrete morphological characters were scored for all 78 adult male specimens. Our sample includes individuals described as H. kubai and H. americanus, with specimens possessing character combinations typical for these species (see description above and Bougie et al. 2021). The remaining specimens are described as H. americanus x H. kubai hybrids. We scored a subset of morphological characters from Bougie et al. (2021), focusing on characters that varied within the Mount Shasta samples. Characters were scored by examining individuals in 100% ethanol under a dissecting microscope. We scored hue (not vibrancy) for some characters, distinguishing red from white from golden/yellow. Acknowledging that color can fade for specimens preserved in ethanol, the specimens that we scored were relatively fresh (2018 collections) and preserved in 100% EToH at -80°C; we found that most of the pigmentation hue was readily scorable for these specimens.

To summarize morphological variation and identify morphological clusters we performed a non-metric multidimensional scaling (NMDS) analysis of Bray-Curtis dissimilarity values, using the metaMDS function in the R package Vegan v2.5-6 (Oksanen et al. 2020).

Morphological Characterization of the Contact Zone. To statistically determine if there are more hybrids than parental morphs, we performed a corrected one-sided one-proportion z-test to compare the proportion of hybrids to a target proportion of 0.5, using the function ‘prop.test()’ in R v4.0.3 (R Core Team 2020). The null hypothesis is that the proportion of hybrids is not different from 0.5; the alternative hypothesis is that the proportion of hybrids is significantly larger than 0.5, indicating more males with mixed characters than males with parental characters. To test if the number of differently colored palped individuals was skewed in any direction, we performed two Χ² goodness of fit tests. The first test included all male individuals with red, white, or yellow palps, while the second test only included hybrid phenotypes – parental phenotypes were excluded. Both tests shared a null hypothesis that there are equal numbers of red, yellow, and white-palped individuals across the study site. To further explore spatial characteristics of the zone we mapped the spatial distribution of both parental forms and each hybrid type. For ease of visualization, we binned hybrid morphs into red, white, or yellow categories.

Since the contact zone that we sampled is oriented along an elevational gradient from the south (2270 m) to the north (2424 m), we tested for a relationship between palp color and mean elevation. We performed a one-way analysis of variance (ANOVA) using the function ‘aov()’ in R v4.0.3 (R Core Team 2020), with elevation as the dependent variable and palp color as the categorical independent variable with red, white, and yellow categories. We also performed a post hoc Tukey HSD test to identify if any comparison of means between specific pairs of palp color groups was statistically significant, using the ‘TukeyHSD()’ function in R v4.0.3 (R Core Team 2020).

To test for a relationship between palp color and habitat type we recorded the number of red, yellow, and white-palped individuals in two (barren and shrub) of four recorded land cover types (barren, evergreen forest, developed, and shrub), classified according to the National Land Cover Database (NLCD) 2019 Land Cover dataset (Wickham et al. 2021). We only used the barren and shrub types because there are very few samples in the evergreen forest and developed types. Additionally, while the land cover shows individuals collected on a developed road, these individuals were not actually collected on the pavement, but nearby. Individuals collected in either evergreen forest or “developed” were re-coded to barren or shrub depending on which of the two land types was closest to them. Since the observations are binary (barren or shrub), we performed a Χ² goodness of fit test with our expected counts of red = 23.5, yellow = 9, and white = 5.5. The null hypothesis was that there is no relationship between palp color and habitat type.

To test for a correlation between morphology and geographical location we performed a Mantel test using the Bray-Curtis dissimilarity matrix and a Haversine geographic distance matrix, using the R v 4.0.3 function ‘mantel()’ with 999 permutations in the Vegan v 2.5-7 package (Oksanen et al. 2020). Haversine distances were calculated using the R v 4.0.3 function ‘distm()’ in the geosphere package, version 1.5.10 (Hijmans 2019).

Molecular Data Collection and Analysis. The molecular sample included 72 adult male samples, two penultimate males, and eight females collected at the focal site. We also included a handful of individuals collected from nearby H. americanus (Gumboot Lake, ~30 kilometers distant) and H. kubai (Grass Lake, ~31 kilometers distant) populations. These non-focal H. kubai and H. americanus populations are distant enough that it is unlikely they have contributed to the current dynamic near Panther Meadows. At both non-focal locations, we only found a single species, without evidence for hybrids, but acknowledge that our sample sizes here (less than 20 total individuals per site) may not be large enough to detect rare hybrids.

Two to three legs were used for DNA extraction, performed using the Qiagen DNeasy Blood & Tissue protocol (Qiagen, Valencia, CA). DNA quality and quantity were evaluated using gel electrophoresis and a Qubit Fluorometer, respectively. We used double digest restriction-site associated DNA sequencing (ddRADseq) to gather genomic-scale data. We used the protocol described in Brelsford et al. (2017), using SbfI and MseI enzymes – a combination that increases sequencing depth while accounting for large Habronattus genome sizes (~5.6 Gb, Gregory and Shorthouse 2003). Sequencing was completed at Novogene using 150PE reads on an Illumina Hiseq 4000 platform.

Raw sequence read data were demultiplexed using Stacks v2.5.0 under default settings. After demultiplexing, the reads were processed using Stacks with the --bound_high flag set to 0.05, --min-maf flag set to 0.05, and --min-samples-overall to 90 to require 90% of individuals across the dataset to process a locus. All other parameters were left as default.

We conducted a Structure 2.3.4 (Pritchard et al. 2000) analysis under the admixture model using 2182 unlinked SNPs, including both focal and non-focal site samples. Structure was run from K 1 to 4, each replicated three times using the package Rrunstruct (Anderson 2015) and function ‘structure_runs()’ in R v4.0.3. To identify possible genetic divergence between different morph types we estimated pairwise F_ST between palp color types using all ddRAD loci (5873 SNPs), conducted in Arlequin v3.5.2.2 (Excoffier and Lischer 2010). We calculated F_ST using different palp-colored individuals as “populations” and a separate analysis using the three different palp-color hybrid groups and the two parental morphs as “populations.” Females, two intermediate-palped males, and individuals from non-focal sites were excluded from these analyses.

We tested for a correlation between Haversine geographical distances and Euclidean genetic distances using a Mantel test, with the R v 4.0.3 function ‘mantel()’ with 999 permutations in the Vegan v 2.5-7 package (Oksanen et al. 2020). Euclidean genetic distances were calculated on a matrix of 6359 unlinked SNPs using the ‘dist()’ function in the R package stats (R Core Team 2020). To test for a correlation between genetic and morphological distance, we performed an additional Mantel test using the same R function as above between the Bray-Curtis dissimilarity matrix of morphological characters and Euclidean genetic distances.