Grain and soil samplings from farmers’ fields were completed in November–December 2017 (for most of Amhara region) and November 2018–January 2019 (Amhara, Oromia and Tigray regions) in Ethiopia, and in April–June 2018 in Malawi. The work involved sampling of grain and soil from farmers’ fields and grain stores with the informed consent of the farmers. The work was conducted under ethical approvals from the University of Nottingham, School of Sociology and Social Policy Research Ethics Committee (REC); BIO-1718-0004 and BIO-1819-001 for Ethiopia and Malawi, respectively. These REC approvals were recognized formally by the Directors of Research at Addis Ababa University (Ethiopia) and Lilongwe University of Agriculture and Natural Resources (Malawi), who also reviewed the study protocols.
The objective of this study was to support the spatial mapping of grain Ca, Fe, Se and Zn concentration of cereal crops. We sought reasonable spatial coverage over the target sample frame and used ‘main-site’ and ‘close-pair’ sampling to support the estimation of parameters of the spatial linear mixed model (LMM)22,43.
In Ethiopia (Amhara, Oromia and Tigray regions), target sample frames were constrained to locations at which the probability of the land being under crop production was ≥0.9 based on predictions produced on a 250-m grid. These predictions were derived from the interpretation of high-resolution satellite imagery by trained observers and by machine learning methods applied to multiple covariates derived from remote sensing data and digital elevation models44,45. The sample frame was further constrained to include only those locations from a 250-m grid that fell within 2.5 km of a known road. A map indicating nodes on a 250-m grid (with the same origin as the agricultural land-use grid) that met this requirement was prepared. These constraints may introduce possible biases into predictions made at locations outside the designed sample frame, however, it would not otherwise have been possible to visit all of the sample locations in the time available. Information on the distribution of roads was taken from OpenStreetMap (www.openstreetmap.org). Of a total land area of around 1.1 million km2 in Ethiopia, the total cropland mask represented 354,325 km2, of which 220,467 km2 was within 2.5 km of a known road. In Malawi, the cropland area was determined from the European Space Agency Climate Change Initiative46. The agricultural area used was defined as including all raster cells that included the category of ‘cropland’ in their description. In Malawi, where road access to cropped areas is generally better than in Ethiopia, no constraint to road distance was imposed on sample locations. The mapped cropland areas are shown in Fig. 1.
In Ethiopia, a total of 1,825 primary sample locations were selected a priori, with each 250-m grid node within the sampling frame allocated an equal prior inclusion probability. This was done using the lcube function from the BalancedSampling library for the R platform47,48. The lcube function implements the cube method, to enable random sampling according to specified inclusion probabilities while aiming for balance and spread with respect to the specified covariates49. Here, inclusion probabilities were uniform across the sample frame and sample locations were selected for spatial balance, which entails that the mean coordinates of sample locations are close to the mean coordinates of all points in the sample frame, and spatial spread, which ensures that the observations are spread out rather than clustered with respect to the spatial coordinates50. A subset of 175 of these locations were selected as close-pair sites at which an additional nearby sample was taken to support the estimation of parameters of the spatial LMM43.
In Malawi, a different sampling method was used to achieve good spatial coverage of a total of 1,710 main-site locations. These included 820 fixed sample points from the 2015/16 Demographic and Health Survey of Malawi24,51. The stratify function in the spcosa library for the R platform52 divides a sampling domain into Delaunay polygons centred on a set of fixed points and with the remaining polygon centroids selected to partition the domain into approximately equal-area regions. The centroids of the polygons were selected as sample points. An additional 890 sample points were found in addition to the 820 fixed ones, with the stratify function. Once these were obtained, a further 190 locations were selected at random as close-pair sites for an additional nearby sample, as in Ethiopia.
Sampling was conducted by teams who were trained in standard procedures and risk assessments. Each team planned to visit five main-site locations per day. Main-site locations were loaded onto a computer tablet and printed on paper maps for each team. A team would navigate to the target location, using a handheld GPS device for the last few kilometres. At each sample location, the team would identify the nearest field with a mature cereal crop within a 1-km radius, and sample grain and soil, subject to farmer consent. If a field with a standing mature cereal crop was not apparent, that is, the crop had been harvested, or a non-cereal crop had been grown, the team would ask the farmer to identify a field from which a cereal crop had recently been harvested and stored, and from which a sample could be obtained. If sampling was not possible, then the team would either look beyond a 1-km radius for an alternative location, or abandon the location. At designated close-pair locations, a second field was identified ideally within around 500 m (range, 100–1,000 m) of the main-site location. If a close-pair location could not be found, then a close-pair location would be selected at the next sample location that was not already earmarked for a close-pair sample.
Within a selected field, samples were taken from a 100 m2 (0.01 ha) circular plot. This was centred as close as practical to the middle of the field unless this area was unrepresentative due to disease or crop damage. Five subsample points were located (Extended Data Fig. 1). The first point was at the centre of the plot. Two subsample points were then selected at locations on a line through the plot centre along the crop rows, and two more points on a line orthogonal to the first through the plot centre. Where possible, the central sampling location was fixed between crop rows, and the long axis of the sample array (with sample locations at 5.64 and 4.89 m) was oriented in the direction of crop rows with the short axis perpendicular to the crop rows. A single soil subsample was collected at each of the five subsample points with a Dutch auger with a flight of length of 150 mm and diameter of 50 mm. The auger was inserted vertically to the depth of one flight and the five subsamples were stored in a single bag. Where a mature or ripe crop was still standing in the field, grain samples were taken close to each augering position by a different operator, to minimize further contamination by dust and soil. For maize, a single cob was taken at each of the five points. Maize kernels were stripped from around 50% of each cob lengthways and composited into a single sample envelope for each location. For smaller-grained crops, sufficient stalks were taken so that approximately 20–50% of the sample envelope was filled (dimensions 15 cm × 22 cm), with samples placed grain‐first into the sample bag and the stalks were twisted off the grain heads and discarded. If a crop was in field stacks, then a subsample, comprising five cobs for maize, or a representative sample for other crops was taken from each available stack, taking material from inside the stack to minimize contamination by dust and soil (Extended Data Fig. 1). If a crop was in a farmers’ store, that is, already averaged from across the field, then a representative sample was taken, while avoiding grain from the store floor if grain was loosely stored and avoiding grain with visible soil or dust contamination.
Photographs at sample locations and of sample bags were recorded for quality assurance along with site GPS locations. In Ethiopia, 1,385 of the 1,389 locations from where grain data are reported had positional uncertainties of ≤8 m as recorded by the GPS. The other four locations had positional uncertainties of 9–16 m. In Malawi, 1,790 of the 1,812 locations had positional uncertainties of ≤9 m. A further 16 locations had positional uncertainties of 10–17 m, and 6 locations had positional uncertainties of 2,900–5,000 m. We took a decision not to exclude any data based on positional uncertainties for this study. We used robust estimators of the variograms (Extended Data Fig. 2), which are resistant to effects of spatial outliers due to a small number of points being in the wrong position53 and these models were validated. Any effect of position error on the broad mapped pattern of long-range spatial variation at these scales will therefore be very limited and localized.
Whole-grain samples were air-dried in their sample bags. Each sample was then ground in a domestic stainless-steel coffee grinder, which was wiped clean before use and after each sample with a non-abrasive cloth. Whole grains are consumed in many settings, although more-refined fractions—with concomitant losses of more micronutrient-rich bran and endosperm fractions—are often consumed. All preparation was done away from sources of contamination by soil or by dust. A 20-g subsample of the ground material was then shipped to the University of Nottingham. Soil samples were oven-dried at 40 °C for 24–48 h depending on the moisture content of the soil. Preparation took place in a soil laboratory to avoid cross-contamination with grain samples. Plant material was removed from each soil sample, which was then disaggregated and sieved to pass 2 mm. This material was then coned and quartered to produce subsample splits. A 150-g subsample of soil was poured into a self-seal bag, labelled and shipped to the UK for analysis in the laboratories at the University of Nottingham and Rothamsted Research. Soil pH(water) and soil organic carbon (SOC) content were measured using standard methods22.
Grain micronutrient analyses
Grain micronutrient analysis methods followed standard methods54. Approximately 0.2 g of each ground sample was weighed and digested using a microwave system. For samples collected in the Amhara region in 2017, a Multiwave 3000 48-vessel MF50 rotor (Anton Paar) was used; digestion vessels were perfluoroalkoxy tubes in polyethylethylketone pressure jackets (Anton Paar). Samples were digested in 2 ml 70% trace-analysis-grade HNO3, 1 ml Milli-Q water (18.2 MΩ cm; Fisher Scientific) and 1 ml H2O2. Settings were: 1,400 W, 140 °C, 2 MPa, for 45 min. For samples collected in 2018–2019, we used a Multiwave Prom 41HVT56 rotor and pressure-activated venting vessels made of modified polytetrafluoroethylene (56-ml ‘SMART VENT’, Anton Paar). Samples were digested in 6 ml of 70% trace-analysis-grade HNO3. Settings were: 1,500 W, 10 min heating to 140 °C, 20 min holding at 140 °C, and 15 min cooling to 55 °C. Two operational blanks were typically included in each digestion run. Duplicate samples of a certified reference material (Wheat flour SRM 1567b, NIST) were included in approximately every fourth digestion run. Following digestion, each tube was made up to a final volume of 15 ml by adding 11 ml Milli-Q water, then transferred to a 25-ml universal tube (Sarstedt) and stored at room temperature. Samples were further diluted 1:5 with Milli-Q water into 13-ml tubes (Sarstedt) before analysis.
Multi-elemental analysis of grain (Ag, Al, As, B, Ba, Be, Ca, Cd, Cr, Co, Cs, Cu, Fe, K, Li, Mg, Mn, Mo, Na, Ni, P, Pb, Rb, S, Sr, Ti, Tl, U, V and Zn) was undertaken using inductively coupled plasma mass spectrometry (ICP–MS; iCAPQ, Thermo Fisher Scientific). The instrument used a helium collision cell with kinetic energy discrimination to reduce polyatomic interferences. Samples were introduced from an autosampler incorporating an ASXpress rapid uptake module (Cetac ASX-520, Teledyne Technologies) through a perfluoroalkoxy Microflow PFA-ST nebuliser (Thermo Fisher Scientific). Internal standards were introduced to the sample stream on a separate line via the ASXpress unit and included Sc (20 μg l–1), Rh (10 μg l–1), Ge (10 μg l–1) and Ir (5 μg l–1) in 2% HNO3 (Primar Plus grade; Fisher Scientific). An external multi-element calibration standard (Claritas-PPT grade CLMS-2; SPEX Certiprep) was used to calibrate Ag, Al, As, B, Ba, Be, Cd, Ca, Co, Cr, Cs, Cu, Fe, K, Li, Mg, Mn, Mo, Na, Ni, P, Pb, Rb, S, Se, Sr, Tl, U, V and Zn, in the range of 0–100 μg l–1 (0, 20, 40, 100 μg l–1). A bespoke external multi-element calibration solution (PlasmaCAL, SCP Science) was used to create Ca, K, Mg and Na standards in the range of 0–30 mg l–1. B, P and S calibration used in-house standard solutions (KH2PO4, K2SO4 and H3BO3); Ti was determined semiquantitatively. Sample processing was undertaken using Qtegra software (Thermo Fisher Scientific) with external cross-calibration between pulse-counting and analogue detector modes when required. Se was determined separately using a triple quadrupole ICP–MS (iCAP TQ; Thermo Fisher Scientific) using an oxygen cell to mass shift the isotope 80Se to m/z 96 (80Se16O) to reduce interference from the 40Ar dimer. Drift correction was achieved using Rh as an internal standard; calibration used the CLMS-2 multi-element standard (Certiprep).
Analyses were conducted in batches of around 240 samples per run on the ICP–MS instrument (Extended Data Table 1). Individual grain concentration data were corrected for run-specific operational blanks and then converted to concentration on a dry-matter basis. Element-specific limits of detection were reported as 3× the standard deviation of the operational blank concentrations, assuming a notional starting dry weight of 0.2 g of sample material (Extended Data Table 1). For samples for which the grain element concentration was less than the limit of detection, data were removed before the statistical analyses. No adjustment was made for potential contamination of grain samples, for example, with soil dust from the field or store using typical markers (for example, Fe, V). Two sorghum samples, taken from a grain store in Malawi, were excluded from the data analysis based on high concentrations of Ca, Mg and other elements that were considered unlikely to have arisen from soil contamination.
Summary statistics were computed for concentrations of all elements in grain and histograms were examined. It is common for geochemical variables to be log-normally distributed, and the coefficient of skewness of the data was examined using octile skewness as a robust measure of asymmetry of the distribution55. The data were analysed on the original scales of measurement (mg kg–1) if the octile skewness was within the interval [−0.2, 0, 2] as described previously56. If the octile skewness fell outside this range (with a positive value), and the absolute value of the octile skewness after loge transformation was smaller than on the original scale, then the data were analysed on the loge scale (Extended Data Tables 1, 2).
Variograms were estimated for each variable using three estimators (Extended Data Table 2): the standard estimator57 and two alternatives58,59. The two alternative estimators are more robust than the standard estimator to outlying data, be these marginal outliers (apparent on the histogram) or spatial outliers (apparent in local spatial context). The variogram estimates were formed from all pairwise comparisons among observations up to a maximum lag distance (330 km in Ethiopia and 100 km in Malawi), and with lag bins with a width of 10 km. The distance between any two locations (specified by latitude and longitude) was computed as the great circle distance using the distVincentySphere function from the geosphere library of the R platform47,60. Exponential variogram functions were then fitted to the estimates by weighted least squares61. The exponential variogram was selected because it ensures positive definite covariance matrices for distances on the sphere62. Each model was then tested by cross-validation. The selected models are shown in Extended Data Fig. 2.
In cross-validation, each observation was removed and predicted from the remaining observations by ordinary kriging. This was done using each variogram model fitted for each variable to the sets of point estimates computed by the three estimators. The standardized squared prediction error (SSPE) was computed for each observation as the squared difference between the observed and predicted values divided by the prediction error variance (kriging variance) as computed in the standard ordinary kriging equations53. The median value of the SSPE has an expected value of 0.455 in the case of a valid underlying variogram model with normally distributed kriging errors53. The standard estimator due to Matheron57 is more statistically efficient than the robust alternatives, so if the model fitted to these estimates appeared correct from the cross-validation results (the median SSPE is within the 95% confidence interval) then the alternatives were not considered. If the SSPE suggests that the model fitted to Matheron estimates are affected by outliers, then the models fitted to robust estimates were also cross-validated, and one was selected on the cross-validation results53.
Once a variogram model was selected it was used to compute predictions of the grain concentration (with or without transformation) on nodes of a fine square grid. This was done by ordinary kriging61. In ordinary kriging, it is assumed that the mean value of the variable of interest is locally constant, but unknown. An estimate is found that is a linear combination of the observations such that the expected squared error of the prediction (the ordinary kriging variance) is minimized. The kriging variance is also computed as a measure of the uncertainty of the prediction. For a given variogram, the kriging variance reflects the proximity of the location at which a prediction is made to observations of the variable of interest. In the case of variables that had been transformed to logarithms, the prediction at each location was computed by exponentiation of the prediction on the loge scale. The resulting back-transformed estimate is a median-unbiased predictor, which is appropriate for variables with a skewed distribution63. Note, the kriging variance cannot be meaningfully back-transformed. However, it still indicates how the uncertainty of the predictions varies in space and so it is mapped here on the loge scale for transformed variables. The kriging variance maps are shown in Extended Data Figs. 3, 4. In Ethiopia, where the sample sites are irregularly distributed over the mapped area, the kriging variance differs markedly depending on the local sample density.
In Ethiopia, we mapped the percentage of dietary requirement potentially met by the intake of wheat and teff (Fig. 3); similarly, in Malawi for maize (Fig. 4). In Ethiopia, 97.6 g per capita per day of wheat and 89.3 g per capita per day for teff (based on the Food and Agriculture Organization database item, ‘other cereals’) were used as reference intakes29; in Malawi, 342.8 g per capita per day was used as a reference maize intake29. Estimated average requirement (EAR) thresholds of 860, 22.4, 45 and 10.2 mg per capita per day for Ca, Fe, Se and Zn, respectively, were chosen as a representative threshold, based on an adult woman aged 18–24 years eating an unrefined (that is, high phytate) diet30. These thresholds are similar to other demographic groups. Because a comparable measure of uncertainty to the kriging variance for this derived variable is not available, we defined a mask for Ethiopia, where the sparsity of sampling leads to larger uncertainty compared with Malawi. We considered the Zn concentration in teff, a variable with spatial dependence over long distances, and identified those areas in which the kriging variance for this variable exceeded 75% of the variance of the variable itself because the sampling was sparse. These areas defined the mask used for the maps of the percentage of dietary requirement for all variables and is shown in grey in Fig. 3.
The maps shown in Figs. 3, 4 are made by point kriging—that is, they are predictions of a measurement at an unsampled site. Ordinary kriging can also be used to predict the mean value of a variable across a region or ‘block’61. To illustrate the communication value of this approach, the mean concentrations of Zn and of Ca in grain (teff and wheat) are shown at the woreda level in the Amhara region (Fig. 5). These were obtained by block kriging of woreda means from the same sample data, and with the same variogram models as used to produce the point kriging predictions.
Determining the links between a dietary biomarker of status and grain micronutrient concentration focused here on Se, for which there is a reliable biomarker of dietary status13; this is not the case for Zn7 and the other micronutrients analysed in this study. Data on the concentration of Se in the blood of women of reproductive age were available from micronutrient surveys in Ethiopia (serum23) and Malawi (plasma24). Mean concentrations were computed for each enumeration area available—321 in Ethiopia and 101 in Malawi—for which the latitude and longitude for each enumeration area centroid were available. For each enumeration area, the nearest grain sample site (regardless of crop) was found. In Ethiopia, the median distance to the nearest grain sample site over all enumeration areas was 17 km. In Malawi, the median distance was much shorter (1.4 km), which is attributable both to the denser crop sampling in Malawi and to the fact that enumeration areas used for the micronutrient survey were targeted for sampling.
The enumeration areas do not comprise a simple independent random sample, so the relation between the serum or plasma Se concentration and the Se concentration in the nearest grain sample could not be quantified by a statistic such as the correlation coefficient. It was therefore studied with an appropriate LMM incorporating a spatially correlated random effect, modelled with a Matérn correlation function64. Exploratory analysis indicated that a loge transformation of all serum or plasma Se concentrations was necessary for the assumption of normal random effects to be plausible. The observed serum or plasma Se concentration was modelled as a linear function of the concentration in grain at the nearest sample location. Residual maximum likelihood was used to estimate the random effects parameters. The fixed effects were then estimated by weighted least squares65 along with their standard errors. The evidence against the null hypothesis that the regression coefficient for serum or plasma Se on grain Se concentrations was zero was tested by a log-likelihood ratio test.
Plots of the serum or plasma Se concentration in women of reproductive age, in Ethiopia and Malawi, respectively, show a positive correlation between the variables (Extended Data Fig. 5), which is supported by the statistical models. For Ethiopia, the estimated regression coefficient (log[ng serum Se per ml] and log[mg of grain Se per kg]) was 0.08 with a standard error of 0.02. The null hypothesis that the coefficient was zero is rejected on the grounds of the log-likelihood ratio statistic (L = 14.48, P = 1.4 × 10−4). If there was no relationship between the biomarker and the grain Se concentrations, the probability of obtaining an L-statistic this large or larger would be very small. Similarly, for Malawi, the estimated regression coefficient was 0.09 with a standard error of 0.03 (L = 11.56, P = 6.7 × 10−4).
Data on the concentration of Se and of Zn in grain (maize in Malawi; teff, wheat and maize in Ethiopia) were extracted along with the corresponding data on soil pH(water) and SOC. Observations for three environmental covariates were also extracted for these same locations: (1) mean annual temperature; (2) mean annual precipitation values from the CHELSA datasets66,67, which are downscaled to a spatial resolution of 30 s; (3) topographic index values from the 30-s resolution MERIT Digital Elevation Model68. The topographic index—which is sometimes called the topographic wetness index—is a measure of the tendency for drainage to accumulate at a site and is widely used as a predictive measure for soil properties. Following exploratory analysis, grain Se concentration data were loge-transformed to make the assumption of normal random effects plausible; this was not necessary for grain Zn. Measurements of SOC from Malawi showed a marked positive skew and were therefore loge-transformed.
Data were analysed using a LMM in which a regression-type function of environmental covariates was considered as a fixed effect along with a spatially autocorrelated random effect, as described for the biomarker data. Random-effect parameters were estimated by maximum likelihood, and then the fixed-effect parameters by weighted least squares. The first model was fitted with a constant mean as the only fixed effect. The second model was fitted with mean annual precipitation as a fixed effect. The evidence for this predictor was evaluated by a log-likelihood ratio test of the second model against the first. The predictor was retained initially if the probability of obtaining a value of L as large or larger than observed under the null hypothesis of no effect of the predictor was <0.05. Additional predictors were then considered in the order of mean annual temperature then topographic index. Finally, the P values were evaluated as the outcome of multiple tests, controlling the false-discovery rate (FDR)69 at 0.05. The FDR is the expected proportion of incorrectly rejected null hypotheses among all rejected ones. For each micronutrient in each crop and country, the predictors were identified that could be regarded as significant with FDR control at 0.05. The same procedure was followed to produce a comparable model based on the two soil properties, considering first soil pH and then SOC.
Plots of grain Se and Zn concentrations against the environmental covariates and soil properties are shown in Extended Data Figs. 6–9. The evidence for environmental covariates or soil properties as predictors of micronutrient content in the grain is summarized in Extended Data Table 3. Extended Data Table 3 also shows the estimated coefficients, and their standard errors, in separate models for each micronutrient; for maize, teff or wheat in Ethiopia and for maize in Malawi. The predictors in these models are only those that were selected with FDR control.
Further information on research design is available in the Nature Research Reporting Summary linked to this paper.