Data acquisition
Our analysis was based on open access crop production data from the United Nation’s Food and Agricultural Organization (FAO) spanning from 1961 to 201718. We extracted data on area harvested (in ha) for 339 FAO-defined crop groups being grown in all UN-recognized countries. Since our research centred on understanding, quantifying, and mapping changes in crop diversity in current agricultural lands, countries that cease to exist (e.g., Yugoslavia) were not included in our analysis, resulting in data for 201 countries (Table S1). Prior to analyses, we adjusted certain crop group listings following our previous analyses of global changes in crop diversity8. Specifically, “Cottonlint” and “Cottonseed” were duplicated in our dataset and were therefore compiled as “Seedcotton”, while “Palmkernels” were renamed as “Oilpalmfruit.” Additionally, “Fruitpomenes”, “Fruitstonenes”, and “Grainmixed” were removed from analysis since these crop groupings are not associated with any specific crop species in the FAO database18. Finally, “Mushroomsandtruffles” were removed since it relates to non-plant species, and “Coir” was removed because it is a plant by-product.
Changes in crop richness over time
All statistical analyses were performed using R version 3.3.3 statistical software (R Foundation for Statistical Computing, Vienna, Austria). The initial step in our analysis was to calculate both crop richness and evenness for each country, at each individual year, using the vegan R package38. Based on these datasets, we then used the analytical framework developed by8 to evaluate how crop species richness and evenness have changed in each individual country across its entire data range.
Specifically, in their analysis Martin et al.8 found that piecewise linear regression models provided the strongest descriptions of crop species richness change over time, across 21 of 22 FAO-defined regions globally. We therefore followed this approach by fitting a piecewise linear regression model for each country individually, that predicts changes in species richness over time. Piecewise model fitting was a two-step process, whereby for each country we first fit a linear regression model of the form:
$$S = a + left( {b times {text{year}}} right)$$
(1)
where a is the intercept and b represents the rate of change in crop group richness (S) through time. This linear model (Eq. 1) was then used as the basis of a piecewise linear regression model, which was fitted in order to estimate breakpoints in the relationship between S and year. Specifically, piecewise models were fit using the segmented function in the segmented R package39, and were of the form:
$$S = a + bleft( {{text{year}}} right) + left( {left( {c({text{year}} -uppsi _{1} } right) times Ileft( {{text{year}} >uppsi _{1} } right)} right) + left( {dleft( {{text{year}} -uppsi _{2} } right) times Ileft( {{text{year}} >uppsi _{2} } right)} right)$$
(2)
where a is as in Eq. (1), and b represents the slope of the S-year relationship prior to the first breakpoint (ψ1). Here, c represents the difference in the slope of the S-year relationship between the first and second piecewise model segments; the c parameter therefore applies only when the first conditional indicator function (denoted by “I”) is true. Similarly, d represents the difference in slopes for the S-year relationship between the first, second, and third segments, which only applies when the second conditional indicator function is true. In sum, the slope of the relationship between S and year is equal to b prior to the ψ1, is equal to b + c between ψ1and ψ2, and is equal to b + c + d after ψ2. Piecewise models were fit with initial starting parameters of 1975 and 2000 for ψ1 and ψ2, respectively. The ψ1 and ψ2 parameters were tuned manually for 29 countries with a shortened data range, following visual inspection of data (see Tables S1 and S2).
Based on this piecewise regression model procedure, we then used parameters from Eq. (2) to determine three key indicator points of crop diversity change through time for each country (displayed visually in Fig. 1). Indicator 1 reflects the onset of diversification in each country, and was calculated as Breakpoint 1 (ψ1) in Eq. (2); this indicator therefore corresponds to the year in which notable changes in species richness began. Indicator 2 reflects the duration of the crop diversification period in each country, and was calculated as the difference between breakpoints 2 and 1 (i.e., ψ2-ψ1 from Eq. 2); this indicator therefore represents the duration of the period when crop prominent changes in crop diversity occurred. Finally, Indicator 3 reflects the rate at which crop diversity changed throughout the diversification period in each country; this indicator was calculated as the rate of crop diversity change (between ψ1 and ψ2), which in our models corresponded to the sum of the slopes (1) prior to the first breakpoint, and (2) between the first and second breakpoints (i.e., corresponding to b + c in Eq. 2). For each indicator we then calculated summary statistics as either mean ± standard deviations or median ± median absolute deviations (m.a.d.), where data was normally or log-normally distributed, respectively. Country values for each indicator were mapped using the mapCountryData function in the rworldmap R package40.
Changes in crop evenness over time
Evaluations of temporal changes in crop evenness at national scales followed this same analytical approach as above. First, for each country-by-year combination we calculated Pielou’s evenness index (J′)—which ranges from 0 to 1, with values closer to 0 indicating less evenness or greater abundance of a few dominant crop groups, and values closer to 1 representing more equitable abundances of crop groups—as:
$$J^{prime} = frac{H^prime }{{ln left( S right)}}$$
(3)
where S is again crop richness, and H′ is the Shannon–Weiner diversity index calculated as:
$$H^prime = – mathop sum limits_{i = 1}^{S} p_{i} ln p_{i}$$
(4)
where pi represents the relative proportion of the ith crop group for a given country-by-year combination. In these evenness calculations, all values of pi were estimated as the relative proportion of agricultural area (measured in ha) occupied by a given crop commodity group, within a country at a given year; this analytical approach was employed by Martin et al.8 when assessing crop group composition at supra-national scales. We then evaluated how J′ values changed in each country through time by replicating our stepwise modelling analyses above, substituting J′ for S in Eqs. (1) and (2), and extracting the same model indicators (Fig. 1). Finally, we calculated summary statistics and mapped each of these indicators, as described above.
Changes in crop composition across countries and over time
We used multivariate analyses to evaluate how temporal changes in S and J′ influenced crop composition across countries and over time. To do so, we created a community composition matrix whereby national-level crop assemblages were estimated for each of the country-by-year combinations. In this matrix, area harvested was taken as an approximation of the abundance of each crop group within each country-by-year combination (again following Martin et al.8). Since these abundances (or area harvested) across country-by-year combinations varied over orders of magnitude, we used non-metric multidimensional scaling (NMDS) to analyze and visualize spatial (country) and temporal (year) differences in crop diversity. Specifically, we used the vegan R package38 to calculate all 58,899,231 Bray–Curtis dissimilarities among all 10,854 data points (i.e., crop group composition in every country-by-year data point), as:
$$BC_{jk} = frac{{sum i left| {x_{ij} – x_{ik} } right|}}{{sum i left( {x_{ij} + x_{ik} } right)}}$$
(5)
where BCjk represents the dissimilarity between the jth and kth community, xij represents the abundance (i.e., area harvested) of crop group i in sample j, and xik represents the abundance of crop group i in sample k. We then used a multivariate analysis of variance (i.e., an Adonis test), to test for significant differences in Bray–Curtis distances as a function of country, year, and a country-by-year interaction. Significance was assessed using a permutation test, with 99 permutations used.
Latitudinal gradients in crop richness
To test our hypotheses surrounding the presence of, and temporal changes in, latitudinal gradients in crop group diversity, we focused on 164 countries for which crop group diversity was available in both 1961 and 2017. For each of these two datasets, we fit a separate linear regression model that predicts crop group richness as a function of latitude (expressed as an absolute value) and a 2nd-order polynomial term for the ‘latitude2’ variable. From both of these models, we extracted and compared latitude value at which crop group richness was estimated/ modelled to peak.
Predictors of change in crop diversity and composition
We tested if Human Development Index (HDI) was correlated with patterns of change in crop diversity and composition. Briefly, the HDI is a composite index of four metrics related to socio-economic status, including life expectancy at birth, expected years of schooling for children at a school-centring age, mean years of schooling for adults ≥ 25 years of age, and log-transformed gross national income per capita. These values are then aggregated on a per country basis, into an HDI index that ranges from 0–1 with higher scores denoting higher performance in these indicators. We employed 2017 HDI values in our analysis here, in order to include the most countries possible in each analysis (since earlier HDI scores are less readily available)41.
We then used linear mixed effects models to test if patterns of change in crop diversity and evenness varied systematically with HDI values. This entailed fitting six linear mixed models, where each of our six indicators (i.e., Indicators 1–3 for both S and J′) were predicted as a function of HDI; these models also accounted for potential spatial autocorrelation in Indicator values by including the FAO-defined continent identity and FAO-defined region identity of each country, as a nested random variable. Models were fit using the lme function in the nlme R package41. We then estimated the proportion of variation in each indicator that is explained by HDI, continent identity, and region identity, using the varcomp function in the ape R package42—which partitioned explained variation across continents and regions—as well as the sem.model.fits function in the piecewiseSEM R package43—which partitioned explained variation across the fixed (i.e., model intercept and HDI) vs. random (i.e., continent and region) effects. Due to differences in HDI data availability and in the number of piecewise models that converged, n = 152 countries for all models of S indicators and n = 139 countries for all models of J′ indicators. Log-transformed values of Indicators were used in these analyses where they better approximated a log-normal distribution, as determined using the fitdistrplus function in the fitdistrplus R package44.
Source: Ecology - nature.com
