Shik et al.

Authors and contact

Corresponding author: Prof. Jonathan Shik
Email: jonathan.shik@bio.ku.dk

Overview

This repository contains the data, derived bootstrap outputs, R scripts, and helper functions used for the analyses in Nutritional niches reveal mutualism resilience in urban fungus-farming ants.

Data files

• FNN_Data_raw.csv — Fundamental Nutritional Niche (FNN) data describing fungal cultivar performance across diets varying in protein and carbohydrate. The response variables include final hyphal growth area (mm²) and staphyla density (number per mm²).
• RNN_Data_raw.csv — Realized Nutritional Niche (RNN) data describing the elemental and macronutrient composition of foraged plant fragments. These data are used to estimate RNN intake targets (RNN ITs) and convex-hull RNN areas.
• Shiketal_data.xlsx — Workbook containing five sheets: pie (data for pie charts), col_hab (foraging data summarized by colony and sampling interval), col_hab_model (foraging data in long format for fragment-type analyses, plus an ancillary wide-format summary block), nutrient_data (elemental and macronutrient composition of plant fragments), and glycogen (glycogen measurements from field-collected fungal cultivar).

Derived bootstrap files

• Shik_Area_boot_RNN_data2.csv — Bootstrapped RNN convex-hull points and areas for the elemental plane.
• Shik_Area_boot_RNN_data2_macro.csv — Bootstrapped RNN convex-hull points and areas for the macronutrient plane.
• Shik_RNN_IT_boot_data2.csv — Bootstrapped RNN intake-target estimates.

R scripts and helper functions

• Shiketal_2026_Code_GFN_Final_V2.R — Main analysis script for constructing FNN performance landscapes, bootstrapping RNN hulls and intake targets, quantifying overlap between habitats, and comparing realized nutrient intake with fungal performance landscapes.
• General_Functions_04-01-2024.R — Helper-function library used for nutritional geometry analyses. These functions are sourced by the main analysis script.
• Shiketal_code.R — Script for analyzing Shiketal_data.xlsx and producing the figures and statistical analyses based on those data.

Repository structure

.
├── README_6Nov2026.md
├── FNN_Data_raw.csv
├── RNN_Data_raw.csv
├── Shiketal_data.xlsx
├── Shik_Area_boot_RNN_data2.csv
├── Shik_Area_boot_RNN_data2_macro.csv
├── Shik_RNN_IT_boot_data2.csv
├── Shiketal_2026_Code_GFN_Final_V2.R
├── Shiketal_code.R
└── General_Functions_04-01-2024.R

Software requirements

R version 4.2 or later is recommended.

Packages used by the main analysis script

dplyr, ggplot2, tidyr, pracma, Hmisc, alphahull, purrr, terra, agricolae, patchwork, lmerTest, stringr, viridis, MuMIn, and effsize.

Packages used by Shiketal_code.R

tidyverse (including ggplot2, dplyr, tidyr, readr, stringr, purrr, and scales), readxl, nlme, car, emmeans, ggh4x, colorspace, ggforce, performance, gridExtra, moments, and lmerTest.

Randomness and reproducibility

The main analysis script sets the random seed to 321321.

Bootstrap settings:

• RNN bootstrap: 1,000 repetitions.
• Additional bootstrap for confidence intervals of mismatch distances: 10,000 repetitions.

Script 1: FNN and RNN analyses

Input data: FNN_Data_raw.csv

Purpose. Raw data used to construct FNN performance landscapes across diets varying in protein and carbohydrate.

Read into R with:

fnndt <- read.csv("FNN_Data_raw.csv", header = TRUE)

Key variables:

• colony_ID: colony or strain identifier, used to assign the duration of data collection.
• dataset: study label used to distinguish the Shik et al. dataset (Shik_NEE) from the other study dataset.
• pctprot: dietary protein (% by mass).
• pctcarb: dietary carbohydrate (% by mass).
• FinalArea: final hyphal growth area (mm²).
• staph_dens: staphyla density (number per mm²).
• prot: protein concentration (g/L), used for supplementary concentration-based figures.
• carb: carbohydrate concentration (g/L), used for supplementary concentration-based figures.

Derived variables.

The script creates days_col from colony_ID as follows: Ac_11 = 44 days, Ac_9 = 73 days, Ac_10 = 71 days, and all other colony IDs = 56 days. It then calculates growth rate as FinalArea / days_col.

Important subsets.

• fnndt_biomass8: records where dataset == "Shik_NEE".
• fnndt_others: records where dataset != "Shik_NEE".
• fnndt_full: the full dataset.

Outputs from the FNN analysis.

The script generates predicted thin-plate spline surfaces (objects named FNN_land_*, typically containing PROT, CARB, and tps columns). PROT and CARB are the protein and carbohydrate coordinates of the prediction grid, and tps is the predicted response. Peak-performance regions are identified using predict_region_GF_lm() and converted to polygons using design_hull(). Peak-region centroids are calculated using calculate_centroid_boot().

Input data: RNN_Data_raw.csv

Purpose. Raw field data used to quantify realized nutritional niches across sites and plant species or tissue types, and to estimate intake targets using dry-mass weighting.

Read into R with:

rnndt <- read.csv(
  "RNN_Data_raw.csv",
  header = TRUE,
  strip.white = TRUE,
  na.strings = "n/a"
)

Pre-processing.

Tissue-type labels are converted to title case. The script also standardizes labels as follows: Leaf_petiole becomes Petiole, Leaf_stipule becomes Stipule, and Fruit_2 becomes Fruit.

Key variables:

• Site: site identifier; GAM = Gamboa (urban) and PLR = Pipeline Road (rainforest).
• type: plant tissue type.
• sp: plant species code or name.
• Tot_fragment_no: total number of fragments recorded for the species.
• Sum_dry_weight..g.: total dry mass (g), as read by R after standardization of the source column name.
• pct_N and pct_C: nitrogen and carbon content (% by mass).
• pct_protein: protein content (% by mass), estimated from nitrogen content.
• pct_TNC: total non-structural carbohydrate content (% by mass), calculated as the sum of starch and soluble sugar.

Derived quantities.

Nutrient values are multiplied by dry mass to obtain mass-weighted nutrient values. These are summed by site, tissue type, and species, then divided by the corresponding total dry mass to obtain mass-weighted nutrient composition for each site–tissue type–species combination. Site-level intake targets (IT_N, IT_C, IT_protein, and IT_TNC) are calculated as the means of these mass-weighted values across tissue-type–species combinations within each site.

Missing values in the RNN data

The pct_TNC column contains 30 blank cells among the 52 records in RNN_Data_raw.csv. These values indicate that a TNC value is unavailable for the corresponding plant-fragment record due to insufficient biomass for wet chemistry analyses (described in detail in the manuscript Methods); they do not represent zero or a value below detection. The same 30 records thus also have missing pct_TNC values in the nutrient_data sheet of Shiketal_data.xlsx. The elemental measurements (pct_N and pct_C) and protein estimates (pct_protein) are present for these records.

Analyses requiring TNC use only records with an available TNC value. In particular, the macronutrient convex-hull bootstrap excludes records with missing values. The missing TNC values are retained as missing in the deposited data rather than being replaced with zero or an estimated value.

Derived bootstrap output: Shik_Area_boot_RNN_data2.csv

Purpose. Bootstrapped convex-hull points and hull areas for each site in the elemental plane (%N versus %C).

Within each bootstrap replicate, records are sampled with replacement within site. Convex-hull points are identified from pct_N and pct_C, and hull area is calculated using terra::expanse() on the convex hull. The output includes site, hull-point coordinates, hull area (area_chull), and bootstrap replicate identifier (rep).

Derived bootstrap output: Shik_RNN_IT_boot_data2.csv

Purpose. Bootstrapped estimates of RNN intake targets for each site.

Within each bootstrap replicate, records are sampled with replacement within site and the intake-target calculation is repeated. The output includes Site, IT_N, IT_C, IT_protein, IT_TNC, their corresponding standard deviations where returned by the function, and the bootstrap replicate identifier (rep).

Derived bootstrap output: Shik_Area_boot_RNN_data2_macro.csv

Purpose. Bootstrapped convex-hull points and hull areas for each site in the macronutrient plane (%protein versus %TNC).

Before sampling, records with missing values are omitted. Convex-hull points are identified from pct_protein and pct_TNC, and hull area is calculated using terra::expanse() on the convex hull. The output includes Site, pct_protein, pct_TNC, area_chull, and rep.

Analysis workflow

A. FNN performance landscapes

The script estimates continuous performance landscapes across protein–carbohydrate space for staphyla density (staph_dens), final hyphal growth area (FinalArea), and growth rate (FinalArea / days_col). It loads the FNN data, derives collection duration and growth rate, and analyzes the Shik et al. dataset, the other study dataset, and the combined dataset. Thin-plate spline surfaces are fitted using predict_surface(). The high-performance region is identified using predict_region_GF_lm() with a 90 % region threshold, and the region is converted to a hull polygon using design_hull(). The surfaces and peak regions are plotted using ggplot2.

Mixed-effects models are also used to test how performance traits vary with diet, with study dataset and colony identity included as random effects.

B. Standardization of FNN traits

FNN response traits are standardized within each study to make performance values comparable across datasets. The standardize_shik() function creates standardized variables (std_staph_dens, std_FinalArea, and std_growth_rate). Thin-plate spline surfaces and peak regions are then calculated for the standardized traits.

C. RNN bootstrapping and overlap

The analysis quantifies realized nutritional niche breadth (convex-hull area) and intake targets for each site, with uncertainty estimated by bootstrap resampling. Elemental hulls are calculated in the %N–%C plane, and macronutrient hulls are calculated in the %protein–%TNC plane using records with available values for both variables.

The analysis calculates the intersection area between the GAM and PLR convex hulls for each bootstrap replicate and expresses overlap relative to each site's hull area. It also calculates the proportion of bootstrapped intake-target estimates that fall within the overlap region between the two sites and within the overlap between an RNN and an FNN high-performance region.

D. FNN–RNN comparison and correspondence distance

The analysis overlays RNN points, hulls, and intake-target estimates on standardized FNN performance landscapes. For each FNN high-performance hull, the script calculates its centroid. It then calculates Euclidean distances between the FNN centroid and bootstrapped RNN points and compares these distances between sites using Wilcoxon or Kruskal–Wallis tests and Cliff's delta effect sizes.

Helper functions

The main analysis script sources General_Functions_04-01-2024.R. The functions called directly include:

• predict_surface()
• predict_region_GF_lm()
• design_hull()

Other helper functions in this library are also used by the analysis workflow. The main analysis script requires this helper-function file to run.

Running Script 1

Place the following files in the same working directory:

• FNN_Data_raw.csv
• RNN_Data_raw.csv
• Shik_Area_boot_RNN_data2.csv (optional if it is to be regenerated)
• Shik_RNN_IT_boot_data2.csv (optional if it is to be regenerated)
• Shik_Area_boot_RNN_data2_macro.csv (optional if it is to be regenerated)
• Shiketal_2026_Code_GFN_Final_V2.R
• General_Functions_04-01-2024.R

Set the R working directory to this folder and run:

source("Shiketal_2026_Code_GFN_Final_V2.R")

The script creates ggplot objects in the R session, including FNN landscape panels for the individual and combined datasets, standardized FNN landscapes, RNN scatterplots and hulls in elemental and macronutrient planes, and combined FNN–RNN overlay and mismatch plots. The script does not automatically save figures to disk; use ggsave() with explicit filenames and dimensions if exported figure files are needed.

Script 2: Analyses of foraging and fungal glycogen data

Workbook: Shiketal_data.xlsx

The workbook contains five sheets.

pie sheet

Purpose. Data used to construct the pie charts (Fig. 2A), showing the relative biomass of foraged plant species within each habitat.

Read into R with:

pie <- readxl::read_excel("Shiketal_data.xlsx", sheet = "pie")

Key variables:

• site: habitat (GAM = Gamboa, urban; PLR = Pipeline Road, rainforest).
• type: plant-fragment category (leaf, flower, or fruit); Ficus sp. also includes petiole and stipule tissues.
• sp: plant species identification, determined by DNA barcoding.
• sp_type: concatenation of the type and sp variables.
• tot_fragment_no: total number of fragments for each species.
• mass_g: total fragment mass (g) for each species.
• percent_mass: relative mass of each species as a proportion of the total foraged mass in that habitat.
• pct_mass_tot: percent_mass multiplied by 100.

The script summarizes species by site and produces one pie chart per site. Wedges are ordered by relative abundance and use habitat-specific color gradients (orange for rainforest and purple for urban). Number labels restart for each site and correspond to Table 1.

col_hab_model sheet

Purpose. Long-format foraging data used to construct the stacked bar chart of plant-fragment type relative abundance (Fig. 2B).

Read into R with:

m <- readxl::read_excel("Shiketal_data.xlsx", sheet = "col_hab_model") |>
  dplyr::select(dplyr::where(~ !all(is.na(.))))

Key variables:

• site: habitat (GAM = Gamboa, urban; PLR = Pipeline Road, rainforest).
• col_id: abbreviated colony identifier.
• day: sampling day; each colony was sampled on two non-consecutive days.
• sampling_period: sampling interval within a day; each colony was sampled during two 30-minute periods separated by a break.
• Type: fragment category (leaf, flower, fruit, or other).
• worker_no: number of workers observed during a 10-minute observation period on the foraging trail before fragment collection.
• frag_no: number of fragments collected during the 30-minute sampling interval.
• frag_mass: total dry mass (g) of fragments collected during the 30-minute sampling interval.

Additional columns and blank cells. This worksheet also contains an ancillary wide-format summary block in columns to the right of the long-format data. The block has eight variable rows (id) and separate columns for individual colonies (gam_* and plr_*). Blank cells outside this eight-row block are structural: those rows do not correspond to entries in the ancillary summary block. The fully empty spacer columns between the long-format data and this block are also blank by design. These cells are not missing observations in the long-format variables used for the analyses.

The script uses the long-format data to produce stacked bar charts comparing the relative biomass of leaf, flower, fruit, and other fragments among colonies (Fig. 2B), and to perform mixed-effects analyses of fragment number and mass by fragment type and habitat.

col_hab sheet

Purpose. Foraging data summarized by colony and sampling interval, used for figures showing fragment number and mass (Fig. S1A–B) and for mixed-effects analyses.

Read into R with:

m <- readxl::read_excel("Shiketal_data.xlsx", sheet = "col_hab")

Key variables:

• site: habitat (GAM = Gamboa, urban; PLR = Pipeline Road, rainforest).
• col_id: abbreviated colony identifier.
• day: sampling day; each colony was sampled on two non-consecutive days.
• sampling_period: sampling interval within a day; each colony was sampled during two 30-minute periods separated by a break.
• date: sampling date (day.month.year).
• field_col_id: detailed field colony identifier (collector initials, date, and number).
• worker_no: number of workers observed during a 10-minute observation period before fragment collection.
• leaf_frag_no, fruit_frag_no, flower_frag_no, and other_frag_no: numbers of leaf, fruit, flower, and non-plant fragments collected during the 30-minute interval.
• tot_frag_no: total number of plant fragments collected during the 30-minute interval.
• flower_mass_g, fruit_mass_g, leaf_mass_g, and other_mass_g: dry mass (g) of the corresponding fragment categories collected during the 30-minute interval.
• tot_mass_g: total dry mass (g) of plant fragments collected during the 30-minute interval.

The other category is excluded from the mixed-effects analyses of plant-fragment number and mass.

nutrient_data sheet

Purpose. Plant-fragment elemental data used for figures showing %C and %N (Fig. S2A–B) and analyses of elemental composition.

Read into R with:

nutrient_data <- readxl::read_excel("Shiketal_data.xlsx", sheet = "nutrient_data")

Key variables:

• site: habitat (GAM = Gamboa, urban; PLR = Pipeline Road, rainforest).
• type: plant-fragment category.
• sp: plant species.
• Tot_fragment_no: total number of fragments collected for each species across sampling intervals.
• mass_g: total dry mass (g) of fragments collected for each species across sampling intervals.
• pct_N: nitrogen content (% by mass).
• pct_C: carbon content (% by mass).
• pct_protein: protein content (% by mass), estimated from nitrogen content.
• pct_TNC: total non-structural carbohydrate content (% by mass), calculated as the sum of starch and soluble sugar.
• N and C: nitrogen and carbon mass values used in the analyses.

The pct_TNC column has 30 blank cells among the 52 records. These indicate that a TNC value is unavailable for the corresponding plant-fragment record due to insufficient biomass for wet chemistry analyses; they are not zero values. The records retain their measured elemental and protein values. Analyses of TNC use only records with available TNC measurements.

glycogen sheet

Purpose. Data used to test for differences in glycogen content of field-collected fungal cultivar between habitats and fungal tissue types (Fig. S3).

Read into R with:

dfr <- readxl::read_excel("Shiketal_data.xlsx", sheet = "glycogen")

Key variables:

• Sample_ID: parent colony and technical-replicate identifier.
• Colony_ab: abbreviated colony identifier, used to determine habitat (Gamboa or Pipeline Road).
• chamber: excavated chamber within the colony from which the fungal sample originated.
• TechRep: technical replicate; each fungal sample was analyzed twice.
• Tissue: fungal tissue type (mycelium or staphyla).
• Person: researcher who generated the data point.
• Tissue_mass: dry mass of the analyzed sample (mg).
• Run: 96-well plate on which the extracted sample was analyzed.
• Gly_per_Tissue: glycogen concentration (µg/mL) normalized by mg dry tissue mass.

Derived variables and processing.

Habitat identity (Gamboa or Pipeline Road) is derived from Colony_ab. A small number of marginally negative glycogen values were set to zero; these occurred when glycogen or sample mass values were low relative to the calibration curves. Gly_per_Tissue was then log10(x + 1) transformed for analysis.

The script fits a mixed-effects model testing for differences in transformed glycogen concentration between habitats and tissue types, with technical replicate nested within chamber and chamber nested within colony as random effects. The results are presented in Fig. S3.

Running Script 2

Place Shiketal_code.R and Shiketal_data.xlsx in the same working directory. Set the R working directory to this folder and run:

source("Shiketal_code.R")

The script creates pie charts (Fig. 2A), stacked bar charts (Fig. 2B), boxplots of fragment number and mass (Fig. S1A–B), boxplots of fragment %C and %N (Fig. S2A–B), and a boxplot of fungal glycogen content (Fig. S3). It also performs the corresponding mixed-effects model and ANOVA analyses.

References

Source studies

Crumière, A. J. J., et al. (2021). The multidimensional nutritional niche of fungus-cultivar provisioning in free-ranging colonies of a neotropical leafcutter ant. Ecology Letters, 24, 2439–2451.

Shik, J. Z., et al. (2021). Nutritional niches reveal fundamental domestication trade-offs in fungus-farming ants. Nature Ecology & Evolution, 5, 122–134.

Software

R Core Team (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Pinheiro, J., et al. (2024). nlme: Linear and nonlinear mixed effects models. R package version 3.1-164.

Lenth, R. V. (2021). emmeans: Estimated marginal means. R package version 1.10.7.

Nychka, D., et al. (2017). fields: Tools for spatial data. R package version 16.3.1.

Kuznetsova, A., Brockhoff, P. B., and Christensen, R. H. B. (2017). lmerTest package: Tests in linear mixed effects models. Journal of Statistical Software, 82, 1–26.

Hijmans, R. J., et al. (2025). terra: Spatial data analysis. R package version 1.8-54.

Borchers, H. W. (2025). pracma: Practical numerical math functions. R package version 2.4.6.

Torchiano, M. (2020). effsize: Efficient effect size computation. R package version 0.8.1.