Pollinators support the nutrition and income of vulnerable communities

Study sites

Field work took place in ten smallholder farming villages (2,400–3,000 metres above sea level, temperate climate) in Jumla District, Nepal (Fig. 1, Supplementary Fig. 2 and Supplementary Method 2). Jumla is a remote mountainous district, situated in the Karnali Province of western Nepal. Rates of poverty, food insecurity and malnutrition are very high in this region, and 80% of the population directly depends on smallholder agriculture⁴⁷. The study populations are typical of many smallholder communities around the world (Supplementary Fig. 1) and are characterized by a heavy reliance on small (<2 ha) family-run farms for both their income and nutrition, which makes them heavily dependent on local ecosystem services²⁷. Each study village comprised a cluster of 100–400 closely spaced households interspersed with small vegetable gardens and livestock enclosures. Villages were surrounded by many small (0.01–0.3 ha) arable fields and apple orchards as well as large areas of steep, heavily grazed pasture and native forest (Supplementary Figs. 3 and 4). More than 50 crops are grown in this region, including many pollinator-dependent species such as apples, beans, pumpkins, mustard and buckwheat (Supplementary Tables 2–4).

Study population

Households were considered eligible for inclusion in the study if they were permanent residents in the community (spent at least 10 months of the year in the village) and had at least three out of the four respondent categories in the household (adult man, adult woman, adolescent girl and child under 5 years of age). On the basis of a full census of each village, we identified eligible households, randomly selected 20 from each village and obtained consent from the head of the household. From each household, we aimed to enrol one adult woman, one adult man, one adolescent girl and one child under the age of 5 years as participants in the study. This selection provided a diverse picture of diets in each household, including two particularly vulnerable subgroups: young children (owing to their rapid growth and high nutrient needs) and adolescent girls (who may soon have their first child and whose pre-conception nutrition strongly influences maternal and child health outcomes). Adolescent boys were not included as a separate respondent category as their diets are assumed to be more similar to adult men and they are considered less nutritionally vulnerable than adolescent girls. Participants were considered eligible for inclusion in the study if they were permanent residents of the household and not suffering with an illness or disability that prevented them eating a usual diet or responding to questions (Supplementary Method 3, selection of study participants). Our final study population (Supplementary Table 11) consisted of 776 individuals (215 adult women, 186 adult men, 190 adolescent girls and 185 children under 5 years of age).

Ethics statement

Ethical approval for this study was obtained from the Ethical Review Board of the Nepal Health Research Council (reference 1709) and the Faculties of Life Sciences and Science Research Ethics Committee at the University of Bristol (reference 102982). All procedures involving human participants were conducted in accordance with the relevant institutional and national ethics guidelines. For all participants over 18 years of age, informed consent (signature or thumb print) was obtained; for all participants under 18 years, consent was provided by their parent or guardian and the adolescent girls also provided assent. Participants’ time was remunerated with mobile phone credit vouchers, a widely used and valued resource in the study region.

Dietary-recall surveys and recipe collection

Dietary-recall surveys were conducted for each participant every 2 weeks for an entire year (November 2021 to November 2022). During each recall event (24 in total), we recorded the identity and quantity of every food item consumed by respondents during the previous 24 h. To minimize recall bias, we used a five-stage multipass food-recall method⁴⁸ (Supplementary Method 3, dietary-recall surveys). Food models were used to estimate portion sizes: cooked rice as a model for rice and other irregular-sized solid foods such as vegetable curries and cooked green leaves; water for portions of dal (soup of beans or pulses), milk or other liquids; playdough for chunks of meat or pieces of fruit; wheat flour for powdered foods like roasted grain flour, salt and chili powder; and dried corn for dried or roasted nuts and grains. Participants were instructed to only report the food that was consumed and not any that was leftover on the plate. Food models were weighed using Salter kitchen scales and then back-converted to estimate the mass of each food item consumed, accounting for differences in food density (Supplementary Method 3, dietary-recall surveys). For foods in packets, or of a standard size, the number of items consumed was recorded. Interviews were conducted by trained data collectors in the Nepali language and information was entered into a customized data-collection form using the cloud-based data-collection platform CommCare (v.2.49; http://www.commcarehq.org/home/) on an Android tablet. All data collection had range checks and internal validity checks built in to help maintain quality control.

To identify the composition of each mixed dish (that is, the quantity of each ingredient in it), we asked the lead cook in replicate local households to prepare the food following their normal recipe (Supplementary Method 4). Each ingredient was weighed as it was added to the cooking dish and then the final weight of the prepared meal was recorded and the proportional composition of each ingredient (by weight) was calculated. For each ingredient in a recipe, we calculated its mean proportional contribution to the recipe (food item) based on the raw weight of edible ingredients across all recipe replicates collected (ten in most cases). For each food item consumed by a participant, we multiplied the grams of food item consumed by the proportional content of each ingredient to determine the grams of each individual ingredient consumed on a given date. We complied a Nepal-specific food composition table from a range of published sources (Supplementary Method 4, food composition data sources), which was used to assign a nutritional composition to each ingredient. Nutrient values for each ingredient were adjusted for cooking losses using the USDA Table of Nutrient Retention Factors (release 6, 2007), applying the appropriate factor for each cooking method (for example, boiled versus fried potato) to ensure that nutrient values reflected the ingredients in the form they were eaten (Supplementary Method 4, nutrient retention factors).

Height and weight measurements

For each study participant, measurements of their height and weight were taken once per month throughout the survey period to assess their nutritional status. Standing height was recorded using Shorr Board Stadiometers for adults, adolescents and children over 2 years old (or >87 cm in height), whereas length (that is, lying down) was recorded for children under 2 years old or <87 cm. Weight was measured using Seca 877 scales. All measurements were taken in duplicate, with a third reading taken if the difference between the first two readings exceeded a set cutoff. The midpoint of duplicate measures was used. Further details are provided in Supplementary Method 5.

Probability of adequacy calculations

To assess dietary adequacy, we calculated probability of adequacy (PA) values for 11 micronutrients of public health concern (calcium, zinc, vitamin C, thiamin (B1), riboflavin (B2), niacin (B3), vitamin B6, folate (B9), vitamin A, vitamin B12 and iron). Usual intakes were estimated as the mean of all available 24-h recalls for each individual. Participants were surveyed fortnightly over a full year, which provided up to 24 recalls for each participant. With this number of repeated measures, the influence of day-to-day variation on mean intake estimates is likely to be small, and the mean was therefore considered an acceptable approximation of usual intake⁴⁹. PA values were derived using the standard probability approach⁵⁰, with estimated average requirement values applied according to age, sex and physiological status (pregnancy, lactation, pre-menarche and post-menarche). Iron was treated separately owing to its non-normal requirement distribution, and low bioavailability assumptions (low absorption by the body) were applied for iron (5%) and zinc (15%) following previously described methods⁵¹. Further details are provided in Supplementary Method 6.

Plant–pollinator visitation surveys

To characterize the plant–pollinator interaction network that underpins the delivery of the pollination service, we conducted plant–pollinator visitation surveys every 2 weeks from 18 April to 4 November 2021 (spring to autumn). These surveys immediately preceded the dietary-recall surveys but were not synchronous owing to COVID-19 restrictions, which prevented in-person interviews at the time of the ecological surveys. However, most of the crop produce that was grown and pollinated during the plant–pollinator surveys was consumed throughout the period of dietary-recall surveys. Moreover, agricultural production and diets remain similar from year to year in this region (Supplementary Table 12).

Plant–pollinator surveys were conducted in a 600 × 600 m (36 ha) sampling area centred on the midpoint of each study village. This area was divided into three habitat categories: village, crop and semi-natural vegetation. In each of these habitats, we randomly located three replicate fixed survey plots of 60 × 60 m (9 plots per village). Every 2 weeks, a timed 40-min survey was conducted for each plot to record the interactions between plants (both crop and non-crop species) and flower-visiting insects (those observed collecting resources from the flower). Surveys were conducted by a team of ten trained data collectors (recruited from within each village) using the cloud-based data-collection platform CommCare (v.2.49; http://www.commcarehq.org/home/). Plants were identified by data collectors using a custom-made Plant Atlas for Jumla District (with scannable QR codes to minimize transcription error), and species identifications were checked for accuracy by a botanist from the national herbarium of Nepal. All insect specimens were identified by insect taxonomists (see Acknowledgements for details) and are stored in the Central Department of Zoology, Tribhuvan University, Nepal. Insect identifications were checked for accuracy by an independent taxonomist from the University of Agricultural Sciences, Bangalore, India. Further details are provided in Supplementary Method 7. Sampling permission was provided by the Ministry of Forest and Environment, Nepal (reference 258).

Estimating pollinator importance

Without detailed studies of the behaviour and morphology of each insect taxon, and on the seed set of plants exposed to different pollinators, we cannot know which flower visitors are truly pollinators⁵². Here we assumed that any insect that visits a flower is a potential pollinator and we estimated pollinator importance on the basis of its visitation frequency, which has been shown by both empirical³⁹ and theoretical⁵³ work to provide a good indication of ecosystem service delivery. However, visitation information alone can overstate the importance of abundant but less effective pollinator taxa such as butterflies⁵⁴. Thus, to estimate the pollination service delivery of a pollinator taxon, we used a modified version of visitation frequency and weighted each visitation event by the pollen transport capacity of the insect (mean number of pollen grains from all plant species carried on its body, as per a previously described method⁴¹). This metric is based on the assumption that insects that frequently visit a focal crop and carry large amounts of pollen on their bodies are likely to be better pollinators than insects that visit less frequently and carry little pollen, consistent with previous studies^55,56. This estimate of pollination service provision, although imperfect, is a necessary simplification for a large-scale study involving multiple plant and pollinator species across multiple sites. It will help identify dominant crop pollinators, even if the values are more relative than absolute. Pollen transport capacity was estimated by swabbing individuals of each insect taxon and calculating the mean number of pollen grains from all plant species carried on their body and available for pollination. Pollen grains were not identified to species (which is impractical for a study of this scale), so our measure of pollen-carrying capacity is limited to one value per insect taxon rather than specific values for each plant–pollinator combination (more details provided in Supplementary Method 8). For each insect, our pollinator importance (PI) metric estimates its proportional contribution to the total pollen transport for a given crop species based on the following equation:

where V_i,c is the visitation frequency of insect i to crop c and P_i is the pollen transport capacity of insect i. Thus, insects that carry a lot of pollen and make up a high proportion of all visits to a crop will receive a high importance score.

Plant–pollinator interactions that underpin human nutrition

By combining participants’ dietary data with published nutrient composition values for each ingredient, we were able to determine the origin of each nutrient in participants’ diets (that is, the plants and animals it was sourced from). We then calculated the proportion of each nutrient obtained from local versus imported ingredients, as well as pollinator-dependent versus pollinator-independent ingredients (see the categorization process below). Once the source of each nutrient was established, we identified the insects that supported its production via their pollination of the nutrient-provisioning crops.

Categorizing foods and calculating their contribution to diets

Each individual food item or ingredient was classified as either locally derived (that is, grown or raised within the isolated district of Jumla) or imported (from outside Jumla) based on responses to the dietary-recall survey, which included a question on food provenance. The isolation of the study villages from regional markets made it relatively straightforward to determine which ingredients were imported; these were typically non-perishable packaged items such as polished white rice, vegetable oil and spices. Ingredients were also classified as either pollinator-dependent (production increased by animal pollination) or pollinator-independent (production not affected by animal pollinators) based on previous studies^11,46. In our study system, farmers are rarely able to access or afford external seed sources; therefore, all households sow their crops with local seeds saved from the previous year (Supplementary Table 10). Thus, we consider crops as pollinator-dependent if their seed production relies on animal pollination, even if the edible portion of the crop does not (for example, carrots and mustard greens).

Calculating the contribution of insect taxa to diets

To provide an aggregate estimate of the contribution of each insect to human nutrition, we multiplied its PI score for a given crop (see the section above) by the proportional contribution of the crop to the supply of a given nutrient. To account for differences in pollinator dependence among crops, the PI scores for each crop were multiplied by the percentage pollinator dependence of the crop defined as the reduction in yield recorded when flower visitors are excluded (values were taken from the literature and our own pollinator exclusion experiments; Supplementary Table 4 and Supplementary Method 9). Thus, for a given insect (i), its pollinator contribution (PC) to the supply of a given nutrient (n) is calculated as:

where PI_i,c is the importance score of pollinator i for crop c (see the equation for PI_i,c above), D_c is the pollinator dependence of crop c and N_c,n is the proportional contribution of crop c to the nutrient of interest (n). Thus, an insect that carries lots of pollen and frequently visits pollinator-dependent crops with a high nutritional value will receive a high score. The score represents the proportion of nutrient n that is directly attributable to pollination by insect i.

Pollination changes and how they affect human welfare

Understanding how declines in pollinator populations affect human health and livelihoods is challenging because direct long-term data are lacking. To address this issue, we used simulations to explore a range of plausible scenarios of pollinator population change and their consequences for diets and incomes. Although such simulations inevitably involve assumptions, their purpose is not to provide precise forecasts but to generate boundary estimates that illustrate the scale of potential impacts without any adaptation or mitigation responses. We combined crop–pollinator interaction data with dietary-recall data and household income data (recorded in questionnaires; Supplementary Method 10) to estimate the consequences of changes in the pollinator community. For each crop consumed by a participant, or sold by a household, we first calculated the contribution of individual pollinator taxa to its yield, based on the crop’s pollinator dependence and the proportion of pollen transport provided by that taxon (see above). This approach enabled us to estimate the portion of each participant’s nutrient intake and the proportion of each household’s farming income attributable to each pollinator species and to alter it according to the simulations described below.

Simulating declines in pollinators

We modelled declines in pollinator populations by simulating stepwise reductions in pollinator abundance from 0 to 100% loss in 10% increments. For each decline level, every pollinator taxon was assigned a species-specific decline rate drawn from a normal distribution, with the target mean decline and a standard deviation of 0.2, to reflect interspecific variation in sensitivity to environmental change. Rates were bounded between 0 and 100%. We then reduced each participant’s nutrient intake and each household’s farming income based on the share of their diet or income that was attributable to the pollinator in question. For each decline level, we repeated simulations 100 times to generate distributions of possible outcomes.

Although our analysis produced results across the full gradient of pollinator decline, in the main text, we focused on two scenarios that provide useful boundary estimates for clarity and policy relevance. Scenario 1 represents a complete loss of local pollinators, an extreme outcome that illustrates the maximum potential impact on nutrition and livelihoods. Scenario 2 represents a business-as-usual trajectory, derived from projected trends in the native honeybee A. cerana reported by beekeepers in our ten study villages²⁹. In the absence of long-term pollinator monitoring data, we used this projection (33% decline from 2021 to 2030) as a proxy for all pollinator taxa, with the notion that locally reported drivers (climate change and floral resource loss) are equally likely to affect the abundance of wild pollinators^57,58 as for the abundance of A. cerana (details in Supplementary Method 11). Although necessarily a simplification, the rate is consistent with modelled estimates for declines in wild insect elsewhere²⁴ and may even be conservative given that wild pollinators typically decline faster than managed honeybees^59,60.

Simulating pollination recovery

Our third modelled scenario aimed to predict the potential nutritional and economic gains if pollination services were actively managed to close the proportion of each crop’s yield gap that is attributed to insufficient pollination (that is, achieving optimal pollination). We first calculated current yield gaps for four key pollinator-dependent crops (Jumli beans (Phaseolus vulgaris), mustard (Sinapis alba), pumpkin (Cucurbita maxima) and apple (Malus domestica)) that together account for >75% of pollinator-dependent food consumption and >90% of farming income. Yield gaps were defined as the difference between the upper-95 percentile yield and the median yield across farms in our study region. We then assumed that a part of this yield gap (median of 24%) could be closed, drawing on data from a global study of comparable smallholder farms in Asia (including Nepal), Africa and Latin America that quantified the fraction of yield gaps attributable to insufficient pollination³⁶. Notably, we only closed this pollination-related share of the yield gap, whereas the larger proportion caused by other agronomic constraints (for example, fertilizer, water and soil quality) was left unchanged. Following a previously described approach¹⁷, we applied this pollination-related yield increase to each crop, recalculated participants’ nutrient intakes and household incomes and compared these to current values (further details in Supplementary Method 12). Our pollination deficit estimates aligned closely with those reported in regional⁶¹ and global³⁰ assessments, which provided confidence in the validity of this scenario.

Key model assumptions

All simulations were based on the following simplifying but evidence-informed assumptions: (1) the pollinator-dependent yield of each crop declines in proportion to the decline in pollinator abundance, as per a previous study²⁴; (2) only the yields of locally produced pollinator-dependent crops are affected (representing the part of the food system under farmers’ control), whereas imported foods are assumed unchanged; (3) changes in food production lead to equivalent changes in food consumption, as this is a strongly food-limited population for which production and consumption are tightly linked (indicated by the strong seasonal synchrony between production and consumption of micronutrient-rich foods; Supplementary Fig. 7); (4) food lost from the diet is replaced by the main staple crop rice (of an equal number of calories), as this is one of the few staple foods routinely imported and affordable year-round and typically substitutes for local ingredients when unavailable; (5) declines in locally produced pollinator-dependent foods are not offset by imports, as market failures and geographical isolation strongly constrain their availability and affordability⁶² (Supplementary Fig. 17); (6) the market price for exported cash crops (for example, apples and beans) remains fixed, as prices are mostly determined by supply and demand in larger regional markets rather than local production changes; and (7) farmers do not alter cropping patterns in response to declines in pollinator-driven yields because high investment costs and low financial capital limit their ability to adapt. Further details and sensitivity analyses are provided in Supplementary Method 1. Although these assumptions do not capture the full complexity of food systems, they are not improbable in the context of Jumla District. By making the assumptions explicit, we highlight where interventions could mitigate these effects, for example, by strengthening market access to improve the availability of micronutrient-rich foods or by providing economic safety nets that enable farmers to adapt to changing conditions.

Pollinator characteristics that make them nutritionally important

To extend the generality of our findings beyond the specific pollinator taxa present in our study system, we investigated the network characteristics of pollinators that could predict their nutritional importance. We quantified the nutritional importance of a pollinator as the proportional reduction in the intake of six key pollinator-dependent nutrients (vitamin A, folate, vitamin E, calcium, iron and vitamin C) across the entire study population when the pollinator is removed from the crop–pollinator network (same modelling framework as described above). We then tested whether variation in nutritional importance among pollinators could be explained by their abundance and a limited set of species-level network descriptors calculated on the full plant–pollinator network (all villages combined) using the R packages bipartite⁶³ and vegan⁶⁴. Metric selection was based on a previous study⁶⁵ and aimed to represent a small number of complementary dimensions of pollinator behaviour and network position that could plausibly influence crop pollination and nutrient provision while remaining relatively straightforward to interpret for agricultural and nutrition specialists and network ecologists.

Network metrics tested

Pollinator abundance was quantified as the total number of flower visits recorded for a species across all plants, which provided a proxy for its local pollination activity. Abundance was included as a baseline predictor in all models as it is a fundamental ecological property that strongly influences both network structure and function^65,66.

To capture how pollinator activity is allocated between crops and non-crops, we calculated a metric we term ‘crop focus’, which was defined as the proportion of total visits by a pollinator (across all plants) that were directed to crop species. This metric distinguishes pollinators that concentrate their activity on crops from those that primarily visit wild plants, independent of overall activity level.

In addition to pollinator abundance, we characterized the interaction breadth and selectivity of each insect using three complementary metrics: (1) degree, which captures partner range (the number of plant species visited) and is interpreted strictly as a measure of interaction breadth rather than interaction strength⁶⁷; (2) Shannon diversity of interactions, which captures how evenly a pollinator distributes its visits across its plant partners, weighting partners by visitation frequency⁶⁸; and (3) Blüthgen’s d′, a specialization metric that quantifies how selectively a pollinator uses its plant partners relative to their availability in the network⁶⁹. Together, these metrics describe how many plant species an insect visits, how evenly it distributes its visits among those plants and whether its interaction pattern reflects opportunistic use of available resources or clear preferences for particular plant species⁷⁰. We expected that insects with high degree and interaction diversity are nutritionally important pollinators because they interact with a wide range of crops. By contrast, insects with higher d′ values may be particularly important for maintaining pollination of key crops owing to their consistent and selective visitation patterns.

To describe the position of each insect in the network as a whole, we calculated weighted closeness centrality, which was defined as the average distance of a pollinator to all other species in the network⁷¹. Closeness centrality is commonly interpreted as a proxy for the potential of a species to influence network-wide dynamics, which is based on the assumption that direct interactions exert stronger effects than indirect ones and that the strength of indirect effects declines as the number interaction steps increases⁷¹. In this context, centrally positioned pollinators have the potential to affect pollination across multiple crops in the network, although high centrality does not necessarily imply strong or frequent interactions with crop plants.

Because network metrics are influenced by abundance and sampling effects⁶⁵, we converted each observed network metric to a Z score based on 99 null-model simulations that preserved the total number of visits per pollinator while randomizing their distribution across plant species. Null models were generated using the nullmodel function in the R package vegan with the c0_ind method⁶⁴. These Z scores indicated whether a pollinator exhibits higher or lower values of a given network property than expected based on its abundance alone. We then used linear models to test whether pollinator characteristics predict nutritional importance, including abundance as a baseline predictor and testing each Z-standardized network metric for additional explanatory power. Model performance was assessed using changes in Akaike information criterion, adjusted R² and the significance of regression coefficients.

Management of pollinators to enhance human nutrition

One of the most successful and widely practised pollination management strategies is the addition of wild flowering plants to farmland to provide food for pollinators^43,44. Focusing our attention on this common pollination-management approach, we predicted the most important non-crop plants for supporting the key crop pollinators at our field sites. We estimated this ‘indirect contribution’ of each wild plant to the pollination service by calculating its visitation by insects and weighting the visitation of each insect by its PC score (see section above). Because PC is calculated separately for each nutrient, we sum values across all nutrients to obtain an overall contribution for insect (PC_i). Thus, for a wild plant p, our estimate of its indirect contribution (IC) to pollination services is calculated as follows:

where PC_i is derived from the equation in the section described above and V_i,p is the visitation frequency of each insect i to plant p.

Thus, a wild plant that is frequently visited by pollinators of crops that provide a large proportion of peoples’ micronutrient supply will receive a high IC score. This approach assumes that the value of a plant to pollinators is proportional to how often it is visited, without accounting for the nutritional quality of floral resources or how much is consumed during visits. Our estimates should therefore be interpreted as indicative rather than precise. To validate this approach, we investigated the phenological pattern of crop and wild plant visitation to check that wild plants are supporting crop pollinators with floral resources outside the main periods of crop flowering and therefore are facilitating crop pollination rather than simply competing with crops for pollinators (Supplementary Fig. 16).

It should be noted that in addition to the high-scoring plants identified through this metric, there may be some less widely used (and therefore lower scoring) plants that fill temporal or nutritional gaps in food supply and are therefore more important than their score would suggest. For the purposes of this analysis, however, we focus on the most widely used plants that are assumed to provide insects with most of their key resources.

Data analysis

Data cleaning and processing were conducted using STATA (v.18), whereas all data analyses, modelling and figure plotting were performed in R (v.4.4.2, R Core Team).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.