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).

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 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 different-level policies in different counties (Formula 1). Policy   Index = ρ 1 Country + ρ 2 Province + ρ 3 County (1) where ρ 1 , ρ 2 , and ρ 3 indicate the weight of different policy levels from the county, province, to the whole country.

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.

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). Y s = possion ( θ s ) θ s =   β D ( s ) + R ( s ) + ξ ( s ) + γ d s + u s (2) 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.

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 log-odds for goat and sheep density, respectively. The developed social-economic 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.

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.

Figures 5, 6 presents the predicted density distribution change of sheep and goats in 2035 and 2050 in China under different climate scenarios. Figure 5 highlights the area at sheep density distribution change of climate change in various climate change scenarios. The maps show where the model predicts sheep density distribution change over the next 15 and 30 years. Climate change would intensify the sheep distribution in Northwest China, especially Xingjiang, Gansu, and Inner Mongolia (Figure 5). The sheep might be decreased based on the climate change under upon North-East of Qinghai province, and sector of Huabei areas to some extent. Compared with the altered sheep distribution, the goat might be more stable under climate change. The goat would concentrate on the North China Plain and South-West area of China. For goat husbandry (Figure 6), the increased space is located in Liaoning province, Hunan province, and Jiangsu province where local goat breeding, including Liuyang Black goat, White fuzz goat of Liaoning, and Haimen goats, is developed.

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 (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.

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 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).

Differences in climate effect on sheep and goat husbandry call for resilience-enhancing actions that are region-specific and context-specific 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.