Climate Change is Likely to Alter Sheep and Goat Distributions in Mainland China

Climate change endangers food security worldwide, especially in developing countries. Livestock husbandry is one of the essential livelihoods for farmers and herders in remote arid and semiarid regions. However, it remains unclear how climate change will impact livestock husbandry in the future. This study collected sheep and goat distributions from the “gridded livestock of the world” (GLW) dataset for 1943 counties in Mainland China. Current climate data include precipitation and temperature from the National Meteorological Information Center (NMIC). We disentangled the effects of precipitation and temperature on current distributions of sheep and goats with the Bayesian Hierarchical Model by Integrated Nest Laplace Approximation (INLA). Further, we forecasted the potential sheep and goat distributions in 2030 and 2050 under Coupled Model Intercomparison Project (CMIP) scenarios. Our result showed that sheep distribution is significantly correlated with elevation, slope, market density, and highway distance, with absolute correlation coefficients ranging from 0.019 to 0.411. In addition to elevation, slope, and market density, goat distribution is also affected by gain production, with a correlation coefficient of 0.055. There is a dynamic correlation of temperature and precipitation with sheep and goat density. The sheep density distribution is predicted to increase in Northwest China, while the goat density distribution might increase in farming areas under climate change. Finally, this study suggests for the sheep and goat breeding industry to respond to climate change.


INTRODUCTION
Climate change threatens human society in the 21st century, especially on food security globally (Nemani, et al., 2003). China is the largest developing country in the world. Ensuring food provision safety is one of the most crucial issues the Chinese government faces (Zhang, et al., 2016). As reported, the per capita consumption of animal food, including mutton, milk, eggs, and beef, has increased by 160% in China, relative to the beginning of 'Reform and Opening up' (Han, et al., 2020). Climate change is another threat to China, where the warming rate is about twice the global average, and the flooding and droughts are becoming more frequent recently (Qian and Zhu, 2001). Livestock husbandry is important in food security and sustainable supply, especially under climate change. Therefore, it needs to disentangle, evaluate, and predict how climate change will affect livestock husbandry in the future.
Identifying and quantifying the climate effects on sheep and goat distributions is a hot topic in recent research. The standard Ricardian model is widely used to investigate climatic influences on agriculture (Mendelsohn, et al., 1994;Wang, et al., 2014;Gbetibouo and Hassan, 2005). Numerous studies have revealed how agricultural production responds to climate change (Thornton and Herrero, 2015), with much emphasis on decadal and interannual climate variability (Marumbwa, et al., 2019) and climate extremes (Moustakis, et al., 2021). Some studies have investigated the impacts of climate change on grassland degradation at the county scale (Gao, et al., 2010). A few studies also focused on the effects of extreme weather, such as snowstorms and droughts, on livestock production (Liu and Wang, 2012). However, it is hard to predict livestock distribution precisely under changing climate.
Several approaches, including linear regression (Leta and Mesele, 2014a), geographically weighted regression (GWR) (Velado-Alonso, et al., 2020), and Bayesian modeling (Rinella, et al., 2011), are developed for spatial simulation in agricultural studies. However, such data-driven methods cannot well address spatial correlations in agricultural systems (Homburger, et al., 2015). Failure to account for spatial uncertainty violates the assumption of independent and identically distributed data and may lead to biased model estimates (Anselin, 2001). Therefore, it is essential to explore uncertainties in spatial modeling to ensure that the final results are reliable (Rong, et al., 2017). Various uncertainties have been contributed to spatial models, input data, and analytical errors (Chen, et al., 2014;Heuvelink, et al., 2006;Nelson, et al., 2011). However, few studies highlight the uncertainty and potential drivers for sheep and goat distributions for sustainable policy-making (Illian, et al., 2013).
Bayesian statistics have been successfully used to explore the uncertainty in spatial modeling (Liang, et al., 2016;Williamson, et al., 2020). Simulation-based approaches such as Markov Chain Monte Carlo (MCMC) simulation are common to obtain the posterior distribution of model parameters (Taylor and Diggle, 2014). However, MCMC is computationally intensive and time-costing. The integrated nested Laplace approximations method (INLA) offers a simple way to compute complicated data across hierarchical scales (Huang, et al., 2017) and explore spatial correlation problems. That is because INLA can construct flexible fields by using stochastic partial differential equations (SPDE) to handle complicated-structured spatial data (Huang, et al., 2017).
This study proposes a hierarchical spatial framework fitted with SPDE and a dataset with sheep and goat distributions of 1943 counties in Mainland China. We aim to identify how climate drivers affect sheep and goat distributions along with social-economic variables. Further, we predict the possible distribution of sheep and goats in China under climate scenarios in 2035 and 2050. Our findings are expected to guide livestock husbandry development and help make climate adaptation and mitigation policies.

Study Area
This study covered 1943 counties in Mainland China. In 2002, per capita mutton consumption reached 1.91 kg in China, surpassing Organisation for Economic Co-operation and Development (OECD) countries ( Figure 1). The mutton industry development plays a vital role in improving residents' dietary structure, increasing farmers' income, and ensuring livestock production and supply. In 2020, there were 0.133 billion for goats and 0.173 billion for sheep in China (Supplementary Figure S1). Sheep mainly distributes in Inner Mongolia (26.92%) and Xinjiang Uygur Autonomous Region (22.83%), while goats in Inner Mongolia (12.26%), Sichuan (10.32%), Henan (10.22%), and Yunnan (8.31%) (see more details in Supplementary Figure S2).

Data Collection
Previous research (Gowane, et al., 2017;Sejian, et al., 2017) shows that climate, especially temperature and precipitation, significantly affects sheep and goat growth, development, and breeding. Therefore, we employed precipitation and temperature as primary explanatory variables for sheep and goat distributions across Mainland China. In addition, geomorphic factors including elevation and slope (Gong, et al., 2017), market factors in terms of market access, density, and influence indices (Mu, et al., 2017), and grain production, as well as grassland area (Aby, et al., 2014), were also considered ( Table 1).
The dependent variable was sheep and goat density in each county, which was derived from the "gridded livestock of the world" dataset (GLW) (Gilbert, et al., 2018). We validated the GLW data for each county with statistic records from Agricultural Bureaus at province and prefecture levels. The result shows that the GLW dataset matched well with statistic records (R-squared: 0.945, RMSE: 23.83, Supplementary Figure  S3, Supplementary Table S1). Then, we calculated the mean density of sheep and goat, respectively, for each of the 1943 counties in Mainland China, for further analyses.
Climate data were collected from National Meteorological Information Center (NMIC), China Meteorological Administration (CMA). Yearly temperature and precipitation records of over 800 meteorological stations across China were downloaded. We interpolated annual precipitation and temperature into rasters with a 1 km spatial resolution using the Kriging method with ArcGIS 10.8 (see isolines in Figure 2). Finally, we extracted climate data for each county to match with other explanatory variables.
Geomorphic factors are also seriously considered in sheep and goat distributions modeling because they can affect the combination of available water, light, and heat, and indirectly limit plant settlement, sheep growth, and breeding (Gong, et al., 2017;Sholagberu, et al., 2017). Exactly, we calculated each county's mean elevation and slope based on digital elevation model (DEM) data from moderate resolution imaging spectroradiometer (MODIS). Land cover can also affect sheep and goat distributions (Qian, et al., 2012). So, we also calculated the grassland area for each county, according to China's National Land Use/Cover Dataset (NLCD).
Market factors are also important in livestock distribution (Mu, et al., 2017;Ou and Mendelsohn, 2017). We used the dataset of Verburg (2011), which includes market access, market density, and market influence, to assess the market influence on sheep and goat distributions (Verburg, et al., 2011). Such data can jointly represent market strength and accessibility and are widely used to Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 748734 highlight market influences in recent studies (Sloat, et al., 2018). Finally, we resampled and averaged the three market indices and for each of the 1943 counties. Further, the mean distance of highways within each county was calculated as another economic-social factor. Eco-economic policies can significantly affect sheep and goat husbandry (Hu, et al., 2019). We collected all policies related to sheep and goat husbandry proclaimed by the government at all levels. We evaluated different weights to calculate the policy index according to the effectiveness of differentlevel policies in different counties (Formula 1).
Finally, we employed the grain production for each county to assess the feed supply ability based on a previous study. Thus, we calculated the county spatial mean for all variables we collected and correlated into the county by the ArcGIS 10.8.
This study also predicted sheep and goat distributions under different climate change scenarios. The sixth phase of Coupled Model Intercomparison Project (CMIP6) was used in this study. Exactly, we employed the future scenarios, including SSP1-2.6, SSP2-4.5, and SSP5-8.5 (Supplementary Table S2). The SSP1-2.6 scenario represented the low ending range of future scenarios measured by its radiative forcing pathway. This scenario will produce a multi-model mean of significantly less than 2.0°C warming by 2,100; thus, it could support the 2.0°C temperature rise target study. The SSP2-4.5 scenario considers  a medium stabilization scenario, while the SSP5-8.5 was a scenario that stabilizes radiative forcing at 8.5 Wm2 in 2,100 and is regarded as a high radiative forcing scenario.

Analyses
We employed the Bayesian hierarchical model to estimate the spatial distributions of sheep and goats, respectively, across 1943 counties in China. Parameter estimates and marginal posterior probability distributions were obtained by using the integrated nest Laplace approximation (INLA) approach. In this study, weakly informative penalized-complexity priors were generated for all regression coefficients (fixed-effect parameters) and hyperparameters. Sheep and goats distributed in nearby regions were often exposed to similar social environments and climate conditions. We used the INLA-SPDE approach to deal with spatial covariance among explanatory variables. The residual errors may reflect the influences of unmeasured or unmeasurable predictors that vary across space. The potential Bayesian hierarchical spatial framework (Morris, et al., 2019) is used to calculate sheep and goat distributions below (Formula 2). Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 748734 4 where the response variable Y s is the density of sheep or goat at the county s, D(s) is a linear combination of fixed covariates, and β is the corresponding coefficient vector. R(s) denotes the spatial random fields which correlated with climate factors. And ξ(s) is the spatial measurement effect by SPDE. In our study, all the covariates are measured by the logarithm form. It is noted that the random fields ξ(s) could be adopted to compensate unavailable or unobserved external factors not included in our model, and R(s) could also be used to capture spatial correlation components.
The spatial effect was assessed by fitting the model with and without random effects. An information criterion based on the most significant length of the triangle edge was used to solve the trade-off. Exactly, we employed the Deviance Information Criterion (DIC) and Watanabe-Akaike Information Criterion (WAIC) to compare the performance of various candidate models (Sutanto et al., 2021). The model with the smallest value of WAIC is preferred as it achieved the most combination of fitness and parsimony.
It is common to split data into training and validation sets, especially with predictive modeling purposes applied in empirical analysis. This study employs K-fold cross-validation to determine if the model correctly estimates the observed data (Wong and Yeh, 2019). The procedure has a single parameter called k that refers to the number of groups that a given data sample is grouping. We conducted a 10-fold cross-validation experiment by randomly removing 10% of the county in 10 independent estimates, respectively (Supplementary Material S5). These accuracy criteria reflect that the prediction in this study based on the INLA approach is maintained accuracy in both sheep and goat distributions estimation. Besides, the Pearson's correlation coefficient and concordance correlation coefficient indicate that the prediction in this study is valid. The outcomes show that the model predicts better than chance alone ( Table 2).
These accuracy criteria reflect that the prediction in this study based on the INLA approach is maintained accuracy in both sheep and goat distributions estimation. Besides, the Pearson's correlation coefficient and concordance correlation coefficient indicate that the prediction in this study is valid. Further, an Area Under Curve (AUC) > 0.5 indicates that the model predicts better than chance alone. The AUC for the validation dataset is 0.792 for goats and 0.917 for sheep.

Fixed Effect Covariates
The mean posterior coefficient of the fixed effect covariates is presented in log-odds and signifies the estimated response to a one standard deviation change in the sheep and goat distributions when all other variables are held constant. Fixed effect coefficients and 95% credible intervals for covariates are provided in Table 3. Both sheep and goat distributions were influenced significantly by geomorphic factors. As expected, a significant association is observed between elevation and both sheep and goat distributions in China. For a one standard deviation increase in elevation, the expected change in goat density log odds is 0.291, and sheep density log odds is 0.223. Further, for a one standard deviation decrease in the average slope of each county, the expected change in goat and sheep density log odds is 0.321 and 0.411, respectively.
Both goat and sheep density distribution were also positively affected by market density (Tables 3). An increase in one standard deviation in market density might increase 0.038 and 0.036 in logodds for goat and sheep density, respectively. The developed socialeconomic condition might restrict sheep distribution by comparative industrial advantage for the county with a higher market influence index. Grain production has a significant positive connection with goat density, while for one standard deviation increase, the expected change in goat density log odds is 0.055. For sheep husbandry, sheep breeding activity might be more suitable for areas where have lower highway density. For one standard deviation decrease, the expected change in sheep density log odds is 0.019.

Random Effects and Hyperparameters
The spatial model component shows how the spatial and climate affect odds of sheep and goat density distribution. The spatial random effect indicates a significant amount of spatial variation in both sheep and goat density distribution (Figure 3). The spatial random effect indicates a significant spatial variation in the sheep and goat density across China. Higher random effect located in the North and North-west of China for the sheep density. Meanwhile, an intensive spatial effect focused on the farming area, including the North China Plain and southwest China, for goat density.
The marginal distribution of scale parameter, variance parameter of the random field, and the practical range are shown in Table 4. The variance of the spatial effect showed a relatively wide posterior distribution, indicating that the variability in sheep and goat density to location is high. The posterior mean of the spatial correlation range was 618.69 km for sheep and 642.89 km for goats. These ranges indicate the approximate distance between counties where sheep and goat density distribution could be considered correlated. The ratio of range parameter indicates the moderation of spatial autocorrelation of sheep and goat density. Further, the uncertainties of these parameters are minor, indicated by their standard deviation and quantiles in Table 4. Finally, Figure 4 reflects the impact of precipitation on sheep and goat density distribution in China. The precipitation might affect sheep distribution dynamically. The summit effect concentrates on around 400 mm annual precipitations, which means that the sheep husbandry industry in China might have an optimal area correlated with precipitation. This could explain the spatial sheep distribution centralized on the North China Plain and Inner Mongolia Plateau and the other regions with proper climate environment. There is a more sensitive effect for temperature on goats distribution than sheep (Figure 4), and a higher annual cumulative temperature could be more favorable for goats relative to sheep husbandry.

Climate Effect on Sheep and Goat Distributions
Climate factors, including precipitation and temperature, are the key to analyzing and predicting sheep and goat density distribution, an association that has been observed in many studies to date (Anya, et al., 2013). A vast amount of literature has concentrated on examining the mutual influences between climate environment and animal husbandry behavior (Batsuuri and Wang, 2017;Duricic, et al., 2019;Petit and Boujenane, 2018;Sloat, et al., 2018). Several mechanisms operate concurrently to the advantage of sheep and goat distributions at moderate climate conditions (Stanimirova, et al., 2019). Precipitation affects sheep husbandry dynamic; this is consistent with studies on a country level, and global scale agrees with animal experiments (Duricic, et al., 2019). A potential explanation is that grassland and feed supply are affected by precipitation (Derbile and Kasei, 2012). The specific precipitation conduct suitable moisture for sheep growth, decreasing the risk of epidemic disease (Rinaldi, et al., 2015). High humidity might increase the probability of blight and other diseases in sheep growth (Qamar, et al., 2009). The optimal annual precipitation is about 400 mm, which could hold the moisture for livestock grazing. The East of Inner Mongolia and North-East China might decrease the sheep husbandry scale upon climate change (Wang, et al., 2021).
For goat husbandry, a higher annual accumulated temperature could positively affect goat production in a particular range. The most suitable environment temperature range for goats is 5-25°C Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 748734 7 (Joyce, et al., 1966). When the environment temperature exceeds 25°C, the goat's body temperature will rise, and their forage intake will decrease (Moustakis, et al., 2021). The high temperature will adversely affect goat growth, such as heat stress and disease risk (Godde et al., 2020;Havlik et al., 2014).
Climate is forecasted to be warmer and wetter with shifting precipitation spatiotemporally soon in China. Therefore, it becomes urgent to analyze the variation of sheep distribution response to precipitation change and formulate sustainable and adaptive policies on animal products industry regulation. Our study confirmed that such a variety of sheep and goat distributions responses to climate change is heterogeneous at the county level. Sheep distribution sensitivity varies along with precipitation values (Castillo, et al., 2021;Epps, et al., 2004). The effect of precipitation on sheep distribution is positive around precipitation of 400 mm. This might be because ecosystems in arid and semiarid areas are fragile.
Further, an excessively high rainfall scale would not help sheep grow and feed. In the spatial realm, sheep husbandry is located in the North China Plain and Inner Mongolia Plateau, within precipitation of 400 mm. The sheep husbandry might be adversely affected by climate change in the north of Qinghai and the east of Inner Mongolia, likely due to increasing precipitation trends ( Figure 5). However, the areas of northwest China might benefit sheep husbandry because of the increasing trend of precipitations, where there are significant advances in soil moisture and livestock grazing. The forecasting of goats distribution alter in China is more stable compared with sheep husbandry. There is no significant change in pasture areas in China. The increased trend concentrates on the areas which maintain the breeds of advantages, including Jiangsu, Hunan, and Liaoning province. The adaptability of the goat could explain this trend. Further, goat husbandry might be developing in the south of China, which mountainous region upon grassland.

Importance of Geomorphic and Social-Economic Factors in Sheep and Goat Distributions
Classic industrial economic models posit a substantial assumption in agricultural spatial distribution with livestock distribution with traffic, social-economic, and natural conditions (Behnke, et al., 2008). As expected, a significant association is observed between traffic conditions and sheep distribution in China. A higher average highway density not Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 748734 only means better transportation conditions but also presents a higher population density. The sheep husbandry has to move into areas far away from the city to maintain the profit. This finding agrees with the classic agricultural location theory (Dawkins, 2003). Compared with other agricultural industries, sheep husbandry might need to be located in areas with an inexpensive cost for higher marginal profits. Grain production posts a strong positive effect on goat husbandry. A potential explanation is the feeding structure of China, which is focused on grain straw (Liu, et al., 2008). This finding might partially explain the concentrated trend of goat husbandry in the North China Plain. The household could feed the goats with straw, as the by-product of grain crops, with a lower marginal cost. Second, our fitting model partially agrees with (Leta and Mesele, 2014b) that market influence is a critical accelerate factor for sheep husbandry. This study found that market density is the key market factor for both sheep and goat distributions in China's mainland, which the index combines the effect of market accessibility, per-GDP, and population density. This might explain the concentration on the North China Plain in sheep and goat distributions for Shandong and Henan province in China. Because these areas were located in high population density areas, they require a large scale of sheep consumption demand.
Third, we found that geographic factors are critical for both sheep and goat husbandry in China. Geomorphic features, including slope and elevation, affect sheep and goat husbandry significantly, being consistent with (Dwyer, 2009). High elevation areas with flat terrain are more suitable for sheep and goat husbandry in China. The sheep and goat both require suitable geographic and geomorphic conditions. Further, it is noticeable that the fixed effect of slope on goats is lower than sheep, which means that sheep husbandry might be more sensitive to geomorphic change. It could be supported by recent animal science literature that compared to sheep, goat is more adaptive in mountainous and hilly regions (Raoult, et al., 2021).

Policy Implication and Limitation
Differences in climate effect on sheep and goat husbandry call for resilience-enhancing actions that are region-specific and contextspecific and guided by principles of equity and fairness. Most urgent are actions and investment for counties that faced compound climate risk across livestock husbandry outcomes, requiring a transformative change to reduce the negative effect of climate change. For example, there might be a tremendous demand for fodder with the alter of sheep and goat density distribution under climate change. The government needs to balance grain security and protein supply to ensure the fodder supply and promote the livestock husbandry industry. Also, sheep distribution prediction under different climate scenarios indicates the vast increase of sheep breeding stock in the northwest of China, with many livestock pollutants to challenge the local ecological carrying capacity. We need to focus on the balance between the ecology effect and livestock pressure. Through the accurate acquisition of livestock waste information, we can adjust herd sizes, use precise grazing, and apply knowledge and lessons learned from ecological livestock breeding. Further, the possible distribution of sheep and goats demands a better mechanism for preventing and controlling animal epidemics and increasing the investment for animal epidemics monitor and warning, enhancing the technology application on animal epidemics prevention. Finally, effort should also be made to promote the breeding management ability of the livestock breeding operators, improve people's livelihoods, and coordinate the development of regional economies. Although we have developed a reliable model for predicting sheep and goat distributions under climate change, this study has several limitations. The GLW dataset in China still limits our predictive ability in small areas. As a result, some counties with large acreage may be underrepresented vs. counties with tiny areas in some variables. This study incorporated several environmental and social-economic covariates known to be associated with sheep and goat distributions. However, our results showed a strong spatial effect, the cause of which is undetermined in the present analysis. Future studies could add more other social and cultural determinants. The market effect is the key to livestock husbandry and grazing distribution. Therefore, we limited the model to data that could be readily available and most relevant to sheep and goat distributions.

CONCLUSION
We investigated the sheep and goat distributions in China and estimated the spatial factor in sheep and goat distributions using the Bayesian hierarchical model with the INLA approach. The findings of this study reveal that the sheep density distribution is highly correlated with elevation, slope, market density, and Highway distance, with absolute correlation coefficients ranging from 0.019 to 0.411. In addition to elevation, slope, and market density, goat distribution is also affected by gain production, with correlation coefficient of 0.055. For precipitation factors to sheep distributions, with a non-linear and dynamic effect, proper annual precipitation around 400 mm would produce a higher positive effect with sheep production. For temperature factors to goat distributions, a higher accumulated temperature means a positive impact on goat husbandry. To predict sheep distribution under climate change by CMIP6, we found the sheep density distribution might be an increase in the northwest of China, and the goat density distribution might increase in farming areas. Based on our comprehensive analysis of sheep and goat distributions in China, in light of the climate effect on sheep and goat density distribution, this study offers several policy suggestions for the sheep and goat breeding industry to respond to climate change.
Further, we forecasted the potential sheep and goat distributions in 2030 and 2050 under Coupled Model Intercomparison Project (CMIP) scenarios. Our result showed that sheep distribution is significantly correlated with elevation, slope, market density, and highway distance, with absolute correlation coefficients ranging from 0.019 to 0.411. In addition to elevation, slope, and market density, goat distribution is also affected by gain production, with correlation coefficient of 0.055. There are dynamic correlation of temperature and precipitation with sheep and goat density. The sheep density distribution is predicted to increase in Northwest China, while the goat density distribution might increase in farming areas under climate change. Finally, this study suggests for the sheep and goat breeding industry to respond to climate change.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
YZ and GW developed the project, gathered data, and conducted the initial round of analysis including determining variable used. SZ and CH provided editing and technical advice during the initial analysis round. YZ and YZ conducted the final analysis. All authors contributed equally to the writing of the final manuscript.

ACKNOWLEDGMENTS
The first author is grateful to Qionghua Li for answering all my questions promptly and wish you could accompany me on my academic journey.