\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,polish,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Global analysis of the influence of environmental variables to explain
distributions and realized thermal niche boundaries of sea snakes}



\author[1]{Carlos Patrón-Rivero}%
\author[2]{Luis Osorio-Olvera}%
\author[3]{Octavio R. Rojas-Soto}%
\author[4]{Xavier Chiappa-Carrara}%
\author[3]{Fabricio Villalobos}%
\author[1]{Kevin Lopez Reyes}%
\author[5]{Brooke L. Bessesen }%
\author[1]{Carlos Yañez-Arenas}%
\affil[1]{Universidad Nacional Autonoma de Mexico Facultad de Ciencias}%
\affil[2]{Universidad Nacional Autonoma de Mexico Instituto de Ecologia}%
\affil[3]{Instituto de Ecologia}%
\affil[4]{Universidad Nacional Autonoma de Mexico Escuela Nacional de Estudios Superiores Unidad Merida}%
\affil[5]{Ecology and Evolutionary Biology,University of Reading, Reading, UK}%


\vspace{-1em}



  
  \date{May 26, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Understanding the factors that affect species' distributions is a
central topic of ecology and biogeography. However, the majority of
research about this topic has focused on species inhabiting terrestrial
environments. At broad scales, abiotic variables consistently serve as
the primary determinants of species' distributions. In this study, we
investigate the explanatory power of different abiotic variables in
determining the distributional patterns of sea snakes on a global scale.
Additionally, as the boundaries of realized thermal niches have
significant implications for the ecology of the species and their
geographic distributions, we also evaluated the asymmetry of realized
thermal limits (i.e., differences in variances between upper and lower
limits of the realized thermal niche). We obtained global environmental
variables and occurrence data for each species across their entire known
geographic range. Using this data, we employed a correlative ecological
niche modeling procedure to analyze the influence of individual
variables in explaining species' distributions. To estimate the realized
thermal limits of each species, we extracted the mean, minimum, and
maximum values of temperature at four depths (superficial, mean benthic,
minimum benthic, and maximum benthic) for each occurrence record of the
species. We then evaluated the asymmetry of the realized thermal niche
by measuring and comparing the variances of the upper and lower limits.
Both analyses (the importance of variables and realized thermal limit
asymmetry) were performed at three taxonomic levels (sea snakes as a
lineage of marine-adapted elapids {[}true sea snakes + sea kraits{]},
subfamily, and genus) and at two spatial resolutions. Overall, we found
that temperature, phosphate, nitrate, salinity and silicate
concentrations were the most influential factors in explaining the
spatial distribution patterns of sea snakes, regardless of the taxonomic
level or spatial resolution. Similarly, we observed that the realized
thermal limits were asymmetric with higher variance in the lower limits,
and that this asymmetry decreases as the taxonomic level and spatial
resolution increased. Finally, our results align with previous findings
regarding patterns of asymmetry in realized thermal limits and the
significant influence of abiotic variables in explaining the
distribution of marine species.%
\end{abstract}%



\sloppy


\textbf{1 INTRODUCTION}

Anthropogenic activities have led to rapid global changes that are
creating unprecedented climatic conditions, which could compromise the
viability of many species (Walther et al. 2002; Walther et al. 2009).
Understanding how environmental conditions affect the survival of
species' populations and shape their distributions is fundamental for
establishing conservation strategies (Bosch et al. 2018). However,
estimating the relationships between species and their environment is a
complex task, as it involves multiple factors (Gaston and Fuller 2009;
Sexton et al. 2009; Wiens 2011). These factors can be summarized into
three categories: 1) limiting or regulatory, which control the
ecophysiological responses of the species, 2) disturbance factors, which
describe the historical modifications in the occupation of habitats by
both natural and anthropogenic processes, and 3) resource factors, which
represent the supplies necessary for the survival of organisms (Guisan
and Thuiller 2005). Some causal relationships between these factors and
species' attributes, such as abundance and distribution, have been
demonstrated through physiological experiments in laboratories and
wildlife demographic studies (Gause 1931; Hutchings and Myers 1994;
Kordas et al. 2011). However, carrying out laboratory experiments or
fieldwork to obtain demographic data has some drawbacks. First, both are
generally costly and time-consuming. Secondly, implementing
physiological experiments is not ethically feasible for many organisms
(Fraser 1999). Finally, some environmental variables are more relevant
at broad spatial scales (\textgreater{}200 km). For instance, limiting
or regulating factors associated with climate are considered to be the
most relevant in explaining the distribution of species across their
entire range (Contreras-D\selectlanguage{ngerman}íaz et al. 2022; Gillespie et al. 2008; Pearson
and Dawson 2003).

An alternative approach to studying species-environment relationships in
a spatial context is the use of correlative methods, known generically
as ecological niche modeling (ENM) and species distribution modeling
(SDM) (Bentlage et al. 2013; Pearson and Dawson 2003; Peterson and
Soberón 2012). The inputs required by correlative methods are easier to
acquire since they are generally available in online repositories. These
inputs are sets of georeferenced localities where the presence of the
species in question has been observed (some methods also require points
that represent absence) and spatial surfaces that describe the
environmental conditions of the area of interest (Randin et al. 2009).
Thus, correlative methods are extensively used to model
species-environment relationships and estimate their actual and
potential distributions, both in terrestrial (Colwell and Rangel 2010;
Eme et al. 2014) and marine ecosystems (Bosch et al. 2018; Bradie and
Leung 2017). Although correlative methods have certain limitations
(e.g., it can be difficult to discern whether their results are mere
statistical associations or causal relationships; see Wiens et al.
2009), they allow us to recognize diverse ecological patterns (Ashcroft
et al. 2011). As a result, these methods have become popular for
analyzing the importance of environmental variables as limiting factors
of distribution ranges and various ecological aspects across multiple
taxa (Bosch et al. 2018; Bradie and Leung 2017; Lee-Yaw et al. 2016).

One of the most studied limiting factors is temperature, as it has a key
relevance for many physiological and ecological processes (Bicego et al.
2007; Clarke 2003; Kordas et al. 2011). Thus, studying the thermal niche
(TN) and realized thermal niche (TNR) of species has been a major focus
of ecological research (Gvo\selectlanguage{polish}ž\selectlanguage{english}d\selectlanguage{ngerman}ík 2018; Lee-Yaw et al. 2016; Magnuson et
al. 1979; Stuart-Smith et al. 2017). According to Gvo\selectlanguage{polish}ž\selectlanguage{english}d\selectlanguage{ngerman}ík (2018), TN
refers to the range of body temperature that individuals of a species
require for the population to experience positive growth, while TNR
refers to the environmental range of temperatures that individuals of a
species are exposed to (Collin et al. 2021; Stuart-Smith et al. 2017).
Based on the concepts of TN and TNR, two hypotheses have been proposed
to explain patterns of temperature tolerance and organismal response to
global warming. The first hypothesis suggests that tropical and polar
terrestrial species have lower TN than those from temperate climates
(Pörtner 2001; Pörtner and Farrell 2008). The second hypothesis is the
realized thermal limit asymmetry, which refers to significant
differences in the variation of upper limits in a species or taxon
compared to the variation of lower limits within the same taxon, with
the variation being greater in the lower limits. This means that under
climate change scenarios, high rates of global warming can severely
impact the conservation of most species (Collin et al. 2021;
Stuart-Smith et al. 2017). Ideally, evaluating these hypotheses requires
direct estimation of the species' TN in question through physiological
population growth-mortality experiments (Dong and Somero 2009; Helmuth
et al. 2006; Sánchez-Fernández et al. 2012; Stillman and Somero 2000).
However, despite the high cost, time and ethical constrains, there are
global initiatives to analyze this patterns (Bennett et al. 2018).
Nonetheless, characterization of the TN for most taxa is still lacking.

Therefore, through the use of correlative methods, it is possible to
characterize a portion of the TN called TNR, which represents the range
of temperatures associated with the presence of a species (Descombes et
al. 2015; Gaüzère et al. 2015; Quintero and Wiens 2013). Since
correlative methods do not require experimentation and management of
organisms, they are a more feasible option for studying species TNs. As
a result, they have been commonly applied to many taxonomic groups (Huff
et al. 2005; Kuo and Sanford 2009; Sánchez-Fernández et al. 2012;
Sánchez-Fernández et al. 2016). Moreover, most TNs analysis have been
focused on species from terrestrial environments (Bennett et al. 2018),
leaving no information on many important marine groups, such as sea
snakes.

Sea snakes are an ecologically diverse and marine-adapted lineage of
elapids that includes two subfamilies, Hydrophiinae (true sea snakes)
and Laticaudinae (sea kraits). The group comprises \textgreater{}60
species widely distributed in the tropical and subtropical regions of
the Indian and Pacific Oceans (Sanders et al. 2013; Wallach et al.
2014). Sea snakes' lineage has evolved with unique physiological
adaptations to survive in marine environments, being their ability to
thermoregulate a critical adaptation to adjust their metabolic rate and
behavior in response to water changes in temperature (Heatwole et al.
2012). Recent research has also shown that temperature plays a critical
role in the diving behavior and metabolism of sea snakes, as individuals
exhibit different diving behaviors and metabolic rates in response to
changes in water temperature (Dabruzzi et al. 2012; Udyawer et al.
2016). These findings highlight the complex ways in which temperature
influences the behavior and physiology of sea snakes and underscore the
importance of understanding the thermal biology of these reptiles for
conservation and management efforts. Nevertheless, most sea snake
research has focused on their toxicity, and their phylogenetic and
evolutionary relationships (Udyawer et al. 2018). However, to date,
there is no information on the factors that determine sea snakes'
distribution, nor the asymmetry of their realized thermal limits. In
this study, we used correlative ecological niche modeling (ENM) to
evaluate the relative importance of environmental variables in
explaining the distribution of almost all sea snake species on a global
scale. We also analyzed whether there is an asymmetry in the dispersion
of their lower and upper realized thermal limits. The results of our
research provide a better understanding of the biogeographical phenomena
that determine the distribution patterns in this group and will allow us
to infer their vulnerability to future ocean warming scenarios.

\textbf{2 MATERIALS AND METHODS}

\textbf{2.1 Presence records}

We searched presence records with geographic coordinates for all sea
snakes in public online repositories. These platforms were the Global
Biodiversity Information Facility (gbif.org; DOI information associated
with data are available in Table SM1), VertNet (vertnet.org), the Ocean
Biogeographic Information System (obis.org), the Atlas of Living
Australia (ala.org.au), the Online Zoological Collections of Australian
Museums (ozcam.org.au) and the UNAM open data portal
(datosabiertos.unam.mx). We integrated the presence records obtained
from the different platforms into a single database per species. Then,
we eliminated duplicate records, uncertain data (e.g., taxonomic
mistakes, records outside the known distribution of the group), and
obvious errors (e.g., records in terrestrial environments). Next, we
eliminated occurrences that might represent sink populations through an
environmental outliers' analysis. In this process, we considered
environmental variables previously identified in the literature as
important to the group (see below). For this, we used the\emph{outliers}
function of the ``biogeo'' package (Robertson et al. 2016) in R (R Core
Team 2013), which searches for presence records that are located at a
distance above 1.5 times the interquartile range of the data.
Additionally, we employed the \emph{gridSample} function of the
``dismo'' package (Hijmans et al. 2017) to minimize data bias towards
more intensely sampled environments. These processed records comprised
the complete database for each species.

\textbf{2.2 Environmental data}

We downloaded environmental variables in raster format, commonly used in
spatial modeling of marine environments from the global repositories
Bio-Oracle (Assis et al. 2018) and MARSPEC (Sbrocco and Barber 2013) at
four depths (superficial, benthic maximum, benthic mean and benthic
minimum). Then, we classified the variables into five groups: Bsurf
(superficial Bio-Oracle; n = 10), Bma (maximum benthic Bio-Oracle; n =
8), Bme (mean benthic Bio-Oracle; n = 8), Bmi (minimum benthic
Bio-Oracle n = 8), and Msurf (superficial MARSPEC; n = 6). This protocol
was used to select the combination repository/depth that maximizes the
explanatory capacity of predictors (Table 1). To standardize the spatial
resolution of the variables, we resampled the MARSPEC variables to match
the 5 arc-minutes cell resolution of Bio-Oracle. Subsequently, we
resampled both repositories at a resolution of 10 arc-minutes with the
same method, to assess the effect of spatial resolution on our results
(Connor et al. 2019).

For each species, we defined a polygon representing a hypothesis of its
historical accessibility (area \textbf{M} ; sensu Soberon and Peterson
2005) based on the marine biogeographical provinces of the world
(Spalding et al. 2007) and ocean current information
(earth.nullschool.net/\#current/ocean/surface/currents/patterson). We
first selected as area \textbf{M} all the provinces with at least one
presence record of the species in question. Then, based on the natural
history of the species and the geographical configuration of its
occurrences, we limiting or expanding it depending on whether we
consider that the currents could serve as physical barriers or dispersal
channels. Then, we also defined a calibration area and masked the
variables by selecting only the areas within the intersection between
the presence records of the species in question and the marine
eco-regions; a finer sub-regionalization than the provinces. Finally, to
reduce collinearity and dimensionality of predictors, we eliminated
individually for each species, variables in each of the five groups that
had a pairwise Pearson correlation \textgreater{}0.8 with the
function\emph{correlation\_finder} in ``ntbox'' package (Feng et al.
2019; Mateo et al. 2013; Osorio-Olvera et al. 2020).

\textbf{2.3 Modeling ecological niches}

We use the Jackknife cross-validation procedure for improving the
accuracy and to minimize the variability of the models in taxa with
fewer than 25 localities. In this method, one presence record is
excluded randomly from modeling and the procedure is repeated as many
times as the data allow; all records are excluded just once. In each
iteration, the \emph{n} -1 data are used as model training information
and the excluded record is used for model evaluation (Shcheglovitova and
Anderson 2013). Finally, in species with more than 25 occurrence
records, we applied the \emph{get.checkerboard1} function in the
``ENMeval'' package (Muscarella et al. 2014). In this method, the
database is divided into two sets using a checkerboard pattern across
the extent of the study area. Occurrences are separated according to
their position on the board. Finally, we defined the set with more
records as training data and the set with less information as evaluation
data (Brown 2014).

To characterize the ecological niches of sea snakes we used the MaxEnt
3.4.1 algorithm (Phillips et al. 2006). A series of candidate models
were calibrated with the training database of each species via the
``kuenm'' package (Cobos\selectlanguage{english} et al. 2019). This package allows varying
parameter combinations of the MaxEnt algorithm. We tested combinations
of seven feature classes in MaxEnt (``l'', ``q'', ``p'', ``lq'',
``lqp'', ``lp'', and ``qp'', where ``l'' is linear, ``q'' is quadratic,
and ''p'' is product), eight regularization multipliers (0.10, 0.25,
0.50, 0.75, 1, 2, 3 and 4), and the five groups of environmental
variables described above (thus, the maximum possible number of
variables used to build any model is full set of ``repository/depth'' if
no correlation between all variables were found). This results in 560
candidate models per species (spp x 2 resolutions x 7 Fs x 8 RMs x 5
groups of variables). Then we selected using the test database the best
subset of models that met the following criteria hierarchically: 1)
statistical significance, where we retained those models that were
better than expected by chance, depending on the proportion of bootstrap
replicates with ratios of the partial area \textgreater{}1 (Peterson et
al. 2008), 2) predictive capacity, where we selected those models with
statistical significance, which were also able to predict at least 90\%
of the evaluation records (i.e., models with an omission rate --OR-- \selectlanguage{english}[?]
0.10), and 3) complexity, through the Akaike information criterion
corrected for small samples (\selectlanguage{greek}Δ\selectlanguage{english}AICc) we selected as best combinations
those [?] 2 units, that are also the models that present a better fit and
fewer parameters (Radosavljevic and Anderson 2014; Warren and Seifert
2011). Based on the best parameter combinations (full set of best
combinations are available in Table SM2), we built a final set of models
using the whole database (training and evaluation data) with the
bootstrap functionality of MaxEnt performing 10 replicates. In each
iteration, we divided presence records randomly into 80\% for training
and 20\% for evaluation; we define Cloglog as output format and 10,000
background points (masked by our calibration area). Finally, we
calculated the median and the range of the predicted values across the
total replicates (10 replicates \selectlanguage{ngerman}× final parameter combination) to
represent the consistency and variation in predictions.

\textbf{2.4 Importance of variables}

To analyze the relative importance of environmental variables in
explaining sea snakes' distributions, we first classified the variables
into ten groups based on their identity. For example, the group
``Temperature (Tem)'' included seven variables that represent the median
of the mean, maximum and minimum values of the sea surface and benthic
temperature (Table 1).

We used two measures to assess the relative importance of environmental
variables: the percentage of contribution (PC) and the permutation of
importance (PI) estimated by MaxEnt. PC measures the gain of a model
built with a single specific variable divided by the total gain of a
model built with all the variables. PI is obtained from randomizing the
values of each variable so that this variable is not informative and
then measures the resulting drop in the area under the curve (Phillips
2011). To identify the consistency or variability in our results, we
performed this procedure at three levels: ``lineage'', when we analyzed
the variables for all species independently of subfamily; ``subfamily'',
where we divided the results into the Hydrophiinae and Laticaudinae
subfamilies; and ``genus'', where we repeated the same analysis for each
of the six genera. We performed this analysis at both resolutions (5 and
10 arc-minutes).

We also separated the identity variables into two groups (high relative
importance and low relative importance) based on their relative
importance using the natural Jenks thresholds calculated with
the\emph{getJenksBreaks} function of the ``BAMMtools'' package (Rabosky
et al. 2014).

To determine if the results depended on spatial resolution, we evaluated
the consistency of the importance of each variable between the two
spatial resolutions for the two MaxEnt metrics (PC and PI) through a
non-parametric Mann-Whitney-Wilcoxon test through
the\emph{wilxicon.test} function in ``stats'' package (R Core Team,
2013). Overall, our approach allowed us to identify the most important
environmental variables driving the distributions of sea snakes and to
assess the consistency of these results across different metrics of
importance, taxonomic levels and spatial resolutions.

\textbf{2.5 Asymmetry of the realized thermal niche}

We extracted the mean, maximum, and minimum temperature values for the
five groups depending on the repository and depth (Bsurf, Bma, Bme, Bmi,
and Msurf) associated with each post processed presence record (see
above in ``2.1 Presence records''\textbf{)} of each species at 5 and 10
arc-minutes of spatial resolution. For Bio-Oracle, we used directly the
mean, maximum, and minimum temperature variables in each depth. In the
case of MARSPEC, we used the surface temperature of the warmest month to
represent the maximum temperature and the surface temperature of the
coldest month to represent the minimum temperature. With those values,
we conserved and grouped the minimum and maximum value from each species
for characterized the lower and upper TN\textsubscript{R} limits
respectively.

Finally, we tested the asymmetry in TN\textsubscript{R} limits
hypothesis (asymmetry in the variances with more variance in lower
limit) by testing the homogeneity of variances between the lower and
upper realized thermal limits. This process was repeated for each
variable group on each temperature measurement at the three taxonomic
levels and two spatial resolution by applying a non-parametric
Fligner-Killeen test (Fligner and Killeen 1976) (5 groups × 3
temperatures × 3 taxonomic levels × 2 resolutions). Particularly, in
genus level, we excluded the genus \emph{Laticauda} due to represent the
same result at the Laticaudinae subfamily, and \emph{Ephalophis}
and\emph{Microcephalophis} as a result of only have data for one
species.

\textbf{3 RESULTS}

We obtained 5,020 spatially unique presence records of 50 of the
\textgreater{}60 recognized sea snakes. Most of the records of the
species were distributed in the Indo-Pacific Ocean and \emph{Hydrophis
platurus} being the only species with representation in the Pacific of
the new world (density map by genus is available in Fig\_SM1.
Additionally, an interactive map can be performed with
``3\_Density\_plot.r'' code in
https://github.com/PatronRiveroC/Ssnakes).

Regarding model evaluation, we obtained a median partial AUC ratio
\textgreater{}1 for all species (i.e., all our models were better than
expected by chance). The ORs in the final models were \selectlanguage{english}[?] 0.10, but three
species when we used the variables at 5 minutes resolution and seven
species at 10 minutes and six for 10 minutes (Table 2). We show species'
environmental suitability maps for each species in Fig\_SM2-SM51.

\textbf{3.1 Importance of variables}

All results of relative importance for each taxonomic level (lineage,
subfamily and genus) are shown in Figure 1 (results by each species are
in Fig\_SM52). Summary of the relatively high important variables are
presented in Figure 2 for the PC and PI estimated by MaxEnt at 5 and 10
arc-minutes of spatial resolution.

\textbf{3.1.1 Lineage}

At a resolution of 5 arc-minutes, the most important variables for the
sea snakes' lineage were phosphates (Pho), salinity (Sal), and
temperature (Tem), which together accounted for the highest PC. However,
when the importance of variables was permuted (PI), nitrates (Nit),
silicates (Sil), and Tem emerged as the most important factors. At 10
arc-minutes for PC and PI, Sil and Tem were important variables, while
Nit was also important in PI. Finally, we only found differences between
the median of Sil when comparing both resolutions in PC (Table SM3).

\textbf{3.1.2 Subfamily}

At the subfamily level, the patterns were different. For the
Hydrophiinae subfamily, Sal and Tem were the most important variables at
5 arc-minutes for PC, while Nit, Sil and Tem were the most important for
PI. At 10 arc-minutes, Sil and Tem were the most important variables for
PC and Nit and Sil for PI. For the Laticaudinae subfamily, the variables
with highest PC at 5 arc-minutes were Nit, Pho and Sal, while the
variables with highest PI included in addition the Tem. At 10
arc-minutes, Sil and Tem were the most important variables for both PC
and PI. Finally, we did not find differences between each pair of the
same variable when comparing resolutions for both subfamilies in PC and
PI (Table SM3).

\textbf{3.1.3 Genus}

At the genus level, for \emph{Aipysurus} at 5 arc-minutes, calcites
(Cal), dissolved molecular oxygen (Doxy), iron (Iro), Pho, Sal and Tem
were the most important variables for PC, while Doxy, Iro, Pho and Tem
were the most important for PI. At 10 arc-minutes, Sil was the most
important variable for PC and Doxy, Pho and Sil for PI.
In\emph{Emydocephalus} at 5 arc-minutes for PC, Doxy, Sal, Sil and Tem
were the most important variables, while Sil and Tem were the most
important for PI. At 10 arc-minutes, currents Velocity (Cvel), Sal, Sil
and Tem were the most important variables for PC and the same variables
omitting Sal for PI. For \emph{Ephalophis} at the 5 arc-minutes for both
PC and PI, Iro, Sal and Tem were the most important variables. While at
10 arc-minutes, Nit and Sal were important for PC and only Nit for PI.

For \emph{Hydrophis} and \emph{Laticauda} , the results at the genus
level were similar to those obtained at the subfamily level due to
contain same all species in case of Laticaudinae (the reason for we are
not presented the same results here) and almost all for Hydrophiinae.
Specifically, for \emph{Hydrophis} at 5 arc-minutes for PC, Nit, Pho,
Sal and Tem were the most important variables, while Nit and Tem were
the most important for PI. At 10 arc-minutes, Nit, Sil and Tem were the
most important variables for both PC and PI, and additionally, Pho was
important only in PC.

Finally, \emph{Microcephalophis} was influenced by Sal and Tem at 5
arc-minutes for PC, while Tem was the only important variable for PI. At
10 arc-minutes, Tem was the most important variable for both PC and PI.
We also did not find any statistical difference between comparison same
variable in both resolutions, importance metric in any genus (Table S3).

\textbf{3.2 Variation in the TN\textsubscript{R} limits}

At the lineage level, our results reveal significant differences in the
variance of the realized lower thermal limit compared to the realized
upper thermal limit across all temperature measurements (mean, maximum,
and minimum), regardless of the group of variables used (Bma, Bme, Bmi,
Bsurf, and Msurf) or cell resolution (Fig 3A-B; Table SM4). The same
pattern was observed in the Hydrophiinae subfamily (Fig 3C-D; Table
SM4), whereas for the Laticaudinae subfamily, we found this pattern only
in the minimum temperature measure at both resolutions (Fig 3E-F; Table
SM4).

At the genus level, we observed higher variances in the lower limits
than the upper limits for \emph{Aipysurus} and \emph{Hydrophis} and all
temperature measures, repository and depth, and spatial resolution
(Figure 3G-H; Table SM4). However, we did not find significant
differences between the lower and upper limits in any temperature,
repository and depth, and spatial resolution for the\emph{Emydocephalus}
genus (Figure 3I-J; Table SM4).

\textbf{4 DISCUSSION}

\textbf{4.1 Environmental variables}

Our study confirms the importance of environmental variables in shaping
the distribution of sea snakes, which is consistent with previous
studies (Bessesen et al. 2023; Udyawer et al. 2014). Environmental
variables play a crucial role in determining the habitat suitability and
area of occupancy of marine species, including sea snakes, and
prioritizing search efforts based on habitat suitability models can help
locate previously unknown populations of endangered marine reptiles
(Udyawer et al. 2020). This aligns with our findings and emphasizes the
importance of our research in generating crucial information that could
guide future conservation and management strategies for protecting a
highly vulnerable and ecologically important group, which is sea snakes.

Our findings are consistent with a meta-analysis on the importance of
variables used in MaxEnt species distribution models, which found that
the most important factors in determining the distribution of marine
organisms were temperature, bathymetry, and precipitation (Bradie and
Leung 2017). Additionally, our study highlights the importance of
nutrient availability in shaping the distribution of marine species.
Those environmental variables were not important in the meta-analysis.
The study by Bosch et al. (2018) found that temperature, bathymetry, and
distance to shore were the most important predictors for marine species
distribution, while salinity and nutrients had less influence. Although
our study's results on the importance of salinity and nutrients in
determining sea snakes' distribution appear to conflict with previous
analysis, it is worth noting that our study specifically focused on sea
snakes, which may have different environmental preferences than other
marine species. Our study's results on the importance of phosphates,
salinity, and temperature in determining the habitat suitability and
area of occupancy of sea snakes are consistent with the findings of
Bessesen et al. (2023) on \emph{H. platurus} and Udyawer et al. (2014)
on sea snakes in the Great Barrier Reef Marine Park.

The combination of climate and nutrient variables are known to have
implications for the morphology, diving behavior, feeding habits of sea
snakes and also coral reefs health (Brischoux et al. 2012; Gherghel et
al. 2018; Heatwole et al. 2012; Sherratt et al. 2018; Veron 2000, 2008).
Temperature, along with salinity, plays a crucial role in the physiology
of sea snakes, affecting their metabolic rates, oxygen uptake, and blood
chemistry (Heatwole et al. 2012; Udyawer et al. 2016). Studies have also
demonstrated that metabolic thermal sensitivity optimizes sea krait
amphibious physiology, enabling them to tolerate fluctuations in
temperature (Dabruzzi et al. 2012; Heatwole et al. 2012). Furthermore,
trophic specialization drives morphological evolution in sea snakes,
while the dependence of dive behaviors and metabolism on temperature has
also been observed (Dabruzzi et al. 2012; Sherratt et al. 2018; Udyawer
et al. 2016). However, plasticity in thermal tolerance has limited
potential to buffer ectotherms from global warming, which may lead to
mosaic patterns of thermal stress with implications for climate change
(Gunderson and Stillman 2015; Helmuth et al. 2006). Phosphates play a
significant role in the distribution of sea snakes as they are essential
for the growth and development of the coral reefs where sea snakes are
commonly found (Veron 2000). Although the importance of nitrogen and
silicates in the distribution of sea snakes is less clear, they are
likely to play a role in the availability of prey species that sea
snakes feed on (Gherghel et al. 2018). Nitrogen is an essential nutrient
for many marine organisms, and its availability can influence the
abundance and distribution of prey species (Brischoux et al. 2012).
Silicates are important in the growth of diatoms, which are a primary
food source for many marine organisms, including some species of sea
snakes (Sherratt et al. 2018).

The sea snakes' distribution may have influenced their diet and feeding
behavior, as species that occur in areas with higher nutrient
availability may have evolved to specialize on prey that require higher
energy requirements, such as fish, while species occurring in areas with
lower nutrient availability may have evolved to specialize on prey that
require less energy, such as crustaceans (Gherghel et al. 2018; Sherratt
et al. 2018).

Our study utilized correlative models to determine the relationships
between sea snakes and environmental variables, but caution must be
exercised in interpreting the results since correlative models cannot
establish causality. Some correlations may reflect direct effects of
environmental variables, while others may represent indirect
relationships. We used MaxEnt, which can handle collinearity among
explanatory variables; however, collinearity may still pose challenges
when estimating the importance of explanatory variables (Smith and
Santos 2020). The programmers of MaxEnt incorporated two different
metrics (PC and PI) to evaluate the contribution of predictors to the
model (Phillips 2011), although there is no consensus in the literature
on which is more informative. Most studies use the PC, possibly because
it appears first in the algorithm's results. Halvorsen (2013) concluded
that the PC is most informative, whereas (Searcy and Shaffer 2016)
suggested that the PI is better for describing the biological properties
of species. Our study provides insights into the environmental factors
important for the distribution of sea snakes, but further research is
necessary to establish causal relationships between these factors and
the distribution of sea snakes, as well as to gain a better
understanding of their ecology and behavior.

In summary, our results on the relative importance of environmental
variables in shaping the distribution and potentially the abundance of
sea snakes, are consistent with previous research on marine organisms
and provide a valuable framework for developing effective conservation
management plans for sea snakes.

\textbf{4.1.1 Subfamily level}

Previous research has identified important variables for the
Laticaudinae subfamily, which can provide context for our results.
Gherghel et al. (2018) found that the most relevant variables for
predicting the presence of the Laticaudinae subfamily were distal
variables such as distance to the coast and bathymetry. We excluded
these distal variables from our analysis because they may not directly
reflect the ecological and physiological requirements of sea snake
species. The inclusion of distal variables is not recommended when an
interpretation of a species' biology is desired (Merow et al. 2014).
Instead, we focused on environmental variables directly related to
physical and resource factors. On the other hand, as there are no
definitive conclusions in the literature regarding the Hydrophiinae
subfamily, which comprises the vast majority of species within the
lineage, a specific discussion of this subfamily would entail repeating
information already covered above. Therefore, we have not included a
dedicated discussion of this subfamily

Overall, it appears that important variables for the Hydrophiinae and
Laticaudinae subfamilies may vary depending on the particular group of
snakes being studied, as well as the specific methods used for analysis.
Our results add to the growing body of literature on the ecology and
distribution of sea snakes, and provide further insights into the
factors that may influence their distribution patterns.

\textbf{4.1.2 Genus level}

Environmental variables play a significant role in shaping the abundance
and distribution of prey species, influencing the availability of
suitable habitat for sea snakes. Studies on \emph{Aipysurus} indicate
that environmental factors play crucial roles in shaping their
distribution and abundance. (Lynch et al. 2023) reported that
environmental factors are important in shaping the behavior and
distribution of \emph{A. laevis} , while (Borczyk et al. 2021) found
sexual dimorphism in skull size and shape for \emph{A. eydouxii} , which
may affect feeding and ecological niche differentiation between males
and females.

Previous studies have also highlighted the importance of environmental
variables in the distribution and population dynamics
of\emph{Emydocephalus} . (Lukoschek et al. 2007) found that \emph{E.
annulatus} exhibited metapopulation dynamics influenced by prey
abundance, habitat quality, and complex reef habitats, while (Goiran and
Shine 2013) reported a decline in \emph{E. annulatus} abundance on a
protected coral reef system in the New Caledonian Lagoon, attributed to
changes in the physical environment resulting from anthropogenic
modification to the system. Our study found that Sal, Sil, and Tem are
significant variables for both PC and PI at different resolutions,
consistent with (Brischoux et al. 2012) finding that oceanic salinity
influences body condition. However, we found a discrepancy in our
results regarding the importance of Tem compared to a previous study
on\emph{E. annulatus} , which reported no significant effects of annual
thermal variation on population densities. Further research is necessary
to clarify these discrepancies and enhance our understanding of the
importance of environmental variables, such as salinity and other
nutrient variables, in the distribution and population dynamics
of\emph{Emydocephalus} .

For \emph{Ephalophis} , our findings align with the idea of a rapid
radiation of sea snakes, driven by adaptations to environmental factors.
The importance of adaptation to the marine environment is highlighted by
the variable importance of environmental factors observed
in\emph{Ephalophis} . This is supported by previous studies on sea
snakes in Australian waters (Sanders et al. 2013), molecular phylogeny
of sea snakes (Lukoschek and Keogh 2006), association between habitat
type and skin thickness in snakes (Shine et al. 2019), presence of
paddle-shaped tails in viviparous sea snakes (Sanders et al. 2012), and
the feeding strategies and foraging ecology of the Hydrophiinae sea
snakes, as analyzed in previous studies (Hampton 2019).

The \emph{Hydrophis} genus is found to have almost identical results to
the whole lineage analyzed, indicating a similar influence of
environmental variables on their distribution and abundance. Meanwhile,
discussing the \emph{Laticauda} genus alone may be redundant as it
represents the whole Laticaudinae subfamily. Further diverse ecological
research is needed to support the findings and test the relative
importance of variables for the \emph{Microcephalophis} genus. Overall,
the importance of environmental factors in shaping the ecology and
evolution of sea snakes is reinforced by the results obtained from our
analysis. We also emphasize the importance of further diverse ecological
research to support the findings and test the relative importance of
variables for this genus.

\textbf{4.2 TN\textsubscript{R} limits}

Sea snakes demonstrate a greater consistency in their upper thermal
limits compared to their lower thermal limits, according to our research
findings. This aligns with previous studies on some terrestrial (Araujo
et al. 2013; Sunday et al. 2011, 2012) and marine ectothermic species
(Grady et al. 2019; Nishizaki et al. 2015; Stuart-Smith et al. 2017).
The results imply that cold tolerance evolves more quickly than heat
tolerance in both endothermic and ectothermic organisms (Bennett et al.
2018). Sea snakes may have developed greater plasticity in their lower
thermal limits through an evolutionary process (Gunderson and Stillman
2015; Verberk et al. 2018), resulting in less variation in their upper
thermal limits. However, if sea temperatures continue to rise due to
climate change, the reduced variability in upper thermal limits could
have significant consequences for sea snakes (Araujo et al. 2013;
Udyawer et al. 2020).

Marine ectotherms have a low physiological acclimatization rate, making
them particularly susceptible to the effects of climate change (Deutsch
et al. 2008; Sunday et al. 2011). Previous studies have indicated that
sea snakes exist near their upper thermal limits (Heatwole et al. 2012),
and that tropical species generally do not surpass the maximum
temperatures of the organisms that have been collected (Deutsch et al.
2008; Huey et al. 2009). Although our study observed this pattern at the
lineage level, we obtained mixed results at the subfamily and genus
level, potentially owing to our limited sample size. Other factors,
including biotic interactions, geographic barriers, and dispersal
processes, may restrict species distribution and impact the capacity to
measure complete thermal niches (Bennett et al. 2018; Sunday et al.
2011, 2012).

The consequences of the observed asymmetric realized thermal limits for
sea snakes' physiological performance and their ecological interactions
with other species in their habitat are significant. Sea snakes are
critical predators in their marine ecosystems, and alterations in their
distribution or abundance could have cascading effects on the food webs
they inhabit (Lukoschek et al. 2013). Thus, the observed thermal limits
of sea snakes have implications not only for their own survival but also
for the health and functioning of marine ecosystems.

To address the conservation implications of our findings, it is crucial
to consider the potential impact of climate change on sea snakes and the
measures needed to protect them. Management and conservation efforts
must recognize the susceptibility of many marine species, including sea
snakes, to climate change (Somaweera et al. 2021). Our research
highlights the importance of conservation efforts that consider the
potential vulnerability of sea snakes to changes in sea temperature,
including identifying areas that provide thermal refugia for sea snakes,
monitoring populations for changes in distribution or abundance, and
implementing adaptive management strategies that can respond to changing
conditions (Udyawer et al. 2020).

While our study offers valuable insights into the realized thermal
limits of sea snakes, there are limitations and gaps in knowledge that
future research could address. For example, our sample size for some
subfamilies and genera was relatively small, which may have affected our
results. Increasing the sample size may enhance the robustness of our
findings. Moreover, future research could explore the potential for sea
snakes to acclimate or adapt to warming temperatures and investigate the
mechanisms underlying the observed pattern of asymmetric realized
thermal limits. Integrating physiological and genetic data could provide
a more comprehensive understanding of the thermal ecology of sea snakes
and their ability to adapt in the face of climate change.

\textbf{5 CONCLUSIONS}

Our study reveals the significant role of temperature, salinity, and
some nutrients (such as nitrates, phosphates and silicates) in shaping
the distribution of sea snakes. These variables are known to have
implications for the morphology, diving behavior, coral reefs health,
and feeding habits of sea snakes Moreover, our findings offer valuable
insights for conservation management strategies.

Furthermore, our study demonstrates the asymmetric realized thermal
limits of sea snakes, with less variation in their upper thermal limits
compared to their lower thermal limits. This pattern is consistent with
previous research on other endothermic and ectothermic species and has
critical implications for sea snakes under scenarios of sea warming.
Suggesting that sea snakes have accumulated greater plasticity in their
lower thermal limits through an evolutionary process.

\textbf{6 REFERENCES}

Araujo MB, Ferri-Yanez F, Bozinovic F, Marquet PA, Valladares F, Chown
SL (2013) Heat freezes niche evolution. Ecology Letters 16: 1206-1219
doi 10.1111/ele.12155

Ashcroft MB, French KO, Chisholm LA (2011) An evaluation of
environmental factors affecting species distributions. Ecological
Modelling 222: 524-531
doi\emph{https://doi.org/10.1016/j.ecolmodel.2010.10.003}

Assis J, Tyberghein L, Bosch S, Verbruggen H, Serrao EA, De Clerck O
(2018) Bio-ORACLE v2.0: Extending marine data layers for bioclimatic
modelling. Global Ecology and Biogeography 27: 277-284 doi
10.1111/geb.12693

Bennett JM, Calosi P, Clusella-Trullas S, Martinez B, Sunday J, Algar
AC, Araujo MB, Hawkins BA, Keith S, Kuhn I, Rahbek C, Rodriguez L,
Singer A, Villalobos F, Angel Olalla-Tarraga M, Morales-Castilla I
(2018) GlobTherm, a global database on thermal tolerances for aquatic
and terrestrial organisms. Scientific Data 5: 180022 doi
10.1038/sdata.2018.22

Bentlage B, Peterson AT, Barve N, Cartwright P (2013) Plumbing the
depths: extending ecological niche modelling and species distribution
modelling in three dimensions. Global Ecology and Biogeography 22:
952-961 doi\emph{https://doi.org/10.1111/geb.12049}

Bessesen BL, Garrido-Cayul C, Gonzalez-Suarez M (2023) Habitat
suitability and area of occupancy defined for rare New World sea snake.
Conservation Science and Practice 5: e12865
doi\emph{https://doi.org/10.1111/csp2.12865}

Bicego KC, Barros RCH, Branco LGS (2007) Physiology of temperature
regulation: Comparative aspects. Comparative Biochemistry and Physiology
Part A: Molecular \& Integrative Physiology 147: 616-639
doi\emph{https://doi.org/10.1016/j.cbpa.2006.06.032}

Borczyk B, Pa\selectlanguage{polish}ś\selectlanguage{english}ko \selectlanguage{polish}Ł, \selectlanguage{english}Kusznierz J, Bury S (2021) Sexual dimorphism and
skull size and shape in the highly specialized snake species, Aipysurus
eydouxii (Elapidae: Hydrophiinae). PeerJ 9: e11311 doi
10.7717/peerj.11311

Bosch S, Tyberghein L, Deneudt K, Hernandez F, De Clerck O (2018) In
search of relevant predictors for marine species distribution modelling
using the MarineSPEED benchmark dataset. Diversity and Distributions 24:
144-157 doi\emph{https://doi.org/10.1111/ddi.12668}

Bradie J, Leung B (2017) A quantitative synthesis of the importance of
variables used in MaxEnt species distribution models. Journal of
Biogeography 44: 1344-1361 doi 10.1111/jbi.12894

Brischoux F, Tingley R, Shine R, Lillywhite HB (2012) Salinity
influences the distribution of marine snakes: implications for
evolutionary transitions to marine life. Ecography 35: 994-1003 doi
10.1111/j.1600-0587.2012.07717.x

Brown JL (2014) SDMtoolbox: a python-based GIS toolkit for landscape
genetic, biogeographic and species distribution model analyses. Methods
in Ecology and Evolution 5: 694-700 doi 10.1111/2041-210x.12200

Clarke A (2003) Costs and consequences of evolutionary temperature
adaptation. Trends in Ecology \& Evolution 18: 573-581
doi\emph{https://doi.org/10.1016/j.tree.2003.08.007}

Cobos ME, Peterson AT, Barve N, Osorio-Olvera L (2019) kuenm: an R
package for detailed development of ecological niche models using
Maxent. PeerJ 7:e6281: 15 doi 10.7717/peerj.6281

Collin R, Rebolledo AP, Smith E, Chan KYK (2021) Thermal tolerance of
early development predicts the realized thermal niche in marine
ectotherms. Functional Ecology 35: 1679-1692
doi\emph{https://doi.org/10.1111/1365-2435.13850}

Colwell RK, Rangel TF (2010) A stochastic, evolutionary model for range
shifts and richness on tropical elevational gradients under Quaternary
glacial cycles. Philosophical Transactions of the Royal Society B:
Biological Sciences 365: 3695-3707 doi 10.1098/rstb.2010.0293

Connor T, Vi\selectlanguage{ngerman}ña A, Winkler JA, Hull V, Tang Y, Shortridge A, Yang H, Zhao
Z, Wang F, Zhang J, Zhang Z, Zhou C, Bai W, Liu J (2019) Interactive
spatial scale effects on species distribution modeling: The case of the
giant panda. Scientific Reports 9: 14563 doi 10.1038/s41598-019-50953-z

Contreras-Díaz RG, Falconi M, Osorio-Olvera L, Cobos ME, Soberón J,
Townsend Peterson A, Lira-Noriega A, Álvarez-Loayza P, Luis Gonçalves A,
Hurtado-Astaiza J, Gonzáles RdPR, Zubileta IS, Spironello WR,
Vásquez-Martínez R (2022) On the relationship between environmental
suitability and habitat use for three neotropical mammals. Journal of
Mammalogy 103: 425-439 doi 10.1093/jmammal/gyab152

Dabruzzi T, F., Sutton M, A., Bennett W, A. (2012) Metabolic Thermal
Sensitivity Optimizes Sea Krait Amphibious Physiology. Herpetologica 68:
218-225 doi 10.1655/HERPETOLOGICA-D-11-00077.1

Descombes P, Wisz MS, Leprieur F, Parravicini V, Heine C, Olsen SM,
Swingedouw D, Kulbicki M, Mouillot D, Pellissier L (2015) Forecasted
coral reef decline in marine biodiversity hotspots under climate change.
Global change biology 21: 2479-2487 doi 10.1111/gcb.12868

Deutsch CA, Tewksbury JJ, Huey RB, Sheldon KS, Ghalambor CK, Haak DC,
Martin PR (2008) Impacts of climate warming on terrestrial ectotherms
across latitude. Proceedings of the National Academy of Sciences 105:
6668 doi 10.1073/pnas.0709472105

Dong Y, Somero GN (2009) Temperature adaptation of cytosolic malate
dehydrogenases of limpets(genus Lottia): differences in stability and
function due to minor changes in sequence correlate with biogeographic
and vertical distributions. Journal of Experimental Biology 212: 169-177
doi 10.1242/jeb.024505

Eme D, Malard F, Colson-Proch C, Jean P, Calvignac S, Konecny-Dupré L,
Hervant F, Douady CJ (2014) Integrating phylogeography, physiology and
habitat modelling to explore species range determinants. Journal of
Biogeography 41: 687-699 doi\emph{https://doi.org/10.1111/jbi.12237}

Feng X, Park DS, Liang Y, Pandey R, Pape\selectlanguage{polish}ş \selectlanguage{english}M (2019) Collinearity in
ecological niche modeling: Confusions and challenges. Ecology and
Evolution 9: 10365-10376 doi\emph{https://doi.org/10.1002/ece3.5555}

Fligner MA, Killeen TJ (1976) Distribution-Free Two-Sample Tests for
Scale. Journal of the American Statistical Association 71: 210-213 doi
10.2307/2285771

Fraser D (1999) Animal ethics and animal welfare science: bridging the
two cultures. Applied Animal Behaviour Science 65: 171-189
doi\emph{https://doi.org/10.1016/S0168-1591(99)00090-8}

Gaston KJ, Fuller RA (2009) The sizes of species' geographic ranges.
Journal of Applied Ecology 46: 1-9 doi 10.1111/j.1365-2664.2008.01596.x

Gause GF (1931) The Influence of Ecological Factors on the Size of
Population. The American naturalist 65: 70-76

Ga\selectlanguage{ngerman}üzère P, Jiguet F, Devictor V (2015) Rapid adjustment of bird
community compositions to local climatic variations and its functional
consequences. Global change biology 21: 3367-3378 doi 10.1111/gcb.12917

Gherghel I, Brischoux F, Pape\selectlanguage{polish}ş \selectlanguage{english}M (2018) Using biotic interactions in
broad-scale estimates of species' distributions. Journal of Biogeography
45: 2216-2225 doi 10.1111/jbi.13361

Gillespie TW, Foody GM, Rocchini D, Giorgi AP, Saatchi S (2008)
Measuring and modelling biodiversity from space. Progress in Physical
Geography: Earth and Environment 32: 203-221 doi
10.1177/0309133308093606

Goiran C, Shine R (2013) Decline in sea snake abundance on a protected
coral reef system in the New Caledonian Lagoon. Coral Reefs 32: 281-284
doi 10.1007/s00338-012-0977-x

Grady JM, Maitner BS, Winter AS, Kaschner K, Tittensor DP, Record S,
Smith FA, Wilson AM, Dell AI, Zarnetske PL, Wearing HJ, Alfaro B, Brown
JH (2019) Metabolic asymmetry and the global diversity of marine
predators. Science 363: eaat4220 doi 10.1126/science.aat4220

Guisan A, Thuiller W (2005) Predicting species distribution: offering
more than simple habitat models. Ecology Letters 8: 993-1009
doi\emph{https://doi.org/10.1111/j.1461-0248.2005.00792.x}

Gunderson AR, Stillman JH (2015) Plasticity in thermal tolerance has
limited potential to buffer ectotherms from global warming. Proceedings
of the Royal Society B: Biological Sciences 282: 20150401 doi
10.1098/rspb.2015.0401

Gvo\selectlanguage{polish}ž\selectlanguage{english}d\selectlanguage{ngerman}ík L (2018) Just what is the thermal niche? Oikos 127: 1701-1710
doi 10.1111/oik.05563

Halvorsen R (2013) A strict maximum likelihood explanation of MaxEnt,
and some implications for distribution modelling. Sommerfeltia 36: 1-132
doi doi:10.2478/v10208-011-0016-2

Hampton PM (2019) Foraging ecology influences the number of vertebrae in
hydrophiine sea snakes. Biological Journal of the Linnean Society 128:
645-650 doi 10.1093/biolinnean/blz115

Heatwole H, Grech A, Monahan JF, King S, Marsh H (2012) Thermal biology
of sea snakes and sea kraits. Integrative and comparative biology 52:
257-273 doi 10.1093/icb/ics080

Helmuth B, Broitman BR, Blanchette CA, Gilman S, Halpin P, Harley CDG,
O'Donnell MJ, Hofmann GE, Menge B, Strickland D (2006) Mosaic patterns
of thermal stress in the rocky intertidal zone: implications for climate
change. Ecological Monographs 76: 461-479
doi\emph{https://doi.org/10.1890/0012-9615(2006)076{[}0461:MPOTSI{]}2.0.CO;2}

Hijmans RJ, Phillips S, Leathwick J, Elith J, Hijmans MRJJC (2017)
Package `dismo'9: 1-68

Huey RB, Deutsch CA, Tewksbury JJ, Vitt LJ, Hertz PE, Álvarez Pérez HJ,
Garland T (2009) Why tropical forest lizards are vulnerable to climate
warming. Proceedings of the Royal Society B: Biological Sciences 276:
1939-1948 doi 10.1098/rspb.2008.1957

Huff DD, Hubler SL, Borisenko AN (2005) Using Field Data to Estimate the
Realized Thermal Niche of Aquatic Vertebrates. North American Journal of
Fisheries Management 25: 346-360 doi 10.1577/M03-231.1

Hutchings JA, Myers RA (1994) Timing of cod reproduction: interannual
variability and the influence of temperature. Marine Ecology Progress
Series 108: 21-32

Kordas RL, Harley CDG, O'Connor MI (2011) Community ecology in a warming
world: The influence of temperature on interspecific interactions in
marine systems. Journal of Experimental Marine Biology and Ecology 400:
218-226 doi\emph{https://doi.org/10.1016/j.jembe.2011.02.029}

Kuo ESL, Sanford E (2009) Geographic variation in the upper thermal
limits of an intertidal snail: implications for climate envelope models.
Marine Ecology Progress Series 388: 137-146

Lee-Yaw JA, Kharouba HM, Bontrager M, Mahony C, Cserg\selectlanguage{polish}ő \selectlanguage{english}AM, Noreen AME,
Li Q, Schuster R, Angert AL (2016) A synthesis of transplant experiments
and ecological niche models suggests that range limits are often niche
limits. Ecology Letters 19: 710-722 doi 10.1111/ele.12604

Lukoschek V, Beger M, Ceccarelli D, Richards Z, Pratchett M (2013)
Enigmatic declines of Australia's sea snakes from a biodiversity
hotspot. Biological Conservation 166: 191-202
doi\emph{https://doi.org/10.1016/j.biocon.2013.07.004}

Lukoschek V, Heatwole H, Grech A, Burns G, Marsh H (2007) Distribution
of two species of sea snakes, Aipysurus laevis and Emydocephalus
annulatus, in the southern Great Barrier Reef: metapopulation dynamics,
marine protected areas and conservation. Coral Reefs 26: 291-307 doi
10.1007/s00338-006-0192-8

Lukoschek V, Keogh JS (2006) Molecular phylogeny of sea snakes reveals a
rapidly diverged adaptive radiation. Biological Journal of the Linnean
Society 89: 523-539 doi 10.1111/j.1095-8312.2006.00691.x

Lynch TP, Alford RA, Shine R (2023) Sex-based divergence in tidal, lunar
and seasonal cycles of activity in the olive sea snake, Aipysurus laevis
(Elapidae, Hydrophiinae). Coral Reefs doi 10.1007/s00338-022-02341-9

Magnuson JJ, Crowder LB, Patricia AM (1979) Temperature as an Ecological
Resource. American Zoologist 19: 331-343

Mateo RG, Vanderpoorten A, Mu\selectlanguage{ngerman}ñoz J, Laenen B, Désamoré A (2013) Modeling
Species Distributions from Heterogeneous Data for the Biogeographic
Regionalization of the European Bryophyte Flora. PLOS ONE 8: e55648 doi
10.1371/journal.pone.0055648

Merow C, Smith MJ, Edwards Jr TC, Guisan A, McMahon SM, Normand S,
Thuiller W, Wüest RO, Zimmermann NE, Elith J (2014) What do we gain from
simplicity versus complexity in species distribution models? Ecography
37: 1267-1281 doi\emph{https://doi.org/10.1111/ecog.00845}

Muscarella R, Galante PJ, Soley-Guardia M, Boria RA, Kass JM, Uriarte M,
Anderson RP (2014) ENMeval: An R package for conducting spatially
independent evaluations and estimating optimal model complexity for
Maxent ecological niche models. Methods in Ecology and Evolution 5:
1198-1205 doi 10.1111/2041-210x.12261

Nishizaki MT, Barron S, Carew E (2015) Thermal stress increases
fluctuating asymmetry in marine mussels: environmental variation and
developmental instability. Ecosphere 6: art85
doi\emph{https://doi.org/10.1890/ES14-00399.1}

Osorio-Olvera L, Lira-Noriega A, Soberón J, Peterson AT, Falconi M,
Contreras-Díaz RG, Martínez-Meyer E, Barve V, Barve N (2020) ntbox: An r
package with graphical user interface for modelling and evaluating
multidimensional ecological niches. Methods in Ecology and Evolution 11:
1199-1206 doi\emph{https://doi.org/10.1111/2041-210X.13452}

Pearson RG, Dawson TP (2003) Predicting the impacts of climate change on
the distribution of species: are bioclimate envelope models useful?
Global Ecology and Biogeography 12: 361-371
doi\emph{https://doi.org/10.1046/j.1466-822X.2003.00042.x}

Peterson A, Soberón J (2012) Species Distribution Modeling and
Ecological Niche Modeling: Getting the Concepts Right. J Natureza \&
Conservação 10: 102-107

Peterson AT, Pape\selectlanguage{polish}ş \selectlanguage{english}M, Sober\selectlanguage{ngerman}ón J (2008) Rethinking receiver operating
characteristic analysis applications in ecological niche modeling.
Ecological Modelling 213: 63-72
doi\emph{https://doi.org/10.1016/j.ecolmodel.2007.11.008}

Phillips SJ (2011) A brief tutorial on Maxent. AT\&T Research,
Princeton, NJ

Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of
species geographic distributions. Ecological Modelling 190: 231-259
doi\emph{https://doi.org/10.1016/j.ecolmodel.2005.03.026}

Pörtner H (2001) Climate change and temperature-dependent biogeography:
oxygen limitation of thermal tolerance in animals. Naturwissenschaften
88: 137-146 doi 10.1007/s001140100216

Pörtner HO, Farrell AP (2008) Physiology and Climate Change. Science
322: 690-692 doi 10.1126/science.1163156 \%J Science

Quintero I, Wiens JJ (2013) Rates of projected climate change
dramatically exceed past rates of climatic niche evolution among
vertebrate species. Ecology Letters 16: 1095-1103
doi\emph{https://doi.org/10.1111/ele.12144}

R Core Team (2013) R: A language and environment for statistical
computing. Foundation for Statistical Computing, Vienna, Austria

Rabosky DL, Grundler M, Anderson C, Title P, Shi JJ, Brown JW, Huang H,
Larson JG (2014) BAMMtools: an R package for the analysis of
evolutionary dynamics on phylogenetic trees. Methods in Ecology and
Evolution 5: 701-707 doi\emph{https://doi.org/10.1111/2041-210X.12199}

Radosavljevic A, Anderson RP (2014) Making better Maxent models of
species distributions: complexity, overfitting and evaluation. Journal
of Biogeography 41: 629-643 doi 10.1111/jbi.12227

Randin CF, Jaccard H, Vittoz P, Yoccoz NG, Guisan A (2009) Land use
improves spatial predictions of mountain plant abundance but not
presence-absence. Journal of Vegetation Science 20: 996-1008 doi
10.1111/j.1654-1103.2009.01098.x

Robertson MP, Visser V, Hui C (2016) Biogeo: an R package for assessing
and improving data quality of occurrence record datasets. Ecography 39:
394-401 doi\emph{https://doi.org/10.1111/ecog.02118}

Sánchez-Fernández D, Aragón P, Bilton DT, Lobo JM (2012) Assessing the
Congruence of Thermal Niche Estimations Derived from Distribution and
Physiological Data. A Test Using Diving Beetles. PLOS ONE 7: e48163 doi
10.1371/journal.pone.0048163

Sánchez-Fernández D, Rizzo V, Cieslak A, Faille A, Fresneda J, Ribera I
(2016) Thermal niche estimators and the capability of poor dispersal
species to cope with climate change. Scientific Reports 6: 23381 doi
10.1038/srep23381

Sanders KL, Lee MSY, Mumpuni, Bertozzi T, Rasmussen AR (2013) Multilocus
phylogeny and recent rapid radiation of the viviparous sea snakes
(Elapidae: Hydrophiinae). Molecular phylogenetics and evolution 66:
575-591 doi\emph{https://doi.org/10.1016/j.ympev.2012.09.021}

Sanders KL, Rasmussen AR, Elmberg J (2012) Independent Innovation in the
Evolution of Paddle-Shaped Tails in Viviparous Sea Snakes (Elapidae:
Hydrophiinae). Integrative and comparative biology 52: 311-320 doi
10.1093/icb/ics066

Sbrocco EJ, Barber PH (2013) MARSPEC: ocean climate layers for marine
spatial ecology. Ecology 94: 979-979 doi 10.1890/12-1358.1

Searcy CA, Shaffer HB (2016) Do Ecological Niche Models Accurately
Identify Climatic Determinants of Species Ranges? The American
naturalist 187: 423-435 doi 10.1086/685387

Sexton JP, McIntyre PJ, Angert AL, Rice KJ (2009) Evolution and Ecology
of Species Range Limits. Annual Review of Ecology, Evolution, and
Systematics 40: 415-436 doi 10.1146/annurev.ecolsys.110308.120317

Shcheglovitova M, Anderson RP (2013) Estimating optimal complexity for
ecological niche models: A jackknife approach for species with small
sample sizes. Ecological Modelling 269: 9-17
doi\emph{https://doi.org/10.1016/j.ecolmodel.2013.08.011}

Sherratt E, Rasmussen AR, Sanders KL (2018) Trophic specialization
drives morphological evolution in sea snakes. Royal Society Open Science
5: 172141 doi 10.1098/rsos.172141

Shine R, Goiran C, Shilton C, Meiri S, Brown GP (2019) The life aquatic:
an association between habitat type and skin thickness in snakes.
Biological Journal of the Linnean Society 128: 975-986 doi
10.1093/biolinnean/blz136

Smith AB, Santos MJ (2020) Testing the ability of species distribution
models to infer variable importance. Ecography 43: 1-13 doi
10.1111/ecog.05317

Soberon J, Peterson AT (2005) Interpretation of Models of Fundamental
Ecological Niches and Species' Distributional Areas. Biodiversity
Informatics 2: 10 doi 10.17161/bi.v2i0.4

Somaweera R, Udyawer V, Guinea ML, Ceccarelli DM, Clarke RH, Glover M,
Hourston M, Keesing J, Rasmussen AR, Sanders K, Shine R, Thomson DP,
Webber BL (2021) Pinpointing Drivers of Extirpation in Sea Snakes: A
Synthesis of Evidence From Ashmore Reef. Frontiers in Marine Science 8

Spalding MD, Fox HE, Allen GR, Davidson N, Ferdaña ZA, Finlayson M,
Halpern BS, Jorge MA, Lombana A, Lourie SA, Martin KD, McManus E, Molnar
J, Recchia CA, Robertson J (2007) Marine Ecoregions of the World: A
Bioregionalization of Coastal and Shelf Areas. BioScience 57: 573-583
doi 10.1641/B570707 \%J BioScience

Stillman JH, Somero GN (2000) A Comparative Analysis of the Upper
Thermal Tolerance Limits of Eastern Pacific Porcelain Crabs, Genus
Petrolisthes: Influences of Latitude, Vertical Zonation, Acclimation,
and Phylogeny. Physiological and Biochemical Zoology 73: 200-208 doi
10.1086/316738

Stuart-Smith RD, Edgar GJ, Bates AE (2017) Thermal limits to the
geographic distributions of shallow-water marine species. Nature Ecology
\& Evolution 1: 1846-1852 doi 10.1038/s41559-017-0353-x

Sunday JM, Bates AE, Dulvy NK (2011) Global analysis of thermal
tolerance and latitude in ectotherms. Proceedings of Royal Society B
Biological Sciences 278: 1823-1830 doi doi:10.1098/rspb.2010.1295

Sunday JM, Bates AE, Dulvy NK (2012) Thermal tolerance and the global
redistribution of~animals. Nature Climate Change 2: 686-690 doi
10.1038/nclimate1539

Udyawer V, Barnes P, Bonnet X, Brischoux F, Crowe-Riddell JM, D'Anastasi
B, Fry BG, Gillett A, Goiran C, Guinea ML, Heatwole H, Heupel MR,
Hourston M, Kangas M, Kendrick A, Koefoed I, Lillywhite HB, Lobo AS,
Lukoschek V, McAuley R, Nitschke C, Rasmussen AR, Sanders KL, Sheehy C,
Shine R, Somaweera R, Sweet SS, Voris HK (2018) Future Directions in the
Research and Management of Marine Snakes. Frontiers in Marine Science 5:
399 doi 10.3389/fmars.2018.00399

Udyawer V, Cappo M, Simpfendorfer CA, Heupel MR, Lukoschek V (2014)
Distribution of sea snakes in the Great Barrier Reef Marine Park:
observations from 10~yrs of baited remote underwater video station
(BRUVS) sampling. Coral Reefs 33: 777-791 doi 10.1007/s00338-014-1152-3

Udyawer V, Simpfendorfer CA, Heupel MR, Clark TD (2016) Coming up for
air: thermal dependence of dive behaviours and metabolism in sea snakes.
Journal of Experimental Biology 219: 3447-3454 doi 10.1242/jeb.146571

Udyawer V, Somaweera R, Nitschke C, d'Anastasi B, Sanders K, Webber BL,
Hourston M, Heupel MR (2020) Prioritising search effort to locate
previously unknown populations of endangered marine reptiles. Global
Ecology and Conservation 22: e01013
doi\emph{https://doi.org/10.1016/j.gecco.2020.e01013}

Verberk WCEP, Calosi P, Spicer JI, Kehl S, Bilton DT (2018) Does
plasticity in thermal tolerance trade off with inherent tolerance? The
influence of setal tracheal gills on thermal tolerance and its
plasticity in a group of European diving beetles. Journal of Insect
Physiology 106: 163-171
doi\emph{https://doi.org/10.1016/j.jinsphys.2017.12.005}

Veron JEN (2000) Corals of the World

Veron JEN (2008) A reef in time: the Great Barrier Reef from beginning
to end. Harvard University Press

Wallach V, Williams KL, Boundy J (2014) Snakes of the world: a catalogue
of living and extinct species. CRC press

Walther G-R, Post E, Convey P, Menzel A, Parmesan C, Beebee TJC,
Fromentin J-M, Hoegh-Guldberg O, Bairlein F (2002) Ecological responses
to recent climate change. Nature 416: 389-395 doi 10.1038/416389a

Walther GR, Roques A, Hulme PE, Sykes MT, Pysek P, Kühn I, Zobel M,
Bacher S, Botta-Dukát Z, Bugmann H, Czúcz B, Dauber J, Hickler T,
Jarosík V, Kenis M, Klotz S, Minchin D, Moora M, Nentwig W, Ott J, Panov
VE, Reineking B, Robinet C, Semenchenko V, Solarz W, Thuiller W, Vilà M,
Vohland K, Settele J (2009) Alien species in a warmer world: risks and
opportunities. Trends Ecol Evol 24: 686-693 doi
10.1016/j.tree.2009.06.008

Warren DL, Seifert SN (2011) Ecological niche modeling in Maxent: the
importance of model complexity and the performance of model selection
criteria. Ecological Applications 21: 335-342 doi 10.1890/10-1171.1

Wiens JA, Stralberg D, Jongsomjit D, Howell CA, Snyder MA (2009) Niches,
models, and climate change: Assessing the assumptions and uncertainties.
Proceedings of the National Academy of Sciences 106: 19729 doi
10.1073/pnas.0901639106

Wiens JJ (2011) The niche, biogeography and species interactions.
Philosophical transactions of the Royal Society of London Series B,
Biological sciences 366: 2336-2350 doi 10.1098/rstb.2011.0059

\textbf{7 TABLE LEGENDS}

\textbf{Table 1.} Variable groups with description of each individual
variable by repository and depth. VCode = variable code. Bma =
Bio-Oracle benthic maximum, Bme = Bio-Oracle benthic mean, Bmi =
Bio-Oracle benthic minimum, Bsurf = Bio-Oracle surface, Msurf = MARSPEC
surface.

\textbf{Table 2.} Summary of inputs and evaluations from ecological
niche modeling. Scode = species code, Tocc = total occurrences, Focc =
final post-processed occurrences, M5 = number of final models at 5
arc-minutes, M10 = number of final models at 10 arc-minutes, AUCR5 =
median of AUC ratios at 5 arc-minutes, AUCR10 = median of AUC ratios at
10 arc-minutes, OR5 = median of omission rate at 5 arc-minutes, OR10 =
median of omission rate at 5 arc-minutes, * = omission rate
\textgreater{} 10\%.

\textbf{8 FIGURE LEGENDS}

\textbf{Figure 1.} All variable contribution estimated by MaxEnt
percentage of contribution (PC) and permutation of importance (PI) for
each taxonomic level at both spatial resolutions. A) = 5 arc-minutes, B)
= 10 arc-minutes. Warmer colors denoted higher contribution while
lighter colors less contribution.

\textbf{Figure 2.} Chord diagram for the relative important variables
under the percentage of contribution (PC) and permutation of importance
(PI) estimated by MaxEnt and their representation in different taxonomic
and spatial resolution levels. More connected variables represent that
are relatively important variables in more taxonomic and spatial
resolution levels comparing to variables less connected. A) = PC at 5
arc-minutes, B) = PI at 5 arc-minutes, C) = PC at 10 arc-minutes, D) =
PI at 10 arc-minutes, Nit = nitrates, Pho = phosphates, Sal = salinity,
Sil = silicates, Tem = temperature.

\textbf{Figure 3.} Lower and upper thermal limits of sea snakes at three
taxonomic levels, three temperatures measures, five variable groups and
two spatial resolutions. A) = Sea snakes' lineage at 5 arc-minutes, B) =
Sea snakes' lineage at 10 arc-minutes, C) = Hydrophiinae subfamily at 5
arc-minutes, D) = Hydrophiinae subfamily at 10 arc-minutes, E) =
Laticaudinae subfamily at 5 arc-minutes, F) = Laticaudinae subfamily at
10 arc-minutes, G) = \emph{Aipysurus} genus at 5 arc-minutes, H)
=\emph{Aipysurus} genus at 10 arc-minutes, I) =
\emph{Emydocephalus}genus at 5 arc-minutes, H) = \emph{Emydocephalus}
genus at 10 arc-minutes, Max = maximum temperature, Mean = mean
temperature, Min = minimum temperature, Bma = maximum benthic
Bio-Oracle, Bme = mean benthic Bio-Oracle, Bmi = minimum benthic
Bio-Oracle, Bsurf = superficial Bio-Oracle, Msurf = superficial MARSPEC.
All \emph{p} values are in Table SM4.

\textbf{9 SUPPORTING INFORMATION}

\textbf{Table SM1.} List of dataset keys (DOIs) and occurrence count for
all GBIF data used.

\textbf{Table SM2.} List of final configuration models by species.

\textbf{Table SM3.} Non-parametric Wilxicon Test comparison between
medians of variables at 5 and 10 arc-minutes of spatial resolution and
at relative importance metric estimated by MaxEnt for lineage, family
and genus level.

\textbf{Table SM4.} Summary of \emph{p} values for homogeneity of
variances by non-parametric Fligner-Killeen test.

\textbf{Figure SM1.} Density map of Sea snakes' lineage by genus. Bigger
points represent more density of presence records and smaller points
less presence records.

\textbf{Figure SM2-51} . Suitability and range of variability maps at 5
and 10 arc-minutes for each species.

\textbf{Figure SM52.} All variable contribution estimated by MaxEnt
percentage of contribution (PC) and permutation of importance (PI) for
each species at both spatial resolutions.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-1/Fig-1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-2/Fig-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-3/Fig-3}

\end{center}
\end{figure}

\textbf{Hosted file}

\verb`Table_1.docx`  available at \url{https://authorea.com/users/627538/articles/708286-global-analysis-of-the-influence-of-environmental-variables-to-explain-distributions-and-realized-thermal-niche-boundaries-of-sea-snakes}

\textbf{Hosted file}

\verb`Table_2.docx`  available at \url{https://authorea.com/users/627538/articles/708286-global-analysis-of-the-influence-of-environmental-variables-to-explain-distributions-and-realized-thermal-niche-boundaries-of-sea-snakes}

\selectlanguage{english}
\FloatBarrier
\end{document}

