Niche conservation in copepods between ocean basins

1 –––––––––––––––––––––––––––––––––––––––– © 2021 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Dan Warren Editor-in-Chief: Miguel Araújo Accepted 30 July 2021 44: 1–12, 2021 doi: 10.1111/ecog.05690 44 1–12


Introduction
The idea that species retain their ecological niches over a wide range of temporal and spatial scales is known as 'niche conservatism' (Harvey and Pagel 1991, Peterson et al. 1999, Wiens and Graham 2005. First demonstrated empirically by Peterson et al. (1999), the prevailing opinion has been that niche conservatism is widespread across time spans ranging from years to thousands of years (Peterson 2011) and across different taxonomic groups, including plants (Petitpierre et al. 2012), mammals (Strubbe et al. 2015) and fish (Lauzeral et al. 2011). Niche conservatism is central to predicting the current and projected ecological niche of species through the use of environmental niche models (ENMs) (Pearman et al. 2010). ENM quantify the 2 relationship between a species occurrence and the underlying environment to determine their ecological niche, which can be used to predict future changes by projecting the established ENM onto future values of these variables (Guisan and Thuiller 2005). The accuracy of ENM projections relies on the assumption that the niche of a species does not change over the time interval of projection (Guisan et al. 2014) and that populations across their range have a constant niche without any significant genetic differences or adaptations. Niche divergence, where intraspecific taxa have adapted to different abiotic and/or biotic conditions, has been detected in several invasive species (Tran et al. 2015) and for microbial taxa such as coccolithophores (Zhang et al. 2014) and Escherichia coli (Deatherage et al. 2017). Identifying if intraspecific taxa exhibit niche conservatism or divergence can help us test our confidence in predicting species' responses to changing environmental and climatic conditions.
Testing the niche conservatism hypothesis is particularly important for marine plankton, where significant changes in abundance and distribution have been linked to the profound effect that climate change has had on the ocean environment, most notably with changes in ocean temperature (Beaugrand et al. 2002, Edwards and Richardson 2004, McGinty et al. 2021a). Marine zooplankton have been identified as sensitive beacons of climate change (Richardson 2008) and biogeographic shifts of taxa have been used as indicators of climate change (Chiba et al. 2018). The strong mixing and dispersal between copepod populations were thought to limit evolutionary divergence and genetic adaptation to local conditions (Norris 2000, Dawson and Hamner 2005, Provan 2013). Combined with their rapid generation times and their ability to closely follow changes to the underlying environment, changes in certain key indicator species have been used to model the effects of climate change through changes in abundance and distribution (e.g. Calanus finmarchicus - Grieve et al. 2017). There is increasing evidence that the implicit assumption of universal niche conservatism might be overly simplistic for plankton and that some species are capable of altering their ecological niche in response to climate change through evolutionary adaptation or genetic variation (Donelson et al. 2019).
Recent genetic studies have revealed that several zooplankton species with a circumglobal distribution are composed of cryptic species complexes or clades (Halbert et al. 2013, Hirai et al. 2015. Copepod species within the genus Pleuromamma exhibit distinct spatial structuring of the populations across major oceanic gyres (Goetze 2011) and show high levels of genetic variability between populations (Goetze et al. 2017). High genetic variability has also been found for the pteropod Limacina bulimoides across the Atlantic subtropical gyre (Choo et al. 2021). Significant niche differentiation was found among morphotypes of the pteropod Cuvierina spp. including those within the same clade and ocean basin (Burridge et al. 2015). Across the same oceanic scales, some species do not display the same genetic variability and high dispersal capabilities between populations as other species (e.g. C. finmarchicus -Provan et al. 2013), suggesting that oceanic barriers to gene flow might be species specific (Goetze 2005).
Field studies have shown several phytoplankton species are capable of altering their ecological niches to changing environmental conditions within a time span of a few years (e.g. Cariaco Basin, Venezuela -Irwin et al. 2015;Port Hacking, Australia -Ajani et al. 2018). Evolutionary studies have shown the adaptive potential of phytoplankton niche characteristics under climactic stressor conditions such as ocean acidification (Lohbeck et al. 2012, Collins et al. 2014. Ward et al. (2021) has shown that the interaction between dispersal and adaptive evolution of marine microbes greatly limits the connectivity between regions (Ward et al. 2021), which could apply to zooplankton also. Studies on the genetic adaptation of marine zooplankton to climate change are comparatively rare, but there is increasing evidence that rapid adaptive evolution is possible within this group also (Colin andDam 2007, Peijnenburg andGoetze 2013). Therefore, in this study, we ask whether niche conservatism is maintained between zooplankton populations despite potential genetic variation between them.
We compare the realised niches of several copepod species with known populations in multiple oceanic basins based on observations from the Continuous Plankton Recorder (CPR) surveys. Niches were compared for each study species in an n-dimensional niche space to determine if: 1) niches remain largely unchanged between populations in distant ocean basins relative to the differences in the mean environmental conditions of both areas (conserved) or 2) niches exhibited differences greater than the mean difference in environmental conditions between two areas (divergent). In this study, we define niches to be conserved if two or more populations of the same species have more similar realised niches than would be expected by chance. Conversely, we define niches to be diverged when the realised niches are significantly different between two or more populations and less similar than would be expected by chance. We quantify the extent of niche conservatism and divergence for 15 species of marine copepods with a circumglobal distribution. The significant divergence between populations will have an effect on how we quantify and predict the responses of zooplankton to climate change, as each population may not have the same sensitivity to environmental pressures (Burridge et al. 2015).

Zooplankton data
The CPR is a plankton sampler towed behind ships of opportunity at a depth of ~10 m. The plankton is captured on a calibrated roll of moving silk and divided into samples that represent approximately 10 nautical miles of towing. Raw counts are converted to semi-quantitative abundance data (number per m 3 ). The survey has collected year-round plankton samples in the North Atlantic since 1958 (Richardson et al. 2006) and has since expanded to other oceans to become one of the largest global marine sampling programmes (Fig. 1). These include surveys within the Southern Ocean since 1991 (SO-CPR -McLeod et al. 2010) and since 2007 around Australia (AusCPR - Eriksen et al. 2019). While the North Pacific CPR has also been collected routinely since 2000, data were not compared here due to the low presences of comparable species (Batten et al. 2006).
Copepods are the only zooplankton routinely identified to species level in all three regional surveys. Data were extracted from the Archive for Marine Species and Habitats Data (dassh.ac.uk) for all surveys collected between 1958 and 2020 and aggregated into monthly 1° × 1° grid cells for each survey (Fig. 1). We selected species with at least 50 presence points across all months and grid cells to allow for niche modelling, reducing sampling bias errors due to restricted or noisy distributions when there are limited data points available. A total of 15 copepod species were found to satisfy our selection criteria of at least 50 presences in two or more surveys with three species found to occur in all three areas (Supporting information). A total of 13 comparisons were made between species found in the North Atlantic and Australia, four comparisons between Australia and the Southern Ocean and five between the Southern Ocean and the North Atlantic. For brevity, we will abbreviate the comparison groups within the text as follows: SO -Southern Ocean, NAtl -North Atlantic and Aus -Australia.

Environmental data
Six environmental variables were selected for niche modelling that are known to either have a direct (i.e. controls physiological processes) or indirect (i.e. controls movement and dispersal) influence on a species. Sea surface temperature (SST, °C) and salinity are known to affect the physiological and metabolic processes of the individual and appear to be key factors driving the change in a species niche, while chlorophyll-a (Chl-a, mg m −3 ) is often an indicator for phytoplankton biomass and is used as a proxy for food availability of zooplankton (Beaugrand et al. 2013). Mixed layer depth (MLD, m) is an indicator for water column stability and has been shown to influence the niche of the copepod C. finmarchicus (Reygondeau and Beaugrand 2011). Wind stress plays a role in surface mixing and turbulence, which can have an effect on predator-prey interactions and niche structure (Barton et al. 2013). Bathymetric depth (m) is can be used to identify niches for species with an oceanic, coastal or shelf edge affinity (Helaouët and Beaugrand 2009). The environmental data were obtained from various sources including in situ measurements from ships and data buoys to create regular grids of monthly climatologies (Supporting information). These data were aggregated into the same 1° grid cells as the CPR data for each month to create 72 (6 variables × 12 months) climatological maps that were used for determining the species niche. The environmental data were linked to the CPR data by matching samples collected within the same month and within the same 1° grid.

Testing for niche conservatism
The ecological niche is a pivotal concept in contemporary ecology and evolutionary ecology. A species fundamental niche is broadly defined as the optimal range of a set of biotic and abiotic conditions that allows for a species to grow or maintain a stable population (Hutchinson 1957). The realised niche, which can be seen as a subset of the fundamental niche, describes the niche that is occupied by the species after accounting for pressures such as competition and predation. Comparing the realised niches of two populations of the same species can be confounded due to the spatial autocorrelation in environmental data (Buermann et al. 2008). Niches may differ by virtue of the fact that the available environmental conditions do not overlap. This can be addressed by first establishing a baseline or a null model that quantifies how far a niche would be expected to differ by chance based on the environmental conditions in both areas (Warren et al. 2008, McCormack et al. 2010. We judge niche conservation or divergence by comparing the difference in the realised niche for the two populations of the same species against the difference in background conditions. Thus, niche conservation is indicated when the distance between the realised niches for each population is less than the null model. Conversely, niches divergence occurs when niches differ significantly from each other and the distance between them is greater than the null model. We followed the procedure outlined by McCormack et al. (2010) to compare the niches between populations of closely related species using two methods. The first method examined the multivariate niche for each species using projections of the realised niche onto principal components (PCA method). The second approach allowed us to synthesise the effects of all variables together, using weighted combinations of the environmental variables with the most variance (ENM method; Warren et al. 2008). For each species, we extracted the monthly climatological conditions where species are present (for species with n > 1000, a random subset of 1000 is selected) and a random subset of 1000 background conditions in each ocean basin. We determined that there was evidence of niche divergence between populations if at least one of the first three PC axes and one of the null-model distributions of the ENM analysis were divergent.

PCA method
We used the first three principal components from a principal component analysis (PCA) ) of the six environmental variables to define the species niche and background conditions from the two regions being compared. A principal component is defined by a linear combination of each variable (the loadings). These linear combinations were used to compute principal component scores of the realised niche for a species and the background conditions in the two regions. The null model is defined using 95% confidence intervals (CI) for the difference in mean background principal component scores which were determined from 1000 jackknife replicates (Fig. 2). If the niche distance is less than the null model, then the niche was judged to be conserved while niches greater than the null model were found to be divergent. Differences were judged on significance using a t-test (p < 0.05).

Environmental niche model (ENM) method
Our second method quantified the differences in realised niche from the background and presence data for each pairwise comparison using Maximum Entropy (MaxEnt) models. The presence locations and environmental data were used to construct an ENM using the MaxEnt framework (Phillips et al. 2006). MaxEnt was used to estimate each population's distribution across the environmental gradients for each species. MaxEnt estimates the logistic probability of a species being present by comparing the environmental conditions where it is found with a subsample of the available environmental conditions. For each population, the model is run on a training dataset containing 75% of the available data using 100 bootstrap resampling runs with all changeable parameters standardised across all species. Threshold responses were disabled as these increase the chances of overfitting. Models were evaluated for accuracy using the remaining 25% of the data using the area under the curve to calculate the rate of true and false positives. To reduce sampling bias caused by restricted or noisy distributions, only taxa with a minimum of 50 present samples from each location were modelled. MaxEnt provides a percentage estimate of each variable's importance in Figure 2. (a) PCA method -The density of points along PC1 for two hypothetical species showing evidence of a diverged and conserved niche. For the diverged species, the mean distance between population A (yellow) and population B (green) is greater than the mean difference in background conditions (black arrow). For the conserved species, the mean distance between population A (yellow) and population 6 defining the realised niche of each population. Average variable importance was calculated by averaging each variable's importance across all population models within each survey.
The niche overlap of a species was calculated from the ENM of each ocean basin using the Schoener's D metric (Schoener 1968), which varies from 0 (no niche overlap) to 1 (perfect niche overlap). To separate the effect of different background conditions on the level of niche overlap (D), a null distribution (H 0 ) was generated for each population by calculating the differences between 100 ENMs generated using the presence data of one population and 100 ENMs generated from random background samples from the other population (Fig. 2). We assessed niche conservation by comparing D with the central 95% of the distribution of each H 0 . If D > H 0 , the niche overlap was judged to be greater than what would be obtained by chance and the niche was classified as conserved. Conversely if D > H 0 , the niche overlap is less than what would be obtained by chance and was classified as divergent. We determine that there was evidence of niche divergence if at least one of the first three PC axes and one of the ENM null distributions were divergent.

Evidence of divergence
A total of 10 of the 21 populations show evidence of niche divergence across the three ocean basins with the remaining populations showing evidence of either a neutral niche or niche conservatism (Table 1). The divergent populations belonged to 7 of the 15 copepod species compared. Divergent niches were found for four species belonging to the genus Pleuromamma (P. piseki, P. robusta, P. gracilis, P. borealis), two belonging to the genus Candacia (C. bipinnata, C. ethiopica) and the remaining divergent populations were species of Subeucalanus crassus. The divergent niches mainly occurred between populations in the North Atlantic and Australia/Southern Ocean (Table 1). P. piseki was the only species present in all three ocean basins that were found to be divergent in all subsequent pairwise comparisons.

PCA method
Niches were conserved for three of four species within the SO -Aus comparison group with P. piseki showing evidence of niche divergence along PC1 (Fig. 2). Six of 12 pairwise species comparisons made between the NAtl -Aus regions showed evidence of niche divergence. Both Candacia spp. show evidence of divergence along PC2 and PC3 (Supporting information) while for three of four Pleuromamma spp. (borealis, gracilis and piseki), there was evidence of divergence along the PC1 axis (Fig. 2) with a subset of these species showing further evidence of divergence across the PC2 and PC3 (Supporting information) axes. Subeucalanus crassus displayed evidence of niche divergence along PC1 (Fig. 2) and PC3 (Supporting information). For the SO -NAtl pairwise Table 1. The collated results of the PCA and ENM methods testing for niche conservatism are shown in Fig. 3, 4. Each cell is coloured to indicate if the pairwise test of niche similarity showed that the niches were conserved (orange), neutral (grey) or diverged (blue). The final column indicates if there is evidence of divergence for the comparison (Div.) based on the criteria of at least 1 PC axes and 1 null-model distribution is diverged. comparisons, we found that two of five species showed evidence of divergence. These species (P. piseki and P. borealis) were two of three species found in all three regions, and for all comparisons, P. piseki was found to be divergent while P. borealis was divergent between NAtl -Aus and conserved between the SO and Aus (Fig. 2).

ENM method
The variable importance for each of the six environmental variables was averaged across all populations within each basin (Table 2). SST was the most important variable in both the NAtl (41%) and SO (48%) and its importance was found to be more than double that of the next most important variable. Secondary variables that were important for ENM in the SO and NAtl include MLD, bathymetry and chl-a. For Aus populations, both wind stress and SST were equally important (~22%) with a much more even distribution of each variable's importance across the different ENM. Conserved niches were found for three of four population comparisons between the SO and Aus with populations of P. piseki found to be the sole species with a divergent niche (Fig. 3). A total of 6 of 12 populations were found to be divergent between the NAtl and Aus. These populations belonged to species within the genus Candacia, Pleuromamma and Subeucalanus (Fig. 3). A total of three of six populations were found to be divergent between NAtl and SO, which all belong to the genus Pleuromamma. P. piseki was found in all three areas and displayed significant niche divergence between all three populations (Fig. 3).

Discussion
The long-standing opinion had previously been that zooplankton and, in particular, copepods exhibit strong niche conservatism across their range (Beaugrand et al. 2013, Figure 3. PCA method -The difference in the distance along PC1 for each species pairwise comparison and the background conditions of both areas. The distance between each pair of background conditions is zero standardised and the 95% confidence interval (shaded grey area). The distance between each pairwise niche for each species is graphed in relation to the standardised background (background -niche) so that negative values indicate conserved niches and positive values beyond the 95% confidence interval are diverged. The colour of each filled circle indicates which samples are being compared, NAtl-SO (yellow), NAtl-Aus (blue), Aus-SO (red). For PC2 and PC3 (Supporting information).
8 Hinder et al. 2014). For copepods, increasing evidence suggests that genetic variability might be much greater than anticipated between separated populations thought to be the same. However, this does not preclude niche conservatism as genetic diversity does not automatically mean a change in niche (Burridge et al. 2015). Our study provides the first large-scale test of niche conservatism between populations of copepods separated across ocean basins. We examined the pairwise differences in the realised niche of 15 copepod species populations found in at least two of the three databases collected in the North Atlantic, Southern Ocean and Australian waters. A total of seven species did not conform to the niche conservatism hypothesis, showing evidence of niche divergence between their populations. Of the seven divergent species, three belong to the genus Pleuromamma, a group typically associated with oceanic warm water assemblages with a circumglobal distribution (Razouls et al. 2015). Because of their ubiquity, they have been the focus of numerous population genetic studies exploring the level of genetic and molecular diversity within this group. High genetic diversity was observed for P. piseki, P. gracilis (Halbert et al. 2013) and P. xiphias (Goetze et al. 2017) with phylogenetic analyses revealing distinct genetic lineages with diverse biogeographic distributions. Up to 18 strongly supported sub-populations were found within these species alone, making this group the most genetically diverse oceanic group of copepods known. The populations also appear to be correlated to the underlying oceanographic conditions with clade transitions associated with known oceanographic fronts or gyre systems which may act as barriers that limit dispersal or gene flow (Goetze 2011, González et al. 2020. The strongest of these dispersal barriers appears to occur within the equatorial region of the Atlantic where there is a significant genetic break for the species P. xiphias that occurs across a short distance between the southern and northern subtropical gyres (Goetze et al. 2017). The equatorial Atlantic has also been identified as a strong ecological barrier for other zooplankton groups such as pteropods (Burridge et al. 2015, Choo et al. 2021. The species Candacia ethiopica and C. bipinnata were also found to have divergent niches between the North Atlantic and Southern Ocean datasets. Like Pleuromamma, these species are also found as part of the oceanic zooplankton assemblage and in higher abundances in the subtropical and lower latitude ranges (Helaouët and Beaugrand 2009). However, evidence of genetic or morphological variation within these groups is limited. Historical observations of the species C. pachydachtyla, a congener found in similar subtropical assemblages, showed significant morphological variation between individuals collected in the Atlantic and Indian oceans (Jones 1966). Finally, significant genetic differences between individuals of S. crassus have been found in the Atlantic and Indian Ocean (Goetze and Ohman 2010), suggesting limited gene flow between the North Atlantic and Southern Ocean for this species.
For the remaining species, we find that their niches were largely conserved across the three regions studied. Information on genetic variability for these species is lacking. What we do find is that species appear to show varying levels of genetic variation within populations. For Mecynocera clausi mitochondrial sequences show high similarity between the tropical North Atlantic and the Indian Ocean (Blanco-Bercial et al. 2014), although this is limited to observations from a single location in each ocean basin. In contrast for Nannocalanus minor, there is high genetic variability within species (Hirai et al. 2013). The observation that genetic variation between populations occurs for species with conserved and diverged niches suggests that genetic responses to dispersal limitations are species specific. It is likely that niche divergence is a consequence of several interacting forces, including evolutionary and demographic history of species (Chivers et al. 2017).
The ability for a species to adapt to changing thermal regimes is still relatively unknown due to a lack of spatial, temporal and phylogenetic coverage with much of the information derived from a few near-shore or intertidal copepods (Sasaki and Dam 2021). Evidence of thermal adaptation between populations of the intertidal copepod Tigriopis californicus has been shown to occur across their latitudinal range (Pereira et al. 2017) with adaptive differentiation found to occur at small scales < 6 km (Leong et al. 2018). Intertidal species are subjected to very different evolutionary and environmental pressures than open ocean copepods with much smaller population sizes living near their thermal tolerance limits for marine copepods; the data show a positive relationship between temperature and thermal tolerance, suggesting the possibility of thermal adaptation, though this might also come at the expense of a species acclimation capacity (Sasaki and Dam 2021). A field-based study of the copepods C. finmarchicus and C. helgolandicus in the North Atlantic showed little evidence of thermal adaptation over the last 50 years (Hinder et al. 2014), suggesting that a species thermal adaptation ability might be impacted by other biotic and abiotic pressures. For Pleuromamma spp. divergence in the multivariate niche appears to be driven not just by thermal adaptation as changes across salinity and chlorophyll-a gradients appear to be just as strong. Copepods demonstrate high evolutionary potential, but the ability for a species to adapt or demonstrate niche divergence is species specific, and the reasons for this are yet to be understood. To advance our knowledge of the potential drivers of niche adaptation and gene flow, it is important that we study a greater diversity of species from a wide variety of ecological niches and varying traits (Dam 2013, Barton et al. 2016. Even closely related species have been shown to respond differently across the same environmental ranges (Blanco-Bercial et al. 2011, Chen andHare 2011). Varied species responses to dispersal and the potential presence of semiisolated genetic sub-populations have important implications for the future responses of copepod communities under changing climatic conditions (Peijnenburg and Goetze 2013). Without a full understanding of how niche divergence develops over spatial and temporal scales, the ability to use zooplankton species as climate change indicators becomes problematic. In the North Atlantic, studies involving the CPR have demonstrated that there is a prevailing northward movement in copepod communities due to ocean warming (Beaugrand et al. 2002, Villarino et al. 2015, Chust et al. 2014, McGinty et al. 2021a. Without an understanding of the population dynamics for the species, several sub-populations that respond differently to changes in the local environmental conditions may be analysed as one unit. Evidence that sub-populations of several copepod species respond differently to changing ocean temperatures shows that regional-scale variability in species responses does occur (McGinty et al. 2011).
An important consideration for understanding niche conservatism is whether the methods and niche predictors used are sufficient to determine if there is a significant difference in the niches between populations of each species. ENM approaches will be driven by the variables that are the most important in structuring the niches of both populations, which may overlook several environmental gradients that are important in driving niche divergence (Liu et al. 2020). The PCA approach addresses these issues by reducing the environmental niche to one or two important axes of niche variation (Broennimann et al. 2012) and, combined, helps to gain a fuller picture of the important ecological differences between two populations (McCormack et al. 2010). We found that many of the important variables in building the ENM were also important correlates with the major PCA axes, suggesting that the most important sources of niche variation for zooplankton were included in our analysis.
The CPR data is a unique dataset that is unrivalled in temporal and spatial scope that opportunistically samples plankton in the upper surface waters. Since only surface water is sampled, we lack detailed information on the diel vertical migration (DVM), where copepods perform daily movements between surface and deeper waters (Hays 2003). DVM behaviour is dependent on a species' size with larger species tending to undergo larger migrations (Brun et al. 2019). DVM behaviour is not fixed for a particular species as it may be adjusted according to a number of local factors such as ocean temperature. Of the species that show niche divergence, the DVM behaviour ranges from strong for Pleuromamma spp. and weak for Candacia spp. and S. crassus (Benedetti et al. 2018). Regardless, we still address the sampling inefficiencies of the CPR by using only presence data pooled across space and time (including different times of day), which will increase the probability of detecting a species as repeated sampling of the CPR in an area is randomly distributed over the course of a day.
Ecological niche models employ a correlative approach to relate the observed spatial patterns of a species with the underlying environment. While it is an invaluable tool for modelling current spatial patterns as well as future projections under climate change these approaches fail to fully capture information on the adaptive evolutionary potential of subpopulations for the species. The inclusion of phylogenetic information into ecological niche models is a burgeoning area of development. Lineage information below and above the species level can be used to produce improved niche estimates through partial pooling or splitting of similar species and their populations which could account for niche differentiation and conservatism at different levels (Smith et al. 2019). Such an analysis could be facilitated through the use of molecular markers in determining differences between species with minor morphological differences despite differences in their ecology (Choquet et al. 2021). Given the potential for both phytoplankton (Listmann et al. 2016) and zooplankton (Peijnenburg and Goetze 2013) to evolve quite rapidly, future studies should incorporate phylogenetic information into niche model analyses for these groups.