Human land use drives microplastic abundance in terrestrial wildlife

Sampling areas

Lowland tapirs

Lowland tapir (T. terrestris) samples were collected across three sites in Brazil. The agricultural mosaic site in the Cerrado was located between the municipalities of Nova Alvorada do Sul and Nova Andradina, in Mato Grosso do Sul State, Brazil (−21.99° N, −53.83° E). It features a mosaic of natural habitats, including Cerradão fragments, gallery forests and marshlands, which constitute approximately 12% of the study area. The remainder is heavily influenced by human activities, including agriculture (for example, sugarcane, soybean and corn), cattle ranching, eucalyptus plantations, highways and rural settlements. Precipitation exhibits marked seasonality, with approximately 70% of the total annual rainfall occurring between November and March⁵¹.

The field site in the Pantanal was situated on Baía das Pedras Ranch, in the Nhecolândia Sub-Region of the Southern Pantanal, Municipality of Aquidauana, Mato Grosso do Sul State, Brazil (−19.33° N, −55.72° E). This private property spans 145 km² and comprises a mosaic of seasonally inundated grasslands, lakes, gallery forests, scrub, and deciduous forests. The Pantanal is a nearly pristine biome with a low human footprint, and is far removed from areas experiencing deforestation and other anthropogenic pressures. The ranch supports abundant wildlife and extensive cattle grazing on native grasses. The average annual temperature is 25 °C, and the annual rainfall is 1,185 mm.

The site in the Amazon was located at Tanguro Farm, a privately owned 103,000 ha property in Querência, Mato Grosso state, that sits at the transition between the Cerrado savanna and the Amazon rainforest along the southern Arc of Deforestation. Owned by AMAGGI, a major actor in Brazil’s soybean industry, the farm comprises 59,000 ha (57%) of croplands (soybean, corn, and cotton) and 44,000 ha (43%) of legally protected forest reserves. This region is the driest portion of the Amazon basin and experiences a prolonged dry season (May–September) with annual rainfall of 1,700–2,200 mm. Forests in this region are shorter (~20 m) than the moist Amazon rainforest (~30 m). About 40% of the farm was cleared for pasture (1976–1992) and later converted to soy cultivation (2003–2010), leaving riparian forests surrounded by croplands.

Giant anteaters

The giant anteater (M. tridactyla) data collection was conducted in the Cerrado savanna of Mato Grosso do Sul state, Brazil, as part of the Anteaters and Highways Project (ICAS, www.icasconservation.org.br), which investigates wildlife-vehicle collisions and their impacts on wildlife populations. Sampling areas were located near two federal roads: BR-267 (−21.64° N, −53.51° E) and BR-262 (−20.86° N, −52.04° E). As described above, the study region is a mosaic landscape dominated by open areas, including grasslands, pastures, and croplands, interspersed with small remnants of native savanna and forest vegetation, as well as eucalyptus plantations⁵².

Giant armadillos

In addition to the tapir and giant anteater data, we also obtained a smaller amount of data from three giant armadillos (P. maximus). The giant armadillo data collection was also conducted at the Baía das Pedras Lodge and Ranch in the Pantanal, which served as the base of operations for the Giant Armadillo Conservation Program (GACP).

Animal captures and GPS data acquisition

Ethics statement

All animal capture, handling, and blood sampling procedures were approved by the Brazilian Ministry of the Environment (MMA) and the Chico Mendes Institute for Biodiversity Conservation (ICMBio) under license numbers 14603-25 (lowland tapirs), 53798-18 (giant anteaters) and 27587-1 (giant armadillos). All procedures were performed in accordance with the terms of these licenses and with relevant national and institutional animal welfare guidelines. Chemical immobilization was performed under veterinary oversight, and all animals were released at their capture sites following full recovery.

Lowland tapirs

Lowland tapirs were captured using a combination of methods, including darting, following physical restraint in box-traps or pitfall traps, or darting from a distance. Box traps and pitfall traps were deployed in areas frequently visited by tapirs. Tapirs captured by these methods were anaesthetized intramuscularly with a combination of butorphanol (0.17 mg kg⁻¹), medetomidine (0.012 mg kg⁻¹) and ketamine (0.7 mg kg⁻¹). Reversal agents, including atipamezole (0.04 mg kg⁻¹) and naltrexone (0.35 mg kg⁻¹), were administered at the conclusion of procedures. For tapirs darted from a distance, anaesthetics included tiletamine-zolazepam (1.25 mg kg⁻¹), medetomidine (0.006 mg kg⁻¹), ketamine (0.6 mg kg⁻¹) and atropine (0.03 mg kg⁻¹). Darts (3 ml or 5 ml) were used to deliver anaesthetics. Midazolam (0.03 mg kg⁻¹) was administered 30 min later to maintain anaesthesia. Captured tapirs underwent a range of procedures, including subcutaneous microchip insertion, morphometric measurements, sex and age determination, physical examinations, blood sample collection and telemetry collar placement for adults.

Tapirs were equipped with GPS collars (Telonics TGW SOB and GPS IRIDIUM models) programmed to obtain hourly location fixes, with a median operational life of 15.4 months. Additional location data were collected from stationary tags placed in the Pantanal (n = 883) and Cerrado (n = 174) to calibrate the measurement error in the GPS datasets.

Giant anteaters

Giant anteaters were captured between February 2021 and November 2022 using motorized active searches between 15:00 and 17:00, when they exhibit higher diurnal activity due to colder temperatures^53,54. Searches were conducted along dirt roads in open areas, with two team members stationed in the back of a truck to improve detectability. Upon spotting an individual, technicians approached the animal on foot and captures were performed using two long-handled dip nets (1.5 m handles, 0.7 m diameter hoops). Physical capture typically lasted around 1 min.

Captured giant anteaters were chemically immobilized, equipped with GPS radiotelemetry harnesses (Telonics) as previously described⁵⁵ and blood samples were collected. In brief, animals were anaesthetized using a protocol of butorphanol tartrate (Zoetis Indústria de Produtos Veterinários; 0.1 mg kg⁻¹), detomidine hydrochloride (Agener-União; 0.1 mg kg⁻¹) and midazolam hydrochloride (Cristália; 0.2 mg kg⁻¹), adapted from ref. ⁵⁵. The induction period averaged 8.27 ± 4.78 min.

Giant armadillos

Giant armadillos were captured between February 2021 and November 2022, using a method tailored to their natural behaviour⁵⁶. Researchers tracked footprints and signs of recent predation, such as termite mounds, to locate burrows. When a burrow was identified and confirmed to contain an armadillo (indicated by the entrance being sealed with sand), an iron funnel trap (jequi-type) was placed at the entrance. Once captured, animals were transported in iron-reinforced wooden crates and held overnight for nocturnal species, with procedures beginning the following morning at the capture site.

Captured giant armadillos were anaesthetized for the collection of blood samples and surgical implantation of intra-abdominal transmitters. Anaesthetic protocols included an intramuscular injection of butorphanol tartrate (0.1 mg kg⁻¹), detomidine hydrochloride (0.1 mg kg⁻¹), and midazolam hydrochloride (0.1 mg kg⁻¹) into the hind limbs. Before surgery, intradermal local anaesthesia with lidocaine hydrochloride (2 mg kg⁻¹, Eurofarma) was infiltrated into the infra-umbilical region. Ketamine hydrochloride (10 mg kg⁻¹) was administered intramuscularly as a supplementary anaesthetic immediately after the animals were completely washed. Transmitters (IMP 310 model, Telonics; 9.4 × 2.3 cm, 42 g, 36-month operational life) were surgically implanted within the abdominal cavity under sterile conditions. Following surgery, anaesthesia was reversed using a combination of naloxone hydrochloride (0.02 mg kg⁻¹), yohimbine hydrochloride (0.125 mg kg⁻¹), and flumazenil (0.01 mg kg⁻¹). Armadillos were released at their capture sites upon recovery, typically remaining in their burrows for 24 to 120 h (average of 48 h) before emerging. Daily monitoring was conducted for 12–22 days post-release using telemetry receivers (model TR-4, Telonics) and VHF antennas (model RA-23K). After the initial monitoring period, armadillos were monitored monthly until the end of the transmitter’s operational life.

