Predicting future distributions of lanternfish, a significant ecological resource within the Southern Ocean

Lanternfish (Myctophidae) are one of the most abundant and ecologically important families of pelagic teleosts, yet how these species will respond to climate change is unclear, especially within polar regions. The aim of this study was to predict the impact of climate change on the distribution of Southern Ocean lanternfish and to relate these predicted responses to species traits.


| INTRODUC TI ON
Despite being largely isolated from the world's human population, the waters of the Southern Ocean (hereafter defined as waters south of the sub-Tropical Front ca. 40-45°S) have experienced anthropogenic pressure as the exploitation of seals and whales began over 200 years ago (Laws, 1977). Within the last 50 years, there have been observations of physical environmental change including regional changes in sea ice extent (Curran, Van Ommen, Morgan, Phillips, & Palmer, 2003;de la Mare, 2009), increased acidification (Turner et al., 2014), freshening of Antarctic bottom water (Rintoul, 2007) and poleward shifts in the position of the Antarctic Circumpolar Current (Sokolov & Rintoul, 2009). These observations have been linked to anthropogenic climate change (Gillett et al., 2008) and are projected to continue in the coming decades (Bopp et al., 2013).
Globally, climate change is triggering a host of marine biological responses including poleward range shifts as species track optimal environmental conditions (Doney et al., 2012;Poloczanska et al., 2013). Polar ecosystems are thought to be particularly vulnerable to change due to the dominance of relatively short food chains (Murphy et al., 2007), the lower acclimation capacity of fauna (Peck, Morley, Richard, & Clark, 2014) and the predicted threat of invasion from subtropical species shifting distributions towards sub-Antarctic latitudes (Cheung et al., 2009;Jones & Cheung, 2015;Molinos et al., 2016). Though the underlying mechanisms determining whether a species range may shift, contract or expand are not well understood, a species response to change will depend upon their physiological sensitivity (especially thermal tolerance), resilience (e.g., ability to disperse) and exposure to climate warming (Constable et al., 2014;Day, Stuart-Smith, Edgar, & Bates, 2018;Sunday et al., 2015). Specific traits including body size (Daufresne, Lengfellner, & Sommer, 2009) and latitudinal range (Sunday et al., 2015) may also affect the direction and magnitude of species responses.
While some iconic species such as penguins and Antarctic krill have been monitored over the last few decades, making species-specific vulnerability assessments for Southern Ocean fauna has been hindered by difficulties in obtaining spatial and temporal coverage of both species records and environmental data (Constable et al., 2014).
Thus, predicting species' responses under future climate conditions, rather than directly monitoring them, can help determine how resilient or vulnerable polar species may be. For example, the modelling undertaken by Cheung, Lam, and Pauly (2008) has suggested that suitable habitat of benthopelagic Antarctic toothfish, Dissostichus mawsoni, will rapidly diminish within three decades, while changes in the distribution of a key phytoplankton species, Fragilariopsis kerguelensis, are predicted to be minimal even under a high emission scenario (Pinkernell & Beszteri, 2014). Byrne, Gall, Wolfe, and Aguera (2016) also predicted favourable conditions for the invasive Arctic seastar, Asterias amurensis, to expand its introduced range in Southern Australia to waters surrounding sub-Antarctic and Antarctic islands. This diverse set of outcomes warrants further investigation into the impact of climate change on other important species groups and into the possible mechanisms driving these different responses.
There has been growing appreciation for the ecological importance of lanternfish (Family Myctophidae), a species-rich and abundant group of fishes that inhabit the mesopelagic zone of the open ocean (200-1,000 m). In the Southern Ocean, they are the most successful pelagic fish species in terms of diversity, biomass and abundance, with over 60 species recorded below the sub-Tropical Front (Duhamel et al., 2014) and an estimated biomass of between 212 and 396 million tonnes (Lubimova, Shust, & Popkov, 1987), although this biomass could be an order of magnitude higher (Kaartvedt, Staby, & Aksnes, 2012). Their numerical dominance is associated with a vital role in ecosystem functioning, particularly as a trophic link between primary consumers (e.g., copepods and euphausiids) and megafauna including flighted seabirds, penguins and pinnipeds (Cherel, Ducatez, Fontaine, Richard, & Guinet, 2008;Saunders et al., 2015). Importantly, in years when krill are scarce, lanternfish have a key role in an alternative trophic pathway which provides a buffer to the Antarctic food web under environmental change and may maintain ecosystem stability in the long term (Murphy et al., 2007). Additionally, most lanternfish species undergo diurnal vertical migrations, moving between deeper depths during the day and shallow waters at night. Globally, this daily behaviour has a key role in the export of carbon from the euphotic zone (Belcher, Saunders & Tarling, 2019). Given their significant ecological importance, any redistribution or loss of these species would likely have consequences for the foraging success of predators, biogeochemical cycling and wider implications for ecosystem functioning (Constable et al., 2014).
Here, we examine the impact of projected climate change on 10 biomass-dominant lanternfish species in the Southern Ocean.
We use ecological niche models (ENMs) that account both for surface and deep water environmental conditions to estimate distributions under current conditions and for projected distributions under multiple future scenarios and time periods. We then relate these responses to three different species traits: thermal tolerance range, latitudinal preference and body size. Recent investigations into this lanternfish community have found that complex macroecological patterns governing their distribution and size structure are likely to affect species responses to change. Saunders and Tarling (2018) show that the majority of Southern Ocean myctophids follow Bergmann's rule, with intraspecific and interspecific body size increasing with increasing latitude and decreasing temperature. Moreover, all but two of the species studied by Saunders, Collins, Stowasser, and Tarling (2017) are suspected to spawn and recruit in regions to the north of the Southern Ocean, with only individuals from older, mature age classes reaching higher latitudes as expatriates. In the light of these findings, we anticipate that small-and large-bodied species will show different responses to climate change as body size is likely to affect latitudinal distribution patterns and, in turn, the environmental conditions experienced now and in the future.

| Species occurrence records
Occurrence records of species from the family Myctophidae living between 35 and 75˚S were downloaded from the Global Biodiversity Information Facility (GBIF). The 10 species with the highest number of records were retained for analyses, these being: Electrona antarctica, Electrona carlsbergi, Gymnoscopelus bolini, Gymnoscopelus braueri, Gymnoscopelus fraseri, Gymnoscopelus nicholsi, Gymnoscopelus opisthopterus, Krefftichthys anderssoni, Protomyctophum bolini and Protomyctophum tenisoni. All occurrence records were cleaned for unreliable data including duplicated records, records with identical latitude and longitude, and records with a latitude and longitude corresponding to a terrestrial location.
All such occurrence records were removed from the dataset.
Only occurrence records from 1960 onwards were used to correspond with the baseline period of environmental predictors. Due to a sampling bias towards the Austral spring and summer seasons (October-March), only records from these months were kept for analyses. Data from both seasons remained pooled to retain highest statistical power because there is little evidence of these species undertaking distribution shifts between these seasons Saunders et al., 2017).

| Environmental predictors
Seven environmental predictors were selected on which to build ecological niche models under present-day conditions. These were sea surface temperature (SST), sea surface salinity, temperature at 200 m, salinity at 200 m, primary productivity, dissolved oxygen and bathymetry. These were chosen based on their physiological importance for marine ectotherms and previous results demonstrating their significance as determinants of marine species distributions (Duhamel et al., 2014;Koubbi et al., 2011;Loots, Koubbi, & Duhamel, 2007). Climatological means for temperature, oxygen and salinity predictors were extracted from the World Ocean Atlas 2013 database at a resolution of 0.25° × 0.25° for the months October-March across the baseline temporal period 1956-2005(Garcia et al., 2014Locarnini et al., 2013;Zweng et al., 2013). Bathymetry data (i.e., maximum water depth) were obtained from the SRTM30 global elevation and bathymetry dataset (Becker et al., 2009) and were resampled to the same resolution as the other variables using the bi-

| Future climate data
Future climate data were derived from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) and are detailed extensively in the Intergovernmental Panel on Climate Change (IPCC) AR5 report (Flato et al., 2013 (RCP 4.5) or have continued to rise under the current trajectory to ~1,370ppm by 2,100 (RCP 8.5) (Moss et al., 2010).
Six future environmental predictors (bathymetry remained unchanged) were extracted for each ESM, RCP and time period using the NOAA/ESRL Physical Sciences Division Climate Change Web Portal (Scott et al., 2016). Only environmental conditions from October to March were used to match present-day predictors, and all future environmental data were downscaled and de-biased using the change-factor protocol described in Tabor and Williams (2010) which accounts, as much as possible, for artefactual differences between simulated and observed present-day climates. A full summary of the ESMs and processing methods used are described in Appendix S1.

| Ecological niche modelling
For each species, occurrence records and environmental predictors were fitted to the presence-only ecological niche modelling algorithm MaxEnt v. 3.3.3 (Elith et al., 2011;Phillips, Anderson, & Schapire, 2006;Phillips & Dudik, 2008). MaxEnt models the environment from a range of locations across the study region ("background sites") to discriminate against the environment at locations where species are known to be present ("presence sites"). In doing so, the model predicts the relative suitability of the environment for each species across the study region. The modelling algorithm MaxEnt was chosen for its repeatedly high performance against other ENM algorithms (Elith et al., 2006;Monk et al., 2010;Ortega-Huerta & Peterson, 2008). Moreover, MaxEnt's capacity to use presence-only data is appropriate given the high potential for errors under a presence-absence approach with mesopelagic species, due to the low sampling effort relative to the potential habitat area available, the spatial bias of sampling across the region and net avoidance behaviour being common in these species of lanternfish (Collins et al., 2008).
All ENMs were run using a 10-k cross-validation method, and 30% of occurrence data were reserved for model testing. Only linear, quadratic and hinge feature classes were selected in order to obtain a model fit that is more reliable under future conditions, following Elith, Kearney, and Phillips (2010). Ten thousand background data points were selected from within 2 decimal degrees of mesopelagic fish records across the study region. This ensures both the background and presence sites have the same spatial and environmental bias; thus, if a species occupies particular habitats within the sampled space, the model will highlight these habitats, rather than just areas that are more heavily sampled (Phillips et al., 2009). All other MaxEnt settings were kept as default.
We acknowledge high correlation between some environmental predictors (Table S2.2 in Appendix S2), and including correlated predictors can make it difficult to assess the relative importance of each due to issues of collinearity. However, there is evidence to suggest that when dealing with correlated, biologically meaningful variables, including all predictors can have a better predictive performance, in addition to a better fit, than a model parameterized using only one of the correlated predictors (Braunisch et al., 2013).
MaxEnt is particularly effective in dealing with collinearity through its iterative model fitting approach, which can consider variables independently, include nonlinear and interactions between variables and has demonstrated a robust ability to rank variables according to their importance (Braunisch et al., 2013;Phillips & Dudik, 2008). We follow the advice of Dormann et al. (2013) and confirm that the pattern and magnitude of correlation between predictors remain stable in future time periods (see Tables S2.3-10 in Appendix S2).

| Model evaluation
Model performance was evaluated using the omission rate and then by the Area Under the receiving operator Characteristic curve (AUC) score. The omission rate was determined by the proportion of presence localities that fall outside of the prediction once converted to a binary one. In this case, the binary threshold used was the minimum training presence which maintains all pixels that are predicted as being at least as suitable as those where a species' presence has been recorded (Pearson, Raxworthy, Nakamura, & Peterson, 2007).
As this threshold has an expected omission rate of zero for training localities, higher omission rates are indicative of model overfitting (Boria, Olson, Goodman, & Anderson, 2017). The AUC score is a widely used, rank-based measure of predictive accuracy that can be interpreted in the context of MaxEnt as the probability that a randomly chosen presence location is ranked higher than a randomly chosen background point (Merow, Smith, & Silander, 2013). A model with no discriminatory power will have an AUC value equal to 0.5 (no better than random), while a model with perfect fit would have an AUC value of 1.0. MaxEnt's "Jackknife test of variable importance" was also included which runs the model multiple times, each time using one of the environmental variables in isolation. This test finds the most effective single variable when predicting the distribution of the occurrence data that was set aside for testing, and gives a reliable estimate of variable importance.

| Predicting distributions under future conditions
Present-day ENMs were used to predict the future distribution of each species using 26 different climate scenarios as input in to the ENM (eight ESMs and two time periods for RCP 8.5, and five ESMs and two time periods for RCP 4.5). MaxEnt's logistic outputs, which give the conditional probability of presence between 0 and 1 for each grid cell in the study region, were then thresholded to create binary presence-absence maps of present-day and potential future distributions.
The threshold used was informed by the "maximum sum of sensitivity and specificity" (MaxSSS) threshold, as recommended by Liu, White, and Newell (2013) who compared 13 different threshold selection methods using presence-only data and concluded that MaxSSS had higher sensitivity than other methods. Furthermore, MaxSSS satisfies the three criteria necessary for sound threshold selection: objectivity, equality and discriminability (Liu, Newell, & White, 2016).
To visualize the spatial variability in future distributions based upon the different ESMs, binary future distribution maps for each species were summed together to create an index of agreement between outputs, which was repeated for each RCP and time period combination. Thus, under RCP 8.5, the summed maps have values ranging from 0 (a grid cell which is predicted to be unsuitable by all ESMs) to 8 (a grid cell which is predicted to be suitable by all ESMs used as input in to the ENM). These decrease to have a range of 0-5 under RCP 4.5 due to the more limited availability of ESM data, as noted above. Similarly, binary outputs of a species' future distribution were each subtracted from its present-day output. The resulting maps of distribution change were summed to visualize the spatial variability in the projected change for each species. Under RCP 8.5, this created an index of agreement ranging from −8 (maximum agreement of a decrease in habitat suitability across all ESMs) to +8 (maximum agreement of an increase in habitat suitability across all ESMs) and from −5 to +5 under RCP 4.5.
Before quantifying changes between present and future distributions, outputs were re-projected to the South Pole Lambert Azimuthal Equal Area projection in order to avoid potential bias of unequal cell sizes (Budic, Didenko, & Dormann, 2016). Two biogeographical metrics, centroid latitude and suitable habitat area, were To better understand the contrasting predictions among species, the mean values for each metric (ΔSHA and ΔCL) under RCP 8.5 were correlated against three species traits-species' minimum latitude of occurrence, realized thermal niche (i.e., thermal tolerance range) and maximum attained size (standard length)-using the nonparametric Spearman's rank correlation coefficient (r) to account for nonlinear relationships. Minimum latitude values were obtained from species occurrence records, after removing the 2.5 and 97.5 percentiles. We recognize that a species' latitude is determined by multiple physiological and ecological characteristics. In the context of this study, its inclusion is solely biogeographical. Occurrence locations were then matched with present-day SST to obtain realized thermal niches defined by Magnuson, Crowder, and Medvick (1979) as the temperature range in which populations persist in the wild.
Maximum attained size data were taken from values reported in Hulley (1990).

| RE SULTS
A total of 2,918 occurrence records were used in analyses (see Appendix S2, Figure S2.3). All MaxEnt models had AUC scores categorized as fair (0.7-0.8) to good (0.8-0.9) and omission rates between 0 and 0.019 (Table 1). SST and temperature at 200 m were the variables of greatest permutation importance for most ENMs, followed by primary productivity and salinity at 200 m. For seven out of 10 species, SST was the variable that gave the highest AUC score when each variable was run independently in the jackknife procedure (Table 1). Dissolved oxygen was the highest predictor for G. braueri, while salinity at 200 m was the highest predictor for G. nicholsi and G. bolini.

| Correlations with species traits
There was a significant negative relationship between the maximum attained size of a species and the minimum latitude of its occurrence (r 2 = −0.65, p = 0.04, n = 10). There was no significant relationship between a species' predicted rate of distribution shift (ΔCL; km) and body size (r 2 = 0.24, p = 0.51, n = 10) or minimum latitude of occurrence (r 2 = 0.19, p = 0.60, n = 10). Results were similarly insignificant between species' change in suitable habitat area (ΔSHA; million km 2 ) with body size (r 2 = −0.25, p = 0.49, n = 10) and minimum latitude of occurrence (r 2 = −0.07, p = 0.84, n = 10). However, ΔSHA was found to be significantly correlated to species' realized thermal niche (r 2 = 0.72, p = 0.02, n = 10; Figure 6a) showing that fish species with narrower thermal range are predicted to have reductions in areas of occurrence in the future. This correlation was similar when sea temperature at 200m depth rather than sea surface temperature was used to represent niche values (r 2 = 0.64, p = 0.04, n = 10).
Additionally, we found a significant negative correlation between a species' currently realized thermal niche and the mean latitude of its occurrence (r 2 = 0.78, p = 0.01, n = 10; Figure 6b) showing that high-latitude species (more southerly located) have narrower realized thermal niches.

| Present-day species distributions
The ecological niche models developed for this study indicate that myctophid fishes have fine-scale patterns of habitat suitability, including affinities or avoidance of shelf regions, and associations with certain water masses. Overall, these patterns are in agreement with previous biogeographical studies of these species (Hulley, 1981(Hulley, , 1990Koubbi et al., 2011;Loots et al., 2007;McGinnis, 1982). Present-day model outputs are largely consistent with previous predictions built using a presence-absence approach and a boosted regression tree (BRT) algorithm (Duhamel et al., 2014;Mormède, Irisson, & Raymond, 2014). This is reassuring, as the choice of environmental parameters and modelling algorithm can give different outputs, making interpretation difficult (Elith et al., 2006;Ortega-Huerta & Peterson, 2008). This high overlap between different modelling algorithms could be due to the presence of many frontal zones in this region which create sharp transitions in environmental conditions and are ultimately the major delimiting factors shaping bioregions and species distributions (Grant, Constable, Raymond, & Doust, 2006;Sutton et al., 2017). One noticeable exception is that our models for G. braueri, K. anderssoni and P. bolini predict distributions that extend closer to the Antarctic continent than those of Duhamel et al. (2014). This could be explained by the absence records used in their BRT approach, many of which were aggregated around the Antarctic continent, possibly exerting a stronger influence on distributions than the background data used by our MaxEnt study. Similarly, our model for G. nicholsi predicts areas of currently suitable habitat along the length of the Antarctic continental shelf, which may be expected as this species is known to be benthopelagic as adults (Hulley, 1981), and is known from habitats as far south as the western Antarctic Peninsula slope (Duhamel et al., 2014).

| Uncertainties and assumptions
The extreme complexity of the natural system results in limitations to our methodology. Our models do not take into consideration interspecific biotic interactions such as the presence of predators or prey, which can have a significant impact on the modelled range of a species (Araujo & Luoto, 2007;Bateman, VanDerWal, Williams, & Johnson, 2012). However, the large spatial scale used in this study, where climatic factors are dominant, is likely to minimize any impact of biotic interactions (Peterson et al., 2011). Our modelled distributions are based upon a simplified 2-D approach (i.e., using temperature and salinity data from multiple depths as separate environmental predictors). This may not be as full a representation of a species' niche as would otherwise be the case (Duffy & Chown, 2017), but the waters between the surface and 200 m are important habitat for these lanternfishes, which either spend the majority of their time in this depth range or migrate to shallow depths each night to feed Duhamel et al., 2014;Duhamel, Koubbi, & Ravier, 2000;Lancraft, Torres, & Hopkins, 1989;Pusch, Hulley, & Kock, 2004).
We assume that no genetic adaptation or evolutionary processes will take place by the end of our study's timeframe (2100) that some phytoplankton species have adapted to certain aspects of their environmental niche, with spatial distributions tracking changes in temperature and irradiance. Tarling, Ward, and Thorpe (2017) have also demonstrated that the distribution of the South Atlantic copepod community has remained largely unchanged over the past 80 years despite a 1°C warming in surface temperatures, which may be explained by thermal acclimation in biomass-dominant species, as well as other constraints to species distributions such as food availability.
Nevertheless, past evidence suggests that, for most species, particularly marine organisms, the dominant response to climate change is shifting distributions rather than evolutionary changes (Parmesan, Root, & Willig, 2000). We also assumed that dispersal is not limiting distributions, which is likely to be the case given genetic evidence of high connectivity of E. antarctica . Lastly, species were modelled as homogenous biomass pools, such that no changes in environmental preferences or dispersal are seen within populations, for example with age, size or density dependence (Cheung et al., 2008). These assumptions greatly simplify known population dynamics of most of these species, as spawning and recruitment occur at lower latitudes before individuals migrate to colder water (Saunders et al., 2017). Future research should go further into investigating the impacts of climate change on different age classes and their respective ranges to anticipate potential disruptions to sensitive life histories.

| Projected future distributions
We have modelled, for the first time, the impact of climate change on the distribution of several dominant Southern Ocean lanternfish species. In line with the response to ocean warming observed in many marine taxa, we found that these species will undergo poleward distributional shifts in accordance with their environmental preferences, the most important of which was temperature. Our prediction of an average range shift of 24.9 ± 13.6 km/decade under the severe emission scenario RCP 8.5, the pathway which global carbon emissions are currently tracking (Sanford, Frumhoff, Luers, & Gulledge, 2014), is at the lower end of previous estimates for marine fishes of 25-59 km per decade by the end of the century (Cheung et al., 2009;Jones & Cheung, 2015). Being mesopelagic, there is the possibility that species will move to deeper depths to compensate for increased sea surface temperature. This seems unlikely for lanternfish, however, given their dependence on productive surface waters, both during diel migration and during pelagic larval stages. Recent evidence further suggests that vertically migratory fauna that form the deep scattering layer of the oceans, to which myctophids are a large contributor, will in fact become shallower by 2100 (Proud, Cox, & Brierley, 2017), as ocean stratification and surface nutrient supply are altered by projected changes in temperature, wind stress and primary productivity.
Despite a collective poleward shift, a gain or loss of suitable habitat varied among species, with results suggesting that there will be both "winners" and "losers" to climate change. Our results indicate that E. antarctica, G. braueri, G. fraseri, G. nicholsi and G. opisthopterus have a higher probability of losing suitable habitat area by the end of the century, while E. carlsbergi, G. bolini, K. anderssoni, P. bolini and P. tenisoni have a higher probability of gaining suitable habitat area. For some of the species investigated, the direction of change was highly dependent on the climate model employed, rather than the emission scenario. Between-model uncertainty has been found to be the dominant source of climate variability in polar regions (Frölicher, Rodgers, Stock, & Cheung, 2016) and was previously found to affect predictions of myctophid species distributions more than other levels of climate uncertainty (Freer, Tarling, Collins, Partridge, & Genner, 2018). Much of the variability in the outcomes of E. antarctica, G. braueri and G. nicholsi can be contributed to the two ESMs from the NOAA Geophysical Fluid Dynamics Laboratory (GFDL), namely models GFDL-ESM2M and GFDL-ESM2G ( Figure S2.10 in Appendix S2) as these are the only ESMs to predict large areas of SST cooling south of 50° latitude. This highlights that simulating future conditions under multiple climate models is important to gain predictions that are robust and informative (Beaumont, Hughes, & Pitman, 2008;Harris et al., 2014).
To better understand the ecological mechanisms behind these opposing responses, we tested for associations between the predicted change in suitable area and species' ecological traits. Contrasting outcomes between species were not explained by differences in body size or latitudinal preferences, as may have been expected from this community which shows a trend of increasing body size with decreasing latitude (Saunders & Tarling, 2018 and this study). Instead, we find that species with a narrow thermal tolerance range are likely to lose suitable habitat, while a wide thermal tolerance range is correlated with a predicted gain in suitable habitat. This is in agreement with studies of coral reef species responses to short-term warming (Day et al., 2018) and is line with the hypothesis that broad ecological tolerances are important for range expansions within the Southern Ocean (Constable et al., 2014) and elsewhere (Sunday et al., 2015).
Importantly, species with narrow thermal niches are also those found at higher latitudes, suggesting some physiological differences between Antarctic and sub-Antarctic species. High-latitude species including E. antarctica, G. opisthopterus and G. braueri may therefore be more vulnerable to climate change as they are restricted both by their low physiological flexibility and by their biogeography. This combination renders them unable to track distributions further south due to the continental mass of Antarctica and less likely to tolerate temperatures above their extremely low optima (ca. −1 to 1°C associated with Antarctic Surface Water). This is not dissimilar to predictions for the deep-living Antarctic toothfish which was estimated to become extinct in 30 years due to its inability to move further south (Cheung et al., 2008). It also corresponds with recent predictions for Southern Ocean benthic fauna by Griffiths, Meijers, and Bracegirdle (2017) who found that endemic Antarctic species had some of the narrowest thermal ranges out of ~1,000 species south of 40° and were also the species most likely to face a future reduction in habitat.
Some of the most extreme changes in suitable area occur in sub-Antarctic species that are rarely found south of the Polar Front (Duhamel et al., 2014). Gymnoscopelus fraseri is found to have a similar thermal tolerance range to that of K. anderssoni and P. bolini yet is predicted to undergo the most severe reduction in area out of all the species analysed. The lower mean latitude of G. fraseri suggests that regions south of the Polar Front will remain unsuitable for this species by the end of the century, and it is therefore unable to expand its range poleward. Gymnoscopelus bolini and P. tenisoni have similarly low latitudes to G. fraseri but are predicted to have the largest increases in area.
These species have large thermal ranges and, unlike G. fraseri, demonstrate an ability to expand their distribution at both trailing and leading edges. Predictions of Southern Ocean macrozooplankton (Mackey et al., 2012) and benthic fauna (Griffiths et al., 2017) also found that unless species reach "gateways" of warmer water via eddy activity or shallow shelf regions, potential future ranges of sub-Antarctic taxa can be limited by steep temperature gradients across the Polar Front.

| Consequences for the Southern Ocean pelagic ecosystem
According to our findings, by the mid-21st century Antarctic waters (i.e., south of the Antarctic Polar Front) will become more favourable for smaller, sub-Antarctic species of lanternfish. This may have a combined effect of increasing the diversity of the mesopelagic fish community at high latitudes and, by increasing the proportion of small-bodied species, cause a community shift in mean body size.
In the Southern Ocean, lanternfish occupy a key trophic position and provide a major link between zooplankton and higher predators (Cherel, Fontaine, Richard, & Labat, 2010;Saunders et al., 2015).
An increase in smaller lanternfish species could therefore alter food web dynamics as most species from the genera Krefftichthys and Protomyctophum largely consume small copepods (Saunders et al., 2015a(Saunders et al., , 2015b while the larger myctophids, expected to decline in range (e.g., E. antarctica, and G. opisthopterus), have a diet dominated by euphausiids including Antarctic krill, Euphausia superba (Hulley, 1990;Saunders et al., 2014Saunders et al., , 2015a. Lanternfish often comprise up to 90% of fish preyed upon by king penguins, Aptenodytes patagonicus (Cherel & Ridoux, 1992;Olsson & North, 1997); southern elephant seals, Mirounga leonine (Daneri & Carlini, 2002); Antarctic fur seals, Arctocephalus gazella (Daneri, Carlini, Hernandez, & Harrington, 2005); and flighted seabirds (Hopkins, Ainley, Torres, & Lancraft, 1993). Thus, poleward range shifts among sub-Antarctic lanternfish may have negative consequences for predators that rely on foraging grounds north of the Polar Front (e.g., colonies on Kerguelen and Crozet islands), while those foraging south of the Polar Front (e.g., colonies on South Georgia) may benefit from the southern movement of their prey (Cristofari et al., 2018;Peron, Weimerskirch, & Bost, 2012). Predators that forage around Antarctic islands or close to pack ice target species such as E. antarctica and may be negatively affected by decreased foraging success rather than increasing foraging distance per se. A detailed investigation which integrates lanternfish ecology with predator foraging ranges at specific breeding locations would aid predictions concerning the fate of Southern Ocean predator colonies.
Overall, we have shown that despite their broad, circumpolar distributions and distance from human centres of population, the biomass-dominant species of Southern Ocean lanternfish are not immune from climate-induced impacts. Species are predicted to experience distribution shifts and changes in their suitable habitat which is likely to alter the community size structure and may have negative consequences for trophic interactions between prey and predators. We find that the direction of a species' response is dependent on the interplay between species' physiology (realized thermal niche) and biogeography (latitudinal preference), though the magnitude and direction of some species' projected responses are also determined by the climate model used to simulate future conditions. Antarctic species with restricted thermal niches and limited available habitat in which to disperse will be the most vulnerable group of Southern Ocean lanternfish in the face of climate change.

AUTH O R CO NTR I B UTI O N S
All authors contributed to the design of the study; JJF collected and analysed the data and led the writing of the manuscript; and MJG, GAT, JCP and MAC contributed critically to the draft.

DATA AVA I L A B I L I T Y
Occurrence records are available from the Global Biodiversity