Data
We reviewed PPS reporting rates of AMR in healthy animals and animal food products in China between 2000 and 2019 (Supplementary Text S1). We focused on three common food animal species: chicken, pigs and cattle. Dairy cattle and meat cattle were pooled in this study, in consistency with the categorization adopted in the maps of livestock created by the Food and Agriculture Organization27. The review focused on four common foodborne bacteria: E. coli, non-typhoidal Salmonella, S. aureus and Campylobacter. We recorded resistance rates reported in PPSs, defined as the percentage of isolates tested resistant to an antimicrobial compound. In addition, we extracted the anatomical therapeutic chemical classification codes of the drugs tested, the year of publication, the guidelines used for susceptibility testing, the latitude and longitude of sampling sites, the number of samples collected and the host animals. We recorded sample types for each survey, including live animals, slaughtered animals, animal products and faecal samples. Each sample was taken from one animal or animal product. These sample types were pooled in the current analysis. In total, 10,747 rates of AMR were extracted from 446 surveys (Supplementary Fig. 8), including 318 surveys from China’s National Knowledge Infrastructure (CNKI), the leading Chinese-language academic search engine. All data extracted in the review are available at https://resistancebank.org.
Two steps were taken to ensure comparability of the resistance rates extracted from the surveys. First, the panel of drug–bacteria combinations extracted from each survey was that recommended for susceptibility testing by the World Health Organization Advisory Group on Integrated Surveillance of Antimicrobial Resistance40. This resulted in the extraction of 6,295 resistance rates for 76 drug–bacteria combinations. Second, resistance rates were harmonized using a methodology4 accounting for potential variations in the clinical breakpoints used for antimicrobial susceptibility testing (Supplementary Text S1). There are two major families of methods used for susceptibility testing in this dataset: diffusion methods (for example, disc diffusion) and dilution methods (for example, broth dilution). Previous works have shown good agreement between the two approaches in measuring resistance in foodborne bacteria4,46. For each family of methods, variations of breakpoints may result from differences between laboratory guidelines systems (European Committee on Antimicrobial Susceptibility Testing vs Clinical and Laboratory Standards Institute), or from variations over time of clinical breakpoints within a laboratory guidelines system (Clinical and Laboratory Standards Institute or European Committee on Antimicrobial Susceptibility Testing). Here we accounted for both situations using distributions of minimum inhibitory concentrations and inhibition zones obtained from eucast.org (Supplementary Text S1).
Trends in AMR
We defined a composite metric of AMR to summarize trends in resistance across multiple drugs and bacterial species. For each survey, we calculated the proportion of antimicrobial compounds with resistance higher than 50% (P50). For each animal–bacteria combination, we assessed the significance of the temporal trends of P50 between 2000 to 2019 using a logistic regression model, weighted by the log10-transformed number of samples in each survey.
For each bacteria–drug (antimicrobial class) combination, we estimated prevalence of resistance by calculating a curve of the distribution of resistance rates across all surveys (Fig. 2). The analysis was conducted for surveys published between 2000 and 2009, and between 2010 and 2019, respectively. The distribution was estimated at 100 equally spaced intervals from resistance rates of 0% to 100%, using kernel density estimation. We used the centre of mass of the density distribution to estimate prevalence of resistance. The calculation was conducted for six animal–bacteria combinations. This included E. coli in chicken, pigs and cattle, Salmonella in chicken and pigs, and S. aureus in cattle. The remaining animal–bacteria combinations were excluded due to limited sample size, only represented in 32 out of 446 PPSs. The analysis was restricted to antimicrobial classes represented by at least 10 resistance rates. In addition, we estimated the association between resistance rates and the ease of obtaining antimicrobials from the market, using data from online stores (Supplementary Text S3).
Geospatial modelling
We interpolated P50 values from the survey locations to create a map of P50 at a resolution of 10 × 10 km across China. The approach followed a two-step procedure47. In step 1, three ‘child models’ were trained using four-fold spatial cross-validation to quantify the relation between P50 and environmental and anthropogenic covariates (Supplementary Text S2 and Supplementary Table 1). In step 2, the predictions of the child models were stacked using universal kriging (Supplementary Text S2). This approach combined the ability of the child models to capture interactions and non-linear relationships between P50 and environmental and anthropogenic covariates, as well as the ability to account for spatial autocorrelation in the distribution of P50.
The outputs of the two-step procedure were a map of P50 (Fig. 3) and a map of uncertainty on the P50 predictions (Supplementary Fig. 9 and Supplementary Text S2). The overall accuracy of the geospatial model was evaluated using the area under the AUC. The contribution of each covariate was evaluated by permuting sequentially all covariates, and calculating the reduction in AUC compared with a full model including all covariates (Supplementary Fig. 4). The administrative boundaries used in all maps were obtained from the Global Administrative Areas database (http://www.gadm.org).
Identifying (optimal) locations for future surveys on AMR
We identified the locations of 50 hypothetical new surveys—the rounded average number of surveys conducted per year (54 surveys per year) between 2014 and 2019 in China. The location of each new survey was determined recursively such that it minimized the overall uncertainty levels on the geographical trends in AMR across the country. This process took into account the locations of existing surveys and the location of each additional hypothetical survey. The objective of this approach was to maximize gain in information about AMR given the resource invested in conducting surveys.
The map of uncertainty consisted of the variance in the child model predictions Var(PBRT,PLASSO–GLM,PFFNN) (step 1) across 10 Monte Carlo simulations, where PBRT, PLASSO–GLM, and PFFNN were the predictions of P50 using boosted regression trees, logistic regression with LASSO regularization, and feed-forward neural network, and the kriging variance VarK (step 2):Vartotal = Var(PBRT,PLASSO–GLM,PFFNN) + VarK
In this study, the location of hypothetical surveys was solely based on VarK, instead of the sum of both terms. This approach was preferred because including both terms would have required to hypothesize P50 values associated with surveys to be conducted in the future, adding an additional source of uncertainty that cannot be quantified. In any case, the uncertainty attributable to VarK was 4.1 times the Var(PBRT,PLASSO–GLM,PFFNN) (Supplementary Text S2).
The allocation of new surveys was based on a map of ‘necessity for additional surveillance’ (NS), defined as:NS = VarK × Wwhere VarK reflects the uncertainty of the spatial interpolation, and W is the log10-transformed population density of humans48, animals27 in total, and in chicken, pigs and cattle, separately, which reflected exposure (Supplementary Fig. 10). Animal population density was calculated here as the sum of population-corrected units of pigs, chicken and cattle, using methods described by Van Boeckel et al.7. We adjusted the values of W such that its density distribution equals that of VarK. Concretely, for each pixel i, we calculated the quantile of Wi on the map of W, and replaced the value by the corresponding value of VarK at the same quantile. VarK and W were both standardized to range [0,1], thus giving each term equal weight in the need for surveillance.
Four approaches were used to distribute 50 surveys across China based on the map of NS. The reduction in uncertainty on AMR level associated with each of the four spatial configurations of the hypothetical surveys was evaluated by calculating the reduction in the mean values of NS across 7,857 possible pixels on the map of China.
First, we used a ‘greedy’ approach where all possible locations for additional surveys were tested. Concretely, the first hypothetical survey was placed at each of the 7,857 possible pixel locations, and a revised map of NS(+1 survey) was calculated for each of the placements. The survey was eventually placed in the pixel that led to the largest reduction in NS(+1 survey). The map of NS was then revised to account for the reduction in uncertainty in the neighbourhood of the new survey. The process was repeated recursively for the next hypothetical surveys (2nd–50th). This approach, by definition, yields the optimal set of locations to reduce uncertainty, but it also bears a considerable computational burden, because every possible location is tested (Npixels = 7,857) by the geospatial model for each hypothetical survey.
The second approach developed was a computational approximation to the greedy approach, hereafter referred to as the ‘overlap approach’. This approach exploits a key feature of the kriging procedure: the decrease of the kriging variance (VarK) with increasing proximity to existing survey locations. Each additional survey reduces the variance of the geospatial model at its own location, but also in its surrounding area (Supplementary Fig. 11). The ‘overlap approach’ selects an optimal set of locations that reflect a compromise between high local NS and distance to other surveys. It iteratively selects new locations based on the highest local NS penalized by the degree of overlap between the hypothetical new surveys and existing surveys (Supplementary Fig. 12). The first survey was placed at the location Xp,Yp with the highest local NS (Supplementary Fig. 12, part 1). Then the value of NS at each pixel location Xi,Yi was recalculated as ({mathrm{NS}}_{{(+1,{mathrm{survey}})}X_i,,Y_i}={mathrm{NS}}_{X_i,,Y_i}times(1-{mathrm{overlap}},{mathrm{area}}/{mathrm{neighborhood}},{mathrm{area}})) (Supplementary Fig. 12, Part 2), where the neighbourhood area was the circular area of decreased kriging variance around a new survey, and its radius was the distance until which NS decreased due to this new survey; the ‘overlap area’ is the shared area of the neighbourhoods of location Xp,Yp and of location Xi,Yi. The radius of the neighbourhood was determined using a sensitivity analysis, optimized by approximate Bayesian computation (sequential Monte Carlo)49 (Supplementary Text S5). The optimal neighbourhood radius was chosen such that it minimizes reduction in NS across all pixels. The procedure (Supplementary Fig. 12, parts 1 and 2) was repeated recursively for the hypothetical surveys (2nd–50th).
The third approach tested consisted of distributing surveys equally between provinces to reflect a common approach to disease surveillance based on equal allocation of resources between administrative entities. Twenty-two provinces with the highest human population were assigned two surveys, and the remaining six provinces were assigned one survey per province. The exact location of each survey was randomly selected inside a province. Finally, all approaches were compared with the fourth approach (the random approach) as a ‘null model’, in which the 50 hypothetical surveys were located randomly across the country without any geographic weighting criteria. The reduction in NS associated with the third and fourth approaches, which was compared to the greedy approach and overlap approach, was the average over 50 simulations.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