Blood sample collection

Lowland tapirs

Approximately 10 ml of blood was collected within 20 min of immobilization via venipuncture of the medial saphenous or cephalic veins using a sterile glass vacutainer system (BD Vacutainer). Tapir captures were conducted in remote forest sites without access to field refrigeration. Blood samples were therefore stored in pre-cleaned glass collection tubes with no preservatives, sealed with clean aluminium foil and a metal cap to prevent contamination during storage, kept shaded at ambient temperature, and transferred to a −20 °C freezer on return to the base camp within 5–6 h.

Giant anteaters and giant armadillos

Blood samples (~5 ml) were collected approximately 20 min after induction from the cephalic or saphenous veins using a sterile vacutainer system with 23G or 25G needles (BD Vacutainer). Samples were transferred to pre-cleaned glass collection tubes with no preservatives, sealed with clean aluminium foil and a metal cap to prevent contamination during storage, kept refrigerated at 4 °C on site, and frozen at −20 °C within 12–24 h for subsequent processing. We note that blood samples were obtained from three individual giant armadillos, however we obtained repeat samples at two time points for two of these animals. Although these blood samples were processed individually, they were pooled (that is, averaged) in our statistical analyses.

Cold storage under the above-mentioned conditions is standard practice in microplastic analyses and does not affect particle abundance or polymer identity. Microplastics are chemically stable over these timescales, and method reviews emphasize that environmental samples containing microplastics are routinely stored at 4 °C for days to weeks prior to processing, or frozen at −20 °C for longer periods, without alteration of polymer composition (reviewed in ref. ⁵⁷).

Microplastic isolation from blood

Blood microplastic isolation was performed following a protocol previously developed by our team^58,59,60. Depending on the sampling date, blood samples were stored at −20 °C for ~1–22 months between collection (Feb 2021–Nov 2022) and digestion (Dec 2022). On the first day of isolation, 1 ml of each blood sample was placed in an Erlenmeyer flask along with 50 ml of 10% KOH solution (sample ratio of 1:50). The mixture was incubated on a shaker at 60 °C and 250 rpm for 24 h. After the initial 24-h incubation, 10 ml of 47% NaClO were added to each sample to achieve a final concentration of 7.7%. The samples were then incubated under the same conditions (60 °C and 250 rpm) for an additional 24 h.

Following the second day of digestion, all samples were filtered through a 47 mm polytetrafluoroethylene (PTFE) membrane with a 0.45 μm pore size (Merck Omnipore) and rinsed 3 times with filtered ultrapure water. The membranes containing the particles were then submerged in 50 ml of 20% HNO₃ and placed in an ultrasonic bath with sweep function at 100% power and 45 kHz frequency for 15 min to transfer the particles from the membrane to the solution. The samples were further incubated for 30 min at 40 °C and 250 rpm.

After removing the membranes, the samples underwent a 24-h incubation at 40 °C and 250 rpm. On the final day, the samples were filtered through a 13 mm PTFE membrane with a 0.45 μm pore size, rinsed three times with filtered ultrapure water, and mounted on glass slides for microscopy and spectroscopy analysis. The resulting filter membranes were then stored in pre-cleaned glass Petri dishes and analysed by Raman spectroscopy over the subsequent ~12 months (Jan–Dec 2023). After digestion, a proportion of non-plastic particulates remained on the membrane, such that microplastics accounted for an average of 16.4% of imaged particles with HQI > 75 (range 0–89.9%).

Quality assurance and quality control methods used throughout the process are detailed in the Supplementary Methods (see also Extended Data Figs. 6 and 7).

Confocal Raman analysis of particles

Microplastic identification and characterization were conducted using a confocal Raman microscope (alpha300 R, WITec) equipped with a UHTS 300 spectrometer (WITec) and a 532-nm laser. The microscope was operated using Control Five software, and the acquired data were processed with Project Five software to compensate for background noise and cosmic radiation signals.

As an initial step, a whole membrane image was captured using a 10× NA 0.25 Zeiss objective with the Area Stitching option. From this image, a specific field was selected for detailed imaging with a 50× NA 0.75 Zeiss objective, employing both the Area Stitching option and Z-stack autofocus. Particle identification was facilitated by the Particle Scout software, which enabled the acquisition of Raman spectra for all identified particles across the entire membrane. Raman spectra were then obtained using a 100× DIC NA 0.9 Zeiss objective, with a laser power of 10 mW, an accumulation of 4, and an integration time of 0.45 s. Spectra autofocus was used throughout the measurements. The acquired spectra were then compared against the S.T. Japan database (S.T. Japan Europe GmbH) and the SLoPP MPs Raman library⁶¹ using the True Match-Integrated Raman Spectra Database Management software. Microplastic characterization was based on an automated, 3-step process that was optimized for spectra acquisition rather than image sharpness. This precluded the manual assessment of particle shape (for example, fibres or fragments) or colour. Only matches with an HQI greater than 75% were considered for further analysis.

HFI differences between 2001 and 2020

To assess changes in the HFI over time, we obtained and processed global HFI rasters at a 300 m resolution for the years 2001 and 2020 from the Wildlife Conservation Society Human Footprint dataset (https://wcshumanfootprint.org/). After ensuring projections and resolutions were consistent between the HFI rasters, a difference raster was generated by subtracting the 2001 HFI from the 2020 HFI to visualize areas of increasing and decreasing human influence. We cropped and masked the HFI rasters to remove areas outside of Brazil and then computed the mean HFI for 2001 and 2020, as well as the absolute differences. These analyses were conducted in R using the terra (v1.8-54)⁶², sf (v1.0-21)^63,64, tmap (v4.0.0.9000)⁶⁵, fields (v16.3)⁶⁶, rnaturalearth (v1.0.1)⁶⁷, and geobr (v1.9.1)⁶⁸ packages.

Animal movement analysis

Movement data pre-processing

Before analysis, we pre-processed the GPS location data in order to calibrate the measurement error and filter outliers. Pre-processing was done via the methods implemented in the R package ctmm (v1.2.1)⁶⁹. For each location estimate, the GPS trackers recorded a unitless Horizontal Dilution of Precision (HDOP), a value which is a measure of the accuracy of each positional fix. We converted the HDOP values into calibrated errors by specifying 20-m error at an HDOP of 1 for 2D fixes, and 10-m error at an HDOP of 1 for the 3D fixes⁷⁰. For each individual dataset, we then removed outliers based on error-informed distance from the median location, and the minimum speed required to explain each location’s displacement. Through this process we identified a total of 16,656 outliers, which represented ca. 3.55% of the total 469,410 locations. After filtering outliers, the final dataset consisted of 343,913 GPS locations across 20 individual giant anteaters, 6,072 GPS locations across 2 individual giant armadillos, and 102,769 GPS locations across 20 individual lowland tapirs.

Movement modelling

For each of the tracked animals, we used a continuous-time approach to estimate their movement rates and home ranges, enabling robust comparisons across individuals⁷¹. We first fitted a series of continuous-time stochastic process movement models to the location data. Parameters were estimated using perturbative-hybrid residual maximum likelihood (pHREML)⁷², and identified the best model via small-sample-sized-corrected Akaike’s information criterion (AICc). From both visual inspections of the semi-variograms, and the parameters of the fitted movement models, all animals occupied stable home ranges throughout the monitoring period with no indication of range shifts or large-scale dispersal events.

From the selected models, we estimated each animal’s diffusion rate (in km² day⁻¹) as a measure of their movement rate. We next estimated each animal’s home range, using autocorrelated kernel density estimation (AKDE)⁷³. AKDE home range estimates were conditioned on the autocorrelation structure of the fitted movement model, and we implemented the small-sample-size bias correction of Fleming & Calabrese (2017)⁷⁴, as well as the optimal weighting of Fleming et al. (2018)⁷⁵ to correct for uneven sampling. Home range areas (in km²) were estimated as the area contained within the polygon delimited by the 95% isopleth of the probability density.

We estimated the mean and maximum HFI within each individual’s home range. HFI is an aggregate measure of human disturbance that ranges from 0 in natural habitats to 50 in the most anthropogenically modified habitats³³. A raster of HFI data with a spatial resolution of 100 m was obtained from Gassert et al.³⁸. The weighted mean HFI within an individual’s home range was estimated by calculating the weighted arithmetic mean of HFI values:

with weights (w_u) determined by the probability density at each location u. Maximum HFI exposure was estimated by identifying the maximum value of HFI within each animal’s home range. We also estimated more nuanced patterns of habitat use via the MapBiomas database³⁴, which details land use throughout Brazil at a spatial resolution of 10 m. We focused on habitats classified as water, development, agriculture, pasture, and native forest, which represented the majority of the land classes within the monitored animals’ home ranges. The proportions of each of these land classes within each animal’s home range were again estimated by calculating the weighted arithmetic mean. Lastly, we used the MapBiomas land classification data, the Rocky Mountain Institute (RMI), and the Waste Methane Assessment Platform 1.5.0, 2025 (https://wastemap.earth) to generate spatial rasters of distance to development, and distance to the nearest waste management facility (both in m). Development was a compound measure that captured the built environment, and was composed of the MapBiomas collection 10 classes: urban areas (24), mining (30), photovoltaic power plants (75), and other non-vegetated areas (25). From these spatial data, we estimated each animal’s mean distance to development and to the nearest waste management facility, estimated by calculating the weighted arithmetic mean.

Statistical modelling

To study how human land-use change shaped blood microplastic concentrations, we fit hierarchical generalized linear regression models using the methods implemented in the mgcv package (v1.9-1)⁷⁶. Because the microplastic concentrations scaled between [0,∞], we used a Tweedie distribution with a log link for the error distribution⁷⁷, with parameters estimated via restricted maximum likelihood (REML). The regression models testing for correlations between individual habitat use and microplastic concentrations included a population-level effect, as well as random intercepts and slopes for each of the species, with the following general form

where β₀ represents the population-level intercept, β is the population-level slope, u_s[i] is the species-specific random intercept, v_s[i] is the species-specific random slope, and s[i] indexes the species of individual i. This structure was used for all of the habitat covariates x_i. From these models, inference was then made on the population-level trend.

All analyses were conducted in R (v4.3.3). Weighted means and proportions of land-use exposure were computed with the weights package (v1.0.4). Figures were produced with ggplot2 (v3.5.2), ggspatial (v1.1.9), tidyterra (v0.7.0), rphylopic (v1.5.0; animal silhouettes), fmsb (v0.7.6) and gridExtra (v2.3).

Statistics and reproducibility

No statistical method was used to predetermine sample size. Sample sizes were determined by the number of animals that could be safely captured during the study period (February 2021 to November 2022) given the logistical constraints of live-capturing large, threatened, free-ranging mammals under national permits, and all animals captured during this period were sampled. This was an observational field study, so randomization was not applicable; animals were not allocated to experimental groups, and exposure to human land use was determined by each individual’s own space use. The investigators were not blinded to species or sampling site during laboratory processing or analysis; however, particle detection and polymer identification were automated (Particle Scout and True Match software) and based on a fixed spectral-match threshold (HQI > 75%), which limits observer bias. Each animal contributed one blood sample, with the exception of two giant armadillos sampled at two time points, whose values were averaged before analysis. Seven animals lacked sufficient GPS tracking data and were excluded from the habitat use models but retained in descriptive analyses. Samples were processed in two batches, each with a procedural blank (Supplementary Methods); all lowland tapir samples were processed in the same batch. The polymer-integrity test shown in Extended Data Fig. 6 was repeated three times independently with similar results.

Reporting summary

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