Evolutionary innovations in Antarctic brittle stars linked to glacial refugia

Abstract The drivers behind evolutionary innovations such as contrasting life histories and morphological change are central questions of evolutionary biology. However, the environmental and ecological contexts linked to evolutionary innovations are generally unclear. During the Pleistocene glacial cycles, grounded ice sheets expanded across the Southern Ocean continental shelf. Limited ice‐free areas remained, and fauna were isolated from other refugial populations. Survival in Southern Ocean refugia could present opportunities for ecological adaptation and evolutionary innovation. Here, we reconstructed the phylogeographic patterns of circum‐Antarctic brittle stars Ophionotus victoriae and O. hexactis with contrasting life histories (broadcasting vs brooding) and morphology (5 vs 6 arms). We examined the evolutionary relationship between the two species using cytochrome c oxidase subunit I (COI) data. COI data suggested that O. victoriae is a single species (rather than a species complex) and is closely related to O. hexactis (a separate species). Since their recent divergence in the mid‐Pleistocene, O. victoriae and O. hexactis likely persisted differently throughout glacial maxima, in deep‐sea and Antarctic island refugia, respectively. Genetic connectivity, within and between the Antarctic continental shelf and islands, was also observed and could be linked to the Antarctic Circumpolar Current and local oceanographic regimes. Signatures of a probable seascape corridor linking connectivity between the Scotia Sea and Prydz Bay are also highlighted. We suggest that survival in Antarctic island refugia was associated with increase in arm number and a switch from broadcast spawning to brooding in O. hexactis, and propose that it could be linked to environmental changes (such as salinity) associated with intensified interglacial‐glacial cycles.


| INTRODUC TI ON
In marine invertebrates, early life-history strategy influences species dispersal potential, and this, in turn, can shape population-level gene flow and long-term evolutionary histories (Hart & Marko, 2010). Although pelagic development with planktotrophic (feeding) larvae has been suggested as the ancestral mode in most marine taxa (Strathamnn, 1985), non-pelagic direct development (brooding) has been linked to evolution under increased offspring provisioning (Wray & Raff, 1991). Brooding is commonly observed as an evolutionary transition from broadcast spawning across broad lineages (Strathamnn, 1985), but contrasting life histories (brooding and broadcast spawning) are also often reported between congeneric species in speciose clades, including echinoderms (Collin & Moran, 2018). The main drivers behind contrasting life histories are often unknown, but in some reported cases, transitions from pelagic to direct development can be linked to ecological and/or environmental changes (Boissin et al., 2011). Contrasting life histories have also been observed in congeneric species in the Antarctic and Southern Ocean (e.g., Jossart et al., 2019). If investigated, they could offer insights into variation in evolutionary processes, constrained within similar environments.
Throughout the Plio-Pleistocene (5,000,000-12,000 years ago), glacial cycles driven by climatic oscillations were significant in structuring past evolutionary histories in the terrestrial and marine realm (Hewitt, 2004;Maggs et al., 2008;Provan & Bennett, 2008). In the Northern Hemisphere, in response to ice sheet expansion and subsequent erosion of habitats, some Arctic and temperate taxa migrated to warmer, lower latitude ice-free areas for refuge, with some persisting in small-scale ice-free refugia (Maggs et al., 2008;Provan & Bennett, 2008). In the Southern Hemisphere, the continental-based Antarctic ice sheet also expanded and eroded most of the continental shelf seafloor habitats in the Southern Ocean (Clarke & Crame, 1992;Thatje et al., 2005). However, migration to lower latitudes appears improbable for Southern Ocean fauna because of the strong Antarctic Circumpolar Current (ACC) with various frontal boundaries surrounding the Antarctic continent (Rintoul et al., 2001;Thatje et al., 2005). These ocean barriers have been suggested to play an important role in isolating Southern Ocean taxa from other ocean basins since the mid-Miocene (~14 million years ago) (Crame, 2018).
Throughout the Pleistocene glacial cycles, the Southern Ocean benthic fauna are hypothesized to have either persisted in limited, isolated ice-free areas on the Antarctic continental shelf, or migrated to, and survived in, the surrounding deep sea or around Antarctic islands off the shelf Convey et al., 2009;Thatje et al., 2005).
Persistence in isolated Southern Ocean refugia has been suggested to favor non-pelagic development due to higher chances of surviving glacial periods, when the Southern Ocean experienced limited habitat availability and low primary productivity linked to permanent ice cover (Convey et al., 2009;Pearse et al., 2009;Poulin et al., 2002;Thatje et al., 2005). However, it has been suggested that pelagic development was also a successful strategy in persisting throughout glacial cycles in the Southern Ocean, and selection likely acted differently on different developmental modes throughout glacial periods (Lau et al., 2020). Survival in allopatric refugia would present unique challenges for fauna to persist within each isolated environment, as well as providing opportunities to drive evolutionary innovations and phenotypic changes (e.g., morphological variation and reproductive specialization) between isolated populations. Furthermore, the Plio-Pleistocene glacial period was also characterized by several, but rare, "warm" climate periods (between 1 and 4°C warmer than the Holocene) (Noble et al., 2020). The environmental fluctuation between extreme conditions (glacial maxima and warm interglacial) could also promote niche diversity and ecological diversification (Clarke & Crame, 1992).
Evolutionary innovations can be represented by new traits, which often opens new ecological niches where further evolutionary changes can unfold (Wagner, 2011). New traits can include trait expression and/or novel function, and these can be a broad range of behavioral, physiological, and morphological characteristics (Love, 2003;Moczek et al., 2011), whereas key innovations are a small subset of these that are invoked as underpinning evolutionary radiations. Events such as increased diversification rate and utilization of new and/or altered habitats have been suggested to be associated with evolutionary innovations (Dumont et al., 2012;Wilson et al., 2013). There is evidence indicating Southern Ocean glacial refugia could have provided opportunities for evolutionary innovations. First, persistence in Southern Ocean refugia has been suggested to have promoted allopatric diversification and subsequent speciation (i.e., the "species pumps" and the "Antarctic biodiversity pump" hypotheses) (Clarke & Crame, 1989;Crame, 1997;Willis & Whittaker, 2000). Cryptic speciation and/or lineage diversification linked to glacial cycles and/or glacial refugia survival has been suggested for many benthic taxa (Allcock et al., 2011;Baird et al., 2011;González-Wevar et al., 2019;Strugnell et al., 2012;Wilson et al., 2009). Second, Southern Ocean benthic taxa experienced repeated migrations to new and/or altered habitats as ice could be linked to environmental changes (such as salinity) associated with intensified interglacial-glacial cycles.

K E Y W O R D S
contrasting life histories, evolutionary innovation, glacial refugia, morphological innovation, population genetics sheets expanded and contracted throughout glacial-interglacial cycles (see Lau et al., 2020 for a review Bell, 1902, O. hexactis E. A. Smith, 1876, and O. taylori McKnight, 1967. Ophionotus victoriae is characterized by five arms, pelagic planktotrophic larvae (Grange et al., 2004), and a widespread distribution across the Southern Ocean at depths ranging from shallow water to 1750 m (this study; specimen IDs: WAMZ88591-WAMZ88594). Ophionotus hexactis is characterized by six arms, brooding larvae, and is mainly distributed around Antarctic islands near the APF at depths ranging from shallow water to 459 m (GBIF.org, 2019;McClintock, 1994;Turner & Dearborn, 1979). However, O. hexactis has also been collected from the Antarctic Peninsula on the Antarctic continental shelf (Hugall et al., 2016). Ophionotus taylori is characterized by five arms, and its occurrence has never been reported since the type specimen was collected from Cape Hallett, Ross Sea. Compared to O. victoriae, O. taylori possesses notably coarser and thicker scales, with other taxonomic features (including arm shape, arm spines, arm plates, oral shields) different in size and shape relative to O. victoriae (McKnight, 1967).
Ophionotus victoriae and O. hexactis are currently recognized as separate species (WoRMS Editorial Board, 2021), distinguished by the number of arms (5 versus 6) and reproductive features (oviparous vs viviparous) (Bell, 1902;Smith, 1876). However, a single five-armed O. hexactis specimen from South Georgia has also been reported to exhibit brooding behavior with six-arm juveniles (Mortensen, 1936), indicating there could be rare biological exceptions. Recent COI and exon capture data also suggest a close genetic distance between O. victoriae and O. hexactis (Galaska et al., 2017a;Hugall et al., 2016)   and O. hexactis (n = 1) from previous studies were also included in the data analysis (Galaska et al., 2017a;Hunter & Halanych, 2010) (Appendix S1 for GenBank Accession numbers). All brittle star samples (n = 935) investigated in this study were collected between the years 2004-2019, from depths of 34-1750 m during expeditions in the Southern Ocean ( Figure 1, details of sampling information and GenBank Accession numbers are presented in Appendix S1).

| Molecular sequencing
Genomic DNA of the collected samples was extracted using DNeasy Blood and Tissue Kit (Qiagen), following the manufacturer's protocol. Partial COI sequences were then amplified using genus-specific primers Op4f (5′-TAGTGACTGCCCATGCCTTC-3′) and COI_op3r

| Network reconstruction and population genetics
A median joining (MJ) haplotype network (Bandelt et al., 1999) with epsilon = 0 was constructed using PopART (Leigh & Bryant, 2015) to visualize the relationships between individual samples within O. victoriae and O. hexactis, as well as relationships between species. We have also explored a MJ network with epsilon = 10 to widen the search for unobserved sequences (i.e., hypothesized haplotypes).
However, the algorithm produced too many unnecessary hypothesized haplotypes and eventually broke down. Therefore, the MJ network with epsilon = 0 has achieved a sufficient level of exploration in reconstructing all possible shortest and least complex phylogenetic trees for discussion. A TCS haplotype network with a default connection limit of 95% was also constructed using PopART to evaluate consistency of results across different network assumptions.
For population genetic analyses, COI sequences were first grouped by species and then further divided into sample localities defined in Table 1 (Figure 1). All the sampled locations are within the APF.
Sample locations on the Antarctic continental shelf were considered "continental shelf," and islands located off the Antarctic continental shelf were considered as "Antarctic islands." Population genetic statistics including genetic diversity (nucleotide and haplotype), number of polymorphic sites, and average number of nucleotide difference were calculated for each sampling locality using Arlequin v3.5 (Excoffier & Lischer, 2010 (Villesen, 2007).
We have also examined the species boundaries and relationships between O. victoriae and O. hexactis using phylogenetic tree reconstructions and species delimitation methods. Our study only contains one genetic marker (mitochondrial COI gene), which may not contain sufficient information to diagnose species status (DeSalle et al., 2005). However, as species delimitation using COI may still be of interest to the wider research community (DeSalle & Goldstein, 2019) and is also useful in for providing hypotheses for future studies employing nuclear data, we have presented the methods and results of species delimitation in Appendix S3. Throughout this study, we view O. victoriae and O. hexactis as two taxonomically recognized separate species.

| Demographic histories
The and Fu's F S ), mismatch distributions (pairwise differences distributions), and past population size changes (Bayesian Skyline Plots; BSP). Tajima's D and Fu's F S were calculated in Arlequin to examine whether data deviated from a neutral evolution model, with significance tested by 1000 permutations. Distributions of pairwise differences to estimate parameters of demographic expansion (mismatch distribution) were calculated using the R package "adegenet" (Jombart & Ahmed, 2011) and "pegas" (Paradis, 2010)  was employed following other analyses of COI data for ophiuroids (Naughton et al., 2014;Sands et al., 2015).

| Spatial genetic variation within O. victoriae
To explore how genetic variation is spatially structured in the

| Isolation-by-environment between species
Isolation-by-environment was also investigated in MEMGENE to de- whereas seafloor temperature and salinity data were estimated based on the data value at the depth interval closest to the maximum depth available for each point. All extracted temperature and salinity data were transformed to single raster layers via triangular interpolation method in QGIS (Interpolation plug-in).
As the raster layer of surface current speed was pre-defined with a cell resolution of 16 km in Mazloff et al. (2010), the resistance surfaces analyzed in this study were interpolated (temperature and salinity) or resampled to reduce resolution (surface and seafloor elevation) to match the extent of the surface current speed layer for subsequent "mgLandscape" (MEMGENE) analysis (Appendix S4). While the interpolation of oceanic conditions may include prediction errors and deviation from true environmental conditions (especially in a heterogeneous environment) (Rellstab et al., 2015), the raster layers capture the overall dynamics of the Southern Ocean and serve as reasonable estimates for analyzing genetic-environmental association at a circumpolar scale. Collinearity between selected environmental variables was checked using a pairwise Pearson correlation analysis using the R package "Raster" (Hijmans, 2016). The resulting correlation coefficients (r between −0.51 and 0.69) were below the threshold of collinearity (r < 0.7) in ecological datasets (Dormann et al., 2013) and were therefore appropriate for subsequent environmental association analysis.
The "mgLandscape" function of "MEMGENE" was used to char-

| Population genetic metrics
Genetic diversity differed between species and sampling locations.

| Past species demography
Overall, evidence for a past population bottleneck and subsequent expansion was inferred in O. victoriae from all sample locations, as seen from a significantly negative Fu's F S value (−23.84, p = .007) and a unimodal mismatch distribution (Table 1) (Table 1).
In O. hexactis, the hypothesis of past population bottleneck and expansion was rejected due to non-significant negative neutrality tests and multimodal mismatch distribution (Table 1). A BSP also indicated an overall stable population size over time in O. hexactis ( Figure 3).

| Genetic differentiation within and between species
When analyzing genetic structure through molecular variance, Islands + Scott Island) (Figure 4a). However, the genetic structure of the continental shelf and Antarctic islands does not appear to be independent of each other as MEMGENE1 also detected genetic similarity among island localities and Prydz Bay (Figure 4a). Samples from Heard Island also showed genetic similarity from continental shelf samples (Figure 4a). MEMGENE2, the variable that explains the second strongest spatial pattern (31.8% of the adjR 2 ), further indicates relatedness between the continental shelf and Antarctic islands (Figure 4b). In particular, a strong regional structure is observed between Amundsen Sea, West Antarctic Peninsula, Scotia Sea, and Bouvet Island (Figure 4b). MEMGENE3, the variable that explains most of the remaining spatial structure in the dataset (4.78% of the adjR 2 ), demonstrated clear spatial structure connecting Scotia Sea, Heard Island, and Prydz Bay (Figure 4c).

| Isolation-by-environment between species
Isolation-by-environment analysis (via analyzing resistance sur-

| Ophionotus victoriae as a single entity based on COI data
Rather than comprising multiple cryptic species as proposed by previous studies (Galaska et al., 2017a;Hunter & Halanych, 2010), the haplotype network of O. victoriae is frequently connected, suggesting this species is one entity with a circumpolar distribution in the Southern Ocean. In the previously published studies utilizing Southern Ocean COI datasets, O. victoriae had been collected from relatively few locations in West Antarctica (Galaska et al., 2017a;Hunter & Halanych, 2010 where "individual clusters" likely represented artifacts driven by isolation-by-distance.

| Genetic relationship between O. victoriae and O. hexactis
Although current taxonomy and previous studies (Galaska et al., 2017a;Hugall et al., 2016)  of eggs) appears to be species-specific in most, but not all, brittle stars (Miller, 1998;Weber et al., 2017). Therefore, as well as the overlapping distribution of O. victoriae and O. hexactis, physiological opportunities enabling intraspecific hybridization may also exist.
Strong incomplete lineage sorting and past hybridization have also been detected among six cryptic brittle stars species Ophioderma spp. with brooding or broadcast spawning strategies (Weber et al., 2017). As we only utilized a single mitochondrial marker (COI) which is a maternally inherited marker, our data could be influenced by selection, as well as bias toward the history of mitochondrial lineages that may be incongruent with species history, and the history of maternal lineages in the event of sex-biased dispersal (Sloan et al., 2017). Overall, we highlight an interesting additional case of possible incomplete lineage sorting or hybridization between two sister taxa with contrasting morphology and life history for future multi-locus studies.

| Contrasting signals of Southern Ocean refugia between species
Evidence of deep-sea refugia in O. victoriae is demonstrated through the overall absence of population bottleneck signatures (summary statistics), combined with signs of population expansion after the LGM (BSP). Results of BSP can be confounded by the effect of population structure (Heller et al., 2013); therefore, results should be interpreted with caution. However, the overall high haplotypic diversity and the "diffused" pattern in the haplotype network also suggest O. victoriae continued to diversify during glacial periods . These patterns of past population size change and population connectivity point toward the key characteristics associated with deep-sea refugia survival (Lau et al., 2020). The deep sea was hypothesized as the only large scale, ice-free habitable area that could support a large population size and continued diversification during glacial periods; in comparison, the Antarctic continental shelf was largely covered in grounded ice (Thatje et al., 2005). While evidence of LGM grounded ice was observed around some Antarctic islands (Elephant Island, Bouvet Island, Heard Island), areas free of grounded ice were also observed around South Georgia and South Sandwich Islands during the LGM (Barnes et al., 2016;Graham et al., 2008;Hodgson et al., 2014). However, the habitable areas around Antarctic Islands are small and restricted relative to the deep sea, especially for steep, volcanic islands. Therefore, is- with brooding characteristic (Helmuth et al., 1994;Leese et al., 2010;Nikula et al., 2010).
Isolation-by-environment analysis also indicated the spatial genetic patterns between O. victoriae and O. hexactis were most associated with geographic distance and water depth, suggesting isolation-by-geographical distance and depth. Isolation-bydistance and depth are expected when populations have been stable over time, with gene flow occurring more often between spatially neighboring populations, as well as selective ecological forces and reproductive barrier between diverging populations (Wright, 1943

| Evolutionary implications of different refugial uses in the Southern Ocean
In the Southern Ocean, glacial cycles have been hypothesized to drive allopatric speciation due to populations being contained within isolated refugia on the continental shelf (i.e., the Antarctic biodiversity pump hypothesis) (Clarke & Crame, 1989, 1992Crame, 1997 (Damerau et al., 2014). Also, a Scotia Sea-Prydz Bay connection pathway has also been previously described in other Southern Ocean taxa including the asteroid Glabraster antarctica (Moore et al., 2018), the octopod Pareledone turqueti , the amphipod Eusirus giganteus (Baird et al., 2011), and the crinoid Promachocrinus phylogroup C and F (Hemery et al., 2012), highlighting a probable seascape corridor that enables gene flow between Antarctic islands and continental shelf in some Southern Ocean benthic taxa. Given that the Scotia Sea, Bouvet Island, Prydz Bay, and Heard Island are separated by a long geographical distance, there are likely unsampled regions between these areas that contribute to this proposed long-distance connectivity. The regional spatial genetic structure detected by MEMGENE2, found between the Amundsen Sea, Antarctic Peninsula, Scotia Sea, and Bouvet Island, also likely reflects the influence of local oceanographic dynamics into and beyond the Scotia Sea including the eastward flowing ACC (Maldonado et al., 2003).

| Life history and morphological innovation in O. hexactis during glacial periods
In echinoderms, brooding has often emerged under environmental stressful conditions during species selection over macroevolutionary timeframes (Lawrence & Herrera, 2000), even though this strategy requires higher maternal investment compared with pelagic larval development (Fernández et al., 2000). Previous studies have also suggested that, after a shift to hyper-oligotrophy in the Eastern Mediterranean region, the brittle star Ophioderma zibrowii with a brooding characteristic emerged from a broadcast spawning lineage (Boissin et al., 2011;Stöhr et al., 2020;Weber et al., 2019), suggesting brooding character can emerge due to historical environmental changes. Furthermore, in the Southern Ocean, brooding as a characteristic has also been hypothesized as a result of selection from nonpresent-day environmental conditions, rather than an adaptation to generic polar conditions (Pearse et al., 2009). However, an increase in arm number in echinoderms has yet to be linked to changes in past environmental conditions. Nonetheless, laboratories experiments have reported exposure to high salinity and low pH can result in arm number changes in the asteroid Echinaster sp. and in the ophiuroid Ophiothrix fragilis, respectively, under experimental conditions (Dupont et al., 2008;Watts et al., 1983). Furthermore, an increase in arm number (more than 5 arms) in ophiuroids is positively correlated to coordinated locomotion (Clark et al., 2019) and can be linked to increasingly random escape patterns (thus non-predictable escape strategies) (Wakita et al., 2020). Therefore, the six-arm innovation  (Clark et al., 2006). Since O. hexactis persisted in Antarctic island refugia in the shallow Southern Ocean, which would have been directly exposed to the prolonged, as well as intensified elements of glaciations and interglacial periods. Additionally, the rapid deglaciation at the beginning of each interglacial cycle should lead to a rapid and steep decline in salinity in the surface Southern Ocean. Ophiontous hexactis around Antarctic islands would have been directly and repeatedly exposed to deglacial meltwater after each glacial maximum. It has also been recently suggested that the surface Southern Ocean consisted of a lower level of salinity during glacial maxima (~33.4 psu relative to ~34.56 psu in the deep sea) (Hasenfratz et al., 2019). Therefore, it is also plausible that the rapid and steep decline in salinity during intensified interglacial cycles could have driven the character changes in O. hexactis, and the overall lower salinity during glacial cycles also maintained such innovations.
Alternatively, an increase in arm number and brooding strategy may not be directly linked to the proposed selective forces during the Pleistocene. The increase in arm number could be linked to ecosystem dynamics around Antarctic islands leading to enhanced coordination or could have simply arisen as a by-product of vicariance. The advantageous nature of brooding during glacial periods, when the Southern Ocean experienced limited habitat availability and low primary productivity, is more widely accepted (Convey et al., 2009;Pearse et al., 2009;Poulin et al., 2002;Thatje et al., 2005).
Nonetheless, the establishment of the morphological difference be- Each innovation leading to ecological success would have been driven by different ecological opportunities, and possibly occurred on independent occasions. Our data also demonstrate distinct genetic clusters between the Antarctic continental shelf and Antarctic islands near the Antarctic Polar Front, coinciding with the frontal boundary of the ACC in the Southern Ocean. However, genetic connectivity between the Scotia Sea and Prydz Bay was also detected, suggesting connectivity between the Antarctic islands and the Antarctic continental shelf is ongoing. Genetic structure observed between the Scotia Sea and neighboring regions (Bellingshausen Sea, Antarctic Peninsula, Weddell Sea, and Bouvet Island) also highlights the role of the ACC and Weddell gyre in structuring regional genetic patterns.

| CON CLUS ION
Finally, our work also discussed that the morphological and lifehistory innovation in O. hexactis (an increase in arm number and brooding as reproductive strategy) can be linked to selection from environmental conditions in the past, which was first proposed by Pearse et al. (2009). Further work examining genetic structure in O. victoriae and O. hexactis using nuclear data should provide a thorough understanding of the genetic relationship between the two species, and signatures of past environmental selection.

ACK N OWLED G EM ENTS
We thank the Western Australian Museum (

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
All newly generated COI sequences were deposited to GenBank at NCBI under accession numbers MZ543435-MZ543949.