\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[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Population and habitat responses to the global climate change in a
widespread species, the Asiatic Toad (Bufo gargarizans)}



\author[1]{Liuyang Yang}%
\author[1]{Huabin Zhang}%
\author[1]{Naijing Zhang}%
\author[1]{Ruiqing Cai}%
\author[2]{Peng Yan}%
\author[1]{En Li}%
\author[1]{Tao Pan}%
\author[1]{Xiaobing Wu}%
\affil[1]{Affiliation not available}%
\affil[2]{Anhui Normal University}%


\vspace{-1em}



  
  \date{June 17, 2024}


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





\selectlanguage{english}
\begin{abstract}
Abstract: Global climate change is gradually changing the distribution
and diversity pattern of species. In order to study the change of
distribution area and distribution barycenter of Bufo gargarizans under
climate change, we adopted the maximum-entropy modeling (Maxent) model
and barycenter migration analysis. And four General circulation models
(GCMs) under four representative concentration pathways (RCPs), 29
environment variables (including bioclimate, topography, habitat and
human impact) were used to predict the current and future (2050 and
2070) distribution of B. gargarizans. The results show that the Mean
temperature of Driest quarter (Bio9), Precipitation of driest month
(Bio14) and population (POP) have greater influence on the distribution
of B. gargarizans, it indicates that climate and anthropogenic factors
have greater influence on the distribution of B. gargarizans than other
environmental factors such as topography and landform. And in the next
30 to 50 years, from the perspective of the large spatial scale of the
whole China, the distribution area of B. gargarizans is in the eastern
part of China, the southeastern part is reduced, and the Midwest and
northwest part is expanded. On the whole, the suitable habitat of B.
gargarizans will be reduced in the future and will migrate to the high
latitude and altitude area. In this study, B. gargarizans was used as a
sentinel species to study the response of amphibians to climate change,
this paper aims to reveal the mechanism of amphibian response to climate
change from the current and future dispersal patterns of B. gargarizans.%
\end{abstract}%



\sloppy


\textbf{Population and habitat responses to the global climate change in
a widespread species, the Asiatic Toad (\emph{Bufo gargarizans})}

Liuyang Yang\selectlanguage{ngerman}\textsuperscript{a,b¶}, Huabin Zhang\textsuperscript{a,c¶},
Naijing Zhang\textsuperscript{a}, Ruiqing Cai\textsuperscript{a}, Peng
Yan\textsuperscript{a}, En Li \textsuperscript{a*}, Tao
Pan\textsuperscript{a *}, Xiaobing Wu \textsuperscript{a*}

\textsuperscript{a}Anhui Province Key Laboratory of the Conservation and
Exploitation of Biological Resources, College of Life Sciences, Anhui
Normal University, Wuhu 241000, Anhui, China;

\textsuperscript{¶} These authors contributed equally to this work.

\textsuperscript{*} Author for correspondence: Xiaobing Wu,
wuxb@ahnu.edu.cn; Tao Pan, pantao@ahnu.edu.cn; En Li,
lien1971@mail.ahnu.edu.cn.

Address of corresponding author (they have the same address): Zhe Shan
Campus, Anhui Normal University, Jinghu District, Wuhu city, Anhui
Province, China.

Running Head: The Asiatic Toad respond to the climate change

\textbf{Abstract:}Global climate change is gradually changing the
distribution and diversity pattern of species. In order to study the
change of distribution area and distribution barycenter of\emph{Bufo
gargarizans} under climate change, we adopted the maximum-entropy
modeling (Maxent) model and barycenter migration analysis. And four
General circulation models (GCMs) under four representative
concentration pathways (RCPs), 29 environment variables (including
bioclimate, topography, habitat and human impact) were used to predict
the current and future (2050 and 2070) distribution of\emph{B.
gargarizans.} The results show that the Mean temperature of Driest
quarter (Bio9), Precipitation of driest month (Bio14) and population
(POP) have greater influence on the distribution of \emph{B.
gargarizans} , it indicates that climate and anthropogenic factors have
greater influence on the distribution of \emph{B. gargarizans} than
other environmental factors such as topography and landform. And in the
next 30 to 50 years, from the perspective of the large spatial scale of
the whole China, the distribution area of \emph{B. gargarizans} is in
the eastern part of China, the southeastern part is reduced, and the
Midwest and northwest part is expanded. On the whole, the suitable
habitat of \emph{B. gargarizans} will be reduced in the future and will
migrate to the high latitude and altitude area. In this study, \emph{B.
gargarizans} was used as a sentinel species to study the response of
amphibians to climate change, this paper aims to reveal the mechanism of
amphibian response to climate change from the current and future
dispersal patterns of \emph{B. gargarizans} .

\textbf{Key words:} \emph{Bufo gargarizans} , Maxent, distribution,
barycenter migration, climate change, sentinel species



\textbf{1. Introduction}

The distribution and interaction of species are closely related to
climate (Nogues-Bravo \& Rahbek, 2011), and global climate change is a
major challenge facing mankind and the biological world in the 21st
century (Bertrand et al., 2011). According to the Intergovernment Panel
on Climate Change (IPCC), average global temperatures have risen by
about 1.5\selectlanguage{english} compared with pre-industrial levels. And the rising global
average temperature poses a serious threat to the sustainability of the
global ecosystem and alters global biodiversity patterns (Dawson,
Jackson, House, Prentice, \& Mace, 2011). If species can find niches to
satisfy their needs, their habitats will expand or shrink as the climate
changes. However, if climate-friendly habitats disappear or are affected
by landscape barriers, some less diffuse species will become extinct as
a result of climate warming.

Amphibians are at higher risk of extinction than other groups of
vertebrates(Stuart et al., 2004; Wake \& Vredenburg, 2008), with 41
percent at risk, and due to their complex life history and unique
physiological structure, their distribution and diffusion are greatly
affected by climate and landscape, and their ability to cope with
climate change is weak, which may be one of the reasons for the decline
of the global amphibians population (Lambers, 2015; Sinervo et al.,
2010).

In order to study the response of species to climate change, SDMs (SDMs)
is usually chosen to simulate the distribution area of species. SDMs can
calculate the relationship between the existence point of species and
local climate variables, perform function fitting through different
algorithms, and then project the function to a specific research area
and time, and finally obtain the distribution of the target species'
adaptive area (Elith, 2009; Guisan \& Thuiller, 2010). Although the SDMs
has some shortcomings, for example, it cannot fully take into account
the interactions between species (Pearson \& Dawson, 2003), SDMs is
still one of the effective methods for predicting species distribution,
assessing the risk of species extinction, planning and construction of
protected areas (Li et al., 2014; Preston, Rotenberry, Redak, \& Allen,
2010). SDMs are often used to predict the distribution of endangered or
narrow species, and there are relatively few studies on the prediction
of widespread species (Yang, Tang, \& Luo, 2020).

\emph{B. gargarizans} is one of the most widely distributed amphibians
in China, and compared with other amphibians, it has a relatively low
sensitivity to environmental factors, so it is a great sentinel species
for revealing amphibian responses to temperature changes in this study.
And the study of the genealogical geography of\emph{B. gargarizans}
shows that the differentiation of the subspecies of \emph{B.
gargarizans} is mainly caused by dominant dispersion rather than
vicariance, and showed a trend of dispersion from west to east (Fu,
Weadick, Zeng, Wang, \& Hu, 2005; Hu et al., 2007). Relevant studies
show that in the future, the ecological niche of a large number of
amphibians in China will migrate to the west or north, and the
fragmentation of the distribution area will be serious year by year.
(Duan, Kong, Huang, Varela, \& Ji, 2016). Although the genealogical
geography of \emph{B. gargarizans} is relatively thorough, its spatial
distribution prediction has not been carried out in depth. Starting from
the distribution area of \emph{B. gargarizans} , this study aims to
clarify the current and future trend of the \emph{B. gargarizans} 's
dispersion under the context of climate change, so as to provide certain
data support for the response model of amphibians to climate change.

The Maxent model was selected to simulate the current and future spatial
distribution of \emph{B. gargarizans} under the background of global
climate change and the migration of the centroid point, in order to
solve the following problems: i) to simulate the current and future
distribution of \emph{B. gargarizans} in China; ii) to predict the
change of suitable habitat of \emph{B. gargarizans} in different parts
of China under climate change; iii) to predict the distribution
barycenter migration of \emph{B. gargarizans} in China under climate
change, and reveal the mechanism of amphibian response to climate
change.

\textbf{2. Materials and methods}

\textbf{2.1} \textbf{Occurrence records of \emph{B. gargarizans.}}

The distribution data of\emph{B. gargarizans} used in this study were
mainly from two sources: i) Data accumulated through years of laboratory
fieldwork, 69 records were recorded; ii) and were obtained from the
Global Biodiversity Information Facility (GBIF), 287 records were
obtained. GBIF is a global scientific organisation that specialises in
providing basic data on biodiversity, the data is obtained by searching
for the keyword \emph{B. gargarizans} (visited in June 2020) (Jiang,
Zhang, Gao, Cai, \& Zhang, 2019). Distance between all species
distribution data \textgreater{}1Km by processing in ArcGis10.2
software. Finally, there are 274 pieces of data for this study (Table
S1).

\textbf{2.2 Environmental variables}

We collected 29 environmental variables for this study, which can be
roughly divided into four categories: bioclimate, topography, habitat,
and human impact, source of all variables show in Table S2. In order to
avoid model overfitting, all collected variables were processed by the
band collection statistics function in ArcGIS 10.2. Remove the
environment variables with high correlation (r\textgreater{}0.7). After
screening, Bio3, Bio5, Bio9, Bio14, Bio15, slope, Veg, NDVI, POP, GDP
and HI were used for this study. Bio3, Bio5, Bio9, Bio14, Bio15, slope,
Veg, NDVI, POP, GDP and HI. Since variables other than future
bioclimatic factors cannot be predicted and change relatively little
over a short period of time, we assume that they remain constant and use
their projections of the future distribution of \emph{B.
gargarizans}(Yang et al., 2020).

In order to reduce the uncertainty of model prediction, we chose four
internationally recognized GCMs (BCC-CSM1-1, HadGEM2-ES, IPSL-CM5A-LR,
MIROC5) and four RCPs (RCP2.6, 4.5, 6.0 and 8.0) published in the IPCC
fifth Assessment report for prediction in 2050 and 2070.

\textbf{2.} \textbf{3 Model optimization and parameter setting}

To establish the most stable and reliable model, we choose model
calibration, which includes regularization multiplier (0.5 to 6,
interval 0.5) and combination of basic feature classes and from one
different sets of all layers\textbf{.} The best parameters were selected
based on statistical significance (partial ROC), predictive ability
(miss rate E=5\%) and complexity level (AICc) by useing R package kuenm
(Cobos, Peterson, Barve, \& Osorio-Olvera, 2019). After the model
parameters were determined, chose the data output format of mathematical
logic to ensure that the output data was between 0 and 1. And the
jackknife method was selected to evaluate the contribution of
environmental variables. Then the distribution data
of\emph{B.gargarizans} were randomly divided into two groups: 25\% were
randomly selected as the test set, and the remaining 75\% were used as
the training set, repeat Bootstrap replicates 10 times in MaxEnt.3.4.1.
Other model parameters are selected by system default. Finally, the
predictive performance of the model was validated by using the area
under the receiver-operating characteristic (ROC) curve (AUC) and the
mean omission error. The higher the AUC value, the more accurate the
model performance. Normally, AUC \textgreater{} 0.9 represents excellent
prediction performance of the model. The ensemble threshold of model was
calculated according to the Maximizing Sensitivity and Specificity (MSS)
by using the dismo package in R. The MSS method is commonly used in
presence only kind of occurrence data(Hijmans, Phillips, Leathwick, \&
Elith, 2017; Liu, White, Newell, \& Pearson, 2013).

\textbf{2.4Distribution barycenter migration under different scenarios
of future global climate change}

We divided the study area into small grids of 0.1degx0.1deg to evaluate the
distribution barycenter migration of \emph{B. gargarizans} in China.
Assuming that our research area is composed of N small grids, the
proportion of the k\textsuperscript{th} grid is
G\textsubscript{K}=P\textsubscript{K}xS\textsubscript{K}, and
P\textsubscript{K} represents the probability of the species appearing
in the k\textsuperscript{th} small grid, and
S\textsubscript{K}represents the area of the k\textsuperscript{th} small
grid. The coordinates of each grid are obtained by Arcgis10.2. The
following formula is used to obtain the coordinates of the barycenter:

X= \(\frac{\sum_{k=1}^{k}G_{k}X_{k}}{\sum_{k=1}^{k}G_{k}}\); Y=\(\frac{\sum_{k=1}^{k}G_{k}Y_{k}}{\sum_{k=1}^{k}G_{k}}\)

Here, X and Y represent the longitude and latitude of the
k\textsuperscript{th} grid, respectively(He et al., 2011).

The results of four GCMs under four RCPS were averaged to obtain
P\textsubscript{f}, using \selectlanguage{english}^P(\selectlanguage{english}^P=P\textsubscript{f}-P\textsubscript{c},
P\textsubscript{c} represents the current distribution probability
of\emph{B. gargarizans} ) to evaluate the change trend of \emph{B.
gargarizans} distribution in China in 2050 and 2070:
-0.1\textless{}\selectlanguage{greek}Δ\selectlanguage{english}P[?]0.1 indicated that the habitat suitability changes
were not obvious. 0.1\textless{}\selectlanguage{greek}Δ\selectlanguage{english}P[?]0.3 , 0.3\textless{}\selectlanguage{greek}Δ\selectlanguage{english}P[?]0.5 and
0.5\textless{}\selectlanguage{greek}Δ\selectlanguage{english}P[?]0.7 represented slight, moderate and severe increasein
in habitat suitability, respectively. -0.3\textless{}\selectlanguage{greek}Δ\selectlanguage{english}P[?]-0.1,
-0.5\textless{}\selectlanguage{greek}Δ\selectlanguage{english}P[?]-0.3 and -0.7 \textless{}\selectlanguage{greek}Δ\selectlanguage{english}P[?] -0.5 represented slight,
moderate and severe declines in habitat suitability,
respectively\textbf{.}

\textbf{3.Results}

\textbf{3.1 Model optimization}

Through different combinations of feature classes and regularized
multipliers, the best model parameter is selected as among all
combinations (feature classes: linear, product, threshold, and hinge;
regularized multipliers=2), the parameters of all candidate models are
shown in Table S3.

In the simulation of the current and future distribution of \emph{B.
gargarizans} , the AUC value indicated that the simulation results of
all models reached an excellent level (0.9202[?] Training AUC[?]0.9306,
0.9001[?] Test AUC[?]0.9165, SD[?]0.014), all AUC values are presented in
Table S4

\textbf{3.2 Current distribution}

After running maxent model for 10 times, the average AUC value was
0.926\textbf{,} Standard deviation is 0.01 (Fig.1),. The model
validation results show that the model prediction results are excellent
and can be used in the following analysis\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Fig. 1. Receiver operating characteristic (ROC) curve and
average test AUC for accuracy analysis of Bufo gargarizan by MaxEnt
model under the current situation.}

We found that the most suitable habitats were mainly in the midwest and
the eastern coastal areas from current distribution of B.
gargarizans(Fig.2). The threshold value of MSS method is 0.33, and
according to the distribution probability of B. gargarizans, China is
divided into four grades: highly suitable habitat (HS, P[?]0.75); medium
suitable habitat (MS, 0.5[?]P\textless{}0.75); poorly suitable habitat(PS,
threshold[?]P\textless{}0.5); unsuitable habitat(US,
threshold\textless{}P[?]0), and the areas of these four types of habitats
were calculated in Arcgis, as shown in Table 1. Among them, the
unsuitable habitat is 7,216,700 square kilometers, accounting for 75.69
percent of Chinese areas.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Fig.2 Current distribution of MaxEnt models for \emph{B.
gargarizans} in China. The color scale from blue to red indicates the
habitat suitability value from 0 to 1.}

\textbf{Table1 Current habitat composition of \emph{B. gargarizans} in
China under MaxEnt model.HS: highly suitable habitat; MS: medium
suitable habitat; PS: poorly suitable habitat;US: unsuitable habitat.}\selectlanguage{english}
\begin{longtable}[]{@{}lllllllll@{}}
\toprule
Species & Area\textbackslash{}10\textsuperscript{4}
km\textsuperscript{2} & Area\textbackslash{}10\textsuperscript{4}
km\textsuperscript{2} & Area\textbackslash{}10\textsuperscript{4}
km\textsuperscript{2} & Area\textbackslash{}10\textsuperscript{4}
km\textsuperscript{2} & Percentage\textbackslash{}\% &
Percentage\textbackslash{}\% & Percentage\textbackslash{}\% &
Percentage\textbackslash{}\%\tabularnewline
\midrule
\endhead
& US & PS & MS & HS & US & PS & MS & HS\tabularnewline
\emph{B. gargarizans} & 721.67 & 134.16 & 67.63 & 30.35 & 75.69 & 14.07
& 7.1 & 3.18\tabularnewline
\bottomrule
\end{longtable}

\textbf{3.3 Contribution of environmental variables}

The Jackknife analyses results showed that the variables with high
contribution to the model were Bio9, Bio14, and POP, with a total
contribution of 70.8\%. Among them, the contribution rate of Bio9 is
38.4\%, that of Bio14 and POP is 16.2\%, detailed analysis results are
shown in the Table S5. Meanwhile, the range, mean value and optimal
range of each variable (the probability of species
existence\textgreater{}0.5) is shown in Table S6 and Fig.3. For
variables with higher contribution (Bio9, Bio14, POP), there were
significant differences in the mean values between the locations of
occurrence and absence.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Fig.3 The distribution of predictor bioclimatic and topographic
variables over the study area on presence and absence (background)
locations.}

\textbf{3.4 Analysis of future habitat changes}

The areas and proportions of the four habitats of \emph{B.
gargarizans}by 2050 and 2070 are shown in the Table S7. Compared with
the current distribution, the future suitable habitat of \emph{B.
gargarizans} will be reduced in 2050, and the ranges of highly suitable
habitat, medium suitable habitat, poorly suitable habitat are -1.3
\textasciitilde{} 0, -3.41 \textasciitilde{} -0.68, and -3.98
\textasciitilde{} 1.32, and the mean values are -0.72, -2.43, and -1.42
respectively. The suitable habitat of \emph{B. gargarizans} will shrink
further in 2070, and the ranges of highly suitable habitat, medium
suitable habitat, poorly suitable habitat are -1.42 \textasciitilde{}
-0.37, -3.71 \textasciitilde{} -1.59, and -3.94 \textasciitilde{} -0.03,
and the mean values are -0.93, -2.71, -2.04 respectively(Fig.4). And the
overall change trend and degree of suitable habitats of \emph{B.
gargarizans} in China is shown in Fig.5 and Table S8. In addition, the
potential distribution area of Chinese toad in the future shows a trend
of fragmentation, which is gradually serious over time, and this trend
is related to RCPs to a certain extent (RCPS8.5 is the most
serious)(Fig. S1)\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Fig.4 The four levels of suitable habitat fluctuation ranges
of\emph{B. gargarizans} under different future climate scenarios in
China. Red, orange, green, and blue represent increases (positive
number) or decreases (negative number) of the proportion of highly
suitable, medium suitable, poorly suitable, and unsuitable habitat,
respectively. BC, HE, IPand MC represent BCC-CSM1-1, HadGEM2-ES,
IPSL-CM5A-LR and MIROC5, respectively, which are different global
climate models. A is in 2050, B is in 2070.}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\selectlanguage{greek}\textbf{Fig.5 Variations in habitat suitability of Bufo gargarizans
under different future climate scenarios in China.
ΔP=P\textsubscript{f}--P\textsubscript{c}, P\textsubscript{c}
represented the probability value of suitability in the current
situation and P\textsubscript{f}represented the average probability
value of suitability under the future situation of four global climate
models. Gray (-0.1\textless{}ΔP[?]0.1) indicated that the habitat
suitability changes were not obvious. Blue (0.1\textless{}ΔP[?]0.3) , pink
(0.3\textless{}ΔP[?]0.5) and red(0.5\textless{}ΔP[?]0.7) represented slight,
moderate and severe increasein in habitat suitability, respectively.
Green (-0.3\textless{}ΔP[?]-0.1), Light green (-0.5\textless{}ΔP[?]-0.3) and
yellow (-0.7 \textless{}ΔP[?] -0.5) represented slight, moderate and
severe declines in habitat suitability, respectively.}

\textbf{3.5 Distribution barycenter migration}

To study the ecological barycenter migration under the context of
climate change, the centroid coordinates and migration distances of the
current and future(2050 and 2070) distribution of \emph{B.
gargarizans}were calculated(Table S9, Table S10). And the results of
distribution barycenter migration show that the distribution center of
\emph{B. gargarizans} will have moved to 107.05\selectlanguage{ngerman}°E \textasciitilde{}
108.94°E, 31.38°N \textasciitilde{} 32.74°N, and the distance of
barycenter migration will reach 38.47km \textasciitilde{} 249.45km. And
by 2070, the distribution of \emph{B. gargarizans} will continue to
change, and its distribution center will migrate to 107.00°E
\textasciitilde{} 108.75°E, 31.95°N \textasciitilde{} 32.37°E, the
migration distance is between 93.15 Km and 270.67 Km. Fig.5 and Fig.6
show that the suitable distribution area of \emph{B. gargarizans} will
decrease in the southeast, and gradually migrate and increase to the
midwest in China.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{Fig. 6 Barycenter transfer of \emph{Bufo gargarizans} in China.
The red hexagon represents the barycenter of a suitable habitat under
current scenarios. The black and green triangles represent the
barycenter of suitable habitats in 2050 and 2070, respectively, using
four global climate models.}

\textbf{4 Discussion}

\textbf{4.1 Maxent model and environmental variables}

Maxent model can connect the distribution data of species with
geographic and ecological information, and calculate the maximum entropy
value of species distribution, so as to predict the potential
distribution area of species(Giannakopoulos et al., 2019). Maxent model
can achieve a good prediction effect only by using limited species
distribution data. Therefore, it is often used to simulate the potential
distribution area of species, drift of potential distribution area,
suitability evaluation and other studies on global climate change
involving species distribution and genetic diversity in conservation
biology and invasion biology(Hannah et al., 2007; Mays, Hung, Shaner,
Denvir, \& Primerano, 2018; Prieto-Torres, Navarro-Sig\selectlanguage{ngerman}üenza,
Santiago-Alarcon, \& Rojas-Soto, 2016).

Environmental variables commonly used for Maxent model analysis include
climate variables, landscape variables such as topography and landform,
and the impact of human factors on species distribution is rarely
considered (Pearson, Raxworthy, Nakamura, \& Townsend, 2010; Penman,
Pike, Webb, \& Shine, 2010). However, the influence of human factors on
ecological environment is extensive and continuous, and its role cannot
be ignored(Li et al., 2014). From the contribution of the four
environmental factors selected in this study, the climate factor is
still the main factor affecting distribution of \emph{B. gargarizans} in
China. Among them, the contribution of Bio9 and Bio14 is relatively
large, the contribution was 38.4\% and 16.2\%, and the optimal
distribution range of \emph{B. gargarizans} was -1.85\selectlanguage{english} \textasciitilde{}
13.01\selectlanguage{english} and 3.59mm \textasciitilde{} 145.20mm, respectively. Secondly,
human influence also has a certain degree of contribution, and the
overall influence trend shows that, within a certain range, the
distribution probability of \emph{B. gargarizans} increases with the
enhancement of human influence (Table S5, Fig.3), this is caused by the
overlap of the ecological niche of the \emph{B. gargarizans} and the
area of human activities and some human infrastructure is more conducive
to the distribution of amphibians(Tytar, Sobolenko, Nekrasova, \&
Mezhzherin, 2015). Compared with other environmental factors, the
overall contribution of topographic and geomorphic factors is relatively
low, and see Fig.7 for details.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{Fig.7 Results of jackknife train of relative importance of
predictor variables for}B. gargarizans

\textbf{4.2 Future distribution change and barycenter migration ofB.
gargarizans}

Many studies suggest that global temperatures are rising and will
continue to do so for some time to come, with an expected increase of
0.6 in the 20th century and 1.4 \textasciitilde{} 5.8 by the 21th
century (Budikova, 2009; Mccarthy, Canziani, Leary, Dokken, \& White,
2007). The rising global temperature will inevitably affect the change
of species distribution, which will eventually lead to two results: one
is the reduction of species distribution, which will eventually lead to
extinction of species; the other is the migration of species to
previously undistributed areas, which will expand the distribution
(Bellard, Bertelsmeier, Leadley, Thuiller, \& Courchamp, 2012;
Morueta-Holme, Camilla, \& Svenning, 2010; Thomas et al., 2004).

In this study, under the background of climate change, the suitable
distribution area of \emph{B. gargarizans} will continue to decrease in
the next 30\textasciitilde{}50 years. In comparison to three other
global Change Scenarios, we can see the smallest reduction in suitable
habitat for \emph{B. gargarizans} under RCP2.6 Scenarios (The most
optimistic scenario is that greenhouse gas emissions peak between 2010
and 2020, then continue to decline, with a negative value by the end of
the century). And the suitable habitat of\emph{B. gargarizans} will be
further reduced Under the scenarios of RCP4.5 and RCP6.0(That represents
a decline after greenhouse gas emissions peaked in 2040 and 2080,
respectively). the suitable habitat of \emph{B. gargarizans} of China
decreased to the maximum extent under the condition of RCP8.5(Greenhouse
gas emissions continue to rise in this century). The results show that
with the increase of greenhouse gas emission concentration, the area of
suitable habitat of \emph{B. gargarizans} decreases more, on the other
hand, it suggests the importance of energy conservation and emission
reduction to the maintenance of ecological balance.

As show in Fig. 4 and Fig. S5, the suitable habitat of\emph{B.
gargarizans} will decrease in the south, southeast coast and central
part of China, and correspondingly increase in the central and western
part of China from now to 2050. And the trend will become more
pronounced by 2070. In the next 30 to 50 years, the distribution centers
of \emph{B. gargarizans}will gradually move from Yunyang County in
northeast Chongqing City to Wanyuan County, Tongjiang County and Xuanhan
County in northeast Sichuan Province and Langao County, Ziyang County,
Zhenba County, Nanzheng County and Chenggu County in southwest Shaanxi
Province. Compared with the present, the distribution focus has
gradually shifted to the northwest, and this trend has not weakened with
time. In addition, it can be seen from Figure 1 that Bio9 and Bio14, two
climatic factors that have a great influence on the distribution of
\emph{B. gargarizans} . are increasing from 2050 to 2070 compared with
the present, and the growth trend is more obvious in eastern and
southern China. And these are the areas where \emph{B. gargarizans} ,
will be less concentrated in the future. And the reason for this
phenomenon is that the average temperature and precipitation in these
regions are higher at present, and the further increase of temperature
and precipitation in the future will make the values of these factors
beyond the range suitable for the distribution of \emph{B. gargarizans}
(Table S6\textbf{)} .

Four GCMs under four RCPs were used in this study, and the results show
that under the global climate change, like other species that have been
reported, the distribution area of \emph{B. gargarizans} will migrate to
the high latitude and high altitude area (China's terrain is low in the
east)(Kelly \& Goulden, 2008; Parmesan \& Yohe, 2003). On the other
hand, this paper is more inclined to think that the current distribution
pattern of \emph{B. gargarizans} , is caused by the dispersion, because
from the perspective of the current distribution, there is no obvious
geographical isolation between the distribution areas of different
regions of \emph{B. gargarizans} ,, and the future distribution
of\emph{B. gargarizans} , showed that the future distribution of
\emph{B. gargarizans} , would gradually break up, and the geographical
separation would gradually form between different regions (Fig.2, Fig.
S1). The results of this study indicated that the distribution area of
\emph{B. gargarizans} would be reduced in the future, and the habitat
connectivity, resilience and gene exchange between the populations would
be reduced with the serious habitat fragmentation year by year, finally
leading to the gradual decline of \emph{B. gargarizans} population(Duan
et al., 2016). And it suggests that some amphibians, which are more
sensitive to environmental changes and narrower than \emph{B.
gargarizans} , the changes to their niches from climate change will be
even more dramatic and even devastating. Due to the influence of the
resolution of the environmental layer, the prediction results cannot be
very accurate. However, we still believe that this study will provide
theoretical support for the direction of conservation and ecological
management of \emph{B. gargarizans} , and it is believed that it can
provide a research direction for other amphibians to cope with climate
change.

\textbf{Data Accessibility:}

Occurrence records of \emph{B. gargarizans} were obtained from
laboratory survey and global biodiversity information platform (GBIF),
and the data finally used for analysis are listed in Table S1 in the
attachment.

Environment variables are derived from official websites such as
Worldclim, Resource and Environment Science and Data Center and Last of
the Wild Data. Please refer to the Table S2 in the attachment for the
specific download sites.

\textbf{References}

Bellard, C., Bertelsmeier, C., Leadley, P., Thuiller, W., \& Courchamp,
F. (2012). Impacts of climate change on the future of biodiversity.
Ecology Letters, 15(4), 365-377.

Bertrand, R., Lenoir, J., Piedallu, C., Riofrio-Dillon, G., Ruffray, P.,
Vidal, C., Gegout, J. C. (2011). Changes in plant community composition
lag behind climate warming in lowland forests. Nature, 479, 517-520.

Budikova, D. (2009). Role of Arctic sea ice in global atmospheric
circulation: A review. Global Planetary Change, 68(3), 149-163.

Cobos, M. E., Peterson, A. T., Barve, N., \& Osorio-Olvera, L. (2019).
Kuenm: An R package for detailed development of ecological niche models
using Maxent. PeerJ, 7(1), e6281.

Dawson, T. P., Jackson, S. T., House, J. I., Prentice, I. C., \& Mace,
G. M. (2011). Beyond redictions: Biodiversity conservation in a changing
climate. 332(6025), 53-58.

Duan, R. Y., Kong, X. Q., Huang, M. Y., Varela, S., \& Ji, X. (2016).
The potential effects of climate change on amphibian distribution, range
fragmentation and turnover in China. PeerJ, 4(10), e2185.

Elith, J., Leathwick, J. R. (2009). Species Distribution Models:
Ecological Explanation and Prediction Across Space and Time. Annual
Review of Ecology Evolution \& Systematics, 40(1), 677-697.

Fu, J. Z., Weadick, C. J., Zeng, X. M., Wang, Y. Z., \& Hu, Y. (2005).
Phylogeographic analysis of the Bufo gargarizans species complex: A
revisit. Molecular PhylogeneticsEvolution, 37(1), 202-213.

Giannakopoulos, A., Vasileiou, N. G. C., Gougoulis, D. A., Cripps, P.
J., Ioannidi, K. S., \& Chatzopoulos, D. C. (2019). Use of geographical
information system and ecological niche modelling for predicting
potential space distribution of subclinical mastitis in ewes. Veterinary
microbiology, 228, 119-129.

Guisan, A., \& Thuiller, W. (2010). Predicting species distribution:
offering more than simple habitat models. Ecology Letters, 8(9),
993-1009.

Hannah, L., Midgley, G., Andelman, S., Araujo, M., Hughes, G.,
Martinez-Meyer, E., . . . Williams, P. (2007). Protected area needs in a
changing climate. Front. Ecol. Environ., 5, 131-138.

He, Y. B., Chen, Y. Q., Tang, H. J., Yao, Y. M., Yang, P., \& Chen, Z.
X. (2011). Exploring spatial change and gravity center movement for
ecosystem services value using a spatially explicit ecosystem services
value index and gravity model. Environmental Monitoring Assessment, 175,
563-571.

Hijmans, R. J., Phillips, S., Leathwick, J., \& Elith, J. (2017). Dismo:
species distribution modeling. R Package Version, 1(4), 1-1.

Hu, Y. L., Wu, X. B., Jiang, Z. G., Yan, P., Su, X., \& Cao, S. Y.
(2007). Population Genetics and Phylogeography of Bufo gargarizans in
China. Biochemical Genetics, 45, 697-711.

Jiang, F., Zhang, J. J., Gao, H. M., Cai, Z. Y., \& Zhang, T. Z. (2019).
Musk deer (Moschus spp.) face redistribution to higher elevations and
latitudes under climate change in China. Science of The Total
Environment, 704, 135335.

Kelly, A. E., \& Goulden, M. L. (2008). Rapid shifts in plant
distribution with recent climate change. Proceedings of the National
Academy of ences of the United States of America, 105(33), 11823--11826.

Lambers, J. H. R. (2015). Ecology. Extinction risks from climate change.
Science, 348, 501-502.

Li, R., Xu, M., Wong, M. H. G., Qiu, S., Sheng, Q., Li, X., \& Song, Z.
(2014). Climate change-induced decline in bamboo habitats and species
diversity: implications for giant panda conservation Diversity \&
Distributions, 21(4), 379--391.

Liu, C., White, M., Newell, G., \& Pearson, R. (2013). Selecting
thresholds for the prediction of species occurrence with presence-only
data. Journal of Biogeography, 40(4), 778-789.

Mays, H. L., Hung, C. M., Shaner, P. J., Denvir, J., \& Primerano, D. A.
(2018). Genomic Analysis of Demographic History and Ecological Niche
Modeling in the Endangered Sumatran Rhinoceros Dicerorhinus sumatrensis.
Current Biology, 28(1), 70-76.

Mccarthy, J. J., Canziani, O. F., Leary, N. A., Dokken, D. J., \& White,
K. S. (2007). Climate Change 2001: Impacts, Adaptation and Vulnerability
-- Contribution of Working Group II to the Third Assessment Report of
IPCC. Contribution of Working Group II to the Third Assessment Report,
19(2), 81-111.

Morueta-Holme, N., Camilla, F., \& Svenning, J. C. (2010). Climate
change risks and conservation implications for a threatened small-range
mammal species. PLoS ONE, 5(4), e10360.

Nogues-Bravo, D., \& Rahbek, C. (2011). Communities Under Climate
Change. Science in China Series C: Life Sciences, 334, 1070-1071.

Parmesan, C., \& Yohe, G. (2003). A globally coherent fingerprint of
climate change impacts across natural systems. Nature, 421, 37--42.

Pearson, R. G., \& Dawson, T. P. (2003). Predicting the impacts of
climate change on the distribution of species: are bioclimate envelope
models useful? Global Ecology Biogeography, 12, 361-371.

Pearson, R. G., Raxworthy, C. J., Nakamura, M., \& Townsend, P. A.
(2010). ORIGINAL ARTICLE: Predicting species distributions from small
numbers of occurrence records: a test case using cryptic geckos in
Madagascar. Journal of Biogeography, 34(1), 102-117.

Penman, T. D., Pike, D. A., Webb, J. K., \& Shine, R. (2010). Predicting
the impact of climate change on Australia's most endangered snake,
Hoplocephalus bungaroides. Diversity \& Distributions, 16(1), 109-118.

Preston, K. L., Rotenberry, J. T., Redak, R. A., \& Allen, M. F. (2010).
Habitat shifts of endangered species under altered climate conditions:
importance of biotic interactions. Global Change Biology, 14(11),
2501-2515.

Prieto-Torres, D. A., Navarro-Sig\selectlanguage{ngerman}üenza, A. G., Santiago-Alarcon, D., \&
Rojas-Soto, O. R. (2016). Response of the endangered tropical dry
forests to climate change and the role of Mexican Protected Areas for
their conservation. Global Change Biology, 22(1), 364-379.

Sinervo, B., Mendez-de-la-Cruz, F., Miles, D. B., Heulin, B., Bastiaans,
E., \& Cruz, M. V. S. (2010). Erosion of Lizard Diversity by Climate
Change and Altered Thermal Niches. \%J Science. Science, 328, 894-899.

Stuart, S. N., Chanson, J. S., Cox, N. A., Young, B. E., Rodrigues, A.
S. L., Fischman, D. L., \& Waller, \selectlanguage{english}R. W. (2004). Status and Trends of
Amphibian Declines and Extinctions Worldwide. Science, 306(5702),
1783-1786.

Thomas, C. D., Cameron, A., Green, R. E., Bakkenes, M., Beaumont, L. J.,
Collingham, Y. C., . . . Hannah, L. (2004). Extinction risk from climate
change. Nature, 427(6970), 145--148.

Tytar, V., Sobolenko, L., Nekrasova, O., \& Mezhzherin, S. (2015). Using
Ecological Niche Modeling For Biodiversity Conservation Guidance In The
Western Podillya (Ukraine): Amphibians. Vestnik Zoologii, 49(6),
551-585.

Wake, D. B., \& Vredenburg, V. T. (2008). Are we in the midst of the
sixth mass extinction? A view from the world of amphibians. Proceedings
of the National Academy of Sciences of the United States of America,
64(1), 11466-11473.

Yang, C., Tang, S. H., \& Luo, Z. H. (2020). Distribution Changes of
Chinese Skink(Eumeces chinensis) in China: the Impacts of Global Climate
Change. Asian Herpetological Research, 11(2), 132-138.

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

