Overview

This post summarizes fungal community data from carrion-associated insects collected during decomposition experiments. Using ITS1 amplicon sequencing, we characterized the mycobiomes of two key decomposer species: the carrion beetle Necrodes surinamensis and the secondary screwworm fly Cochliomyia macellaria.

Key findings:

  • Distinct fungal communities between beetle and fly species
  • Body-region specific differences within beetles (abdomen vs head-thorax)
  • Identification of decomposer-associated fungal genera
  • Differential abundance patterns suggesting species-specific fungal associations

Table of Contents

Methods Summary Quality & Processing Alpha Diversity Beta Diversity Taxonomic Composition Differential Abundance Decomposer Fungi Interpretation Data Access Next Steps References

Methods Summary

Samples:

  • Necrodes surinamensis (n=39): Abdomen and head-thorax samples
  • Cochliomyia macellaria (n=19): Whole specimens

Sequencing: ITS1 amplicon sequencing (fungi-specific), rarefied to 500 reads per sample

Analyses:

  • Alpha diversity (Shannon, Observed richness, etc.)
  • Beta diversity (Bray-Curtis and Jaccard PCoA)
  • Differential abundance (DESeq2)
  • Prevalence analysis of targeted decomposer fungi

QIIME2 Processing Pipeline

Click to view full QIIME2 workflow

1. Import EMP Format Data

# Activate QIIME2 environment
conda activate qiime2-amplicon-2025.10

# Import Run 1 (EMP paired-end format)
qiime tools import \
  --type EMPPairedEndSequences \
  --input-path Run1 \
  --input-format EMPPairedEndDirFmt \
  --output-path emp_paired_end_run1.qza

# Import Run 2
qiime tools import \
  --type EMPPairedEndSequences \
  --input-path Run2 \
  --input-format EMPPairedEndDirFmt \
  --output-path emp_paired_end_run2.qza

2. Demultiplex by Barcodes

# Demultiplex Run 1 using EMP barcodes
qiime demux emp-paired \
  --i-seqs emp_paired_end_run1.qza \
  --m-barcodes-file demux.tsv \
  --m-barcodes-column BarcodeSequence \
  --p-no-golay-error-correction \
  --o-per-sample-sequences demux_paired_run1.qza \
  --o-error-correction-details error_correction_run1.qza \
  --verbose

# Demultiplex Run 2
qiime demux emp-paired \
  --i-seqs emp_paired_end_run2.qza \
  --m-barcodes-file demux.tsv \
  --m-barcodes-column BarcodeSequence \
  --p-no-golay-error-correction \
  --o-per-sample-sequences demux_paired_run2.qza \
  --o-error-correction-details error_correction_run2.qza \
  --verbose

# Visualize demux quality for both runs
qiime demux summarize \
  --i-data demux_paired_run1.qza \
  --o-visualization demux_summary_run1.qzv

qiime demux summarize \
  --i-data demux_paired_run2.qza \
  --o-visualization demux_summary_run2.qzv

3. Trim ITS1 Primers with Cutadapt

# ITS1F: CTTGGTCATTTAGAGGAAGTAA
# ITS2: GCTGCGTTCTTCATCGATGC

# Trim Run 1
qiime cutadapt trim-paired \
  --i-demultiplexed-sequences demux_paired_run1.qza \
  --p-front-f CTTGGTCATTTAGAGGAAGTAA \
  --p-front-r GCTGCGTTCTTCATCGATGC \
  --p-adapter-f GCATCGATGAAGAACGCAGC \
  --p-adapter-r TTACTTCCTCTAAATGACCAAG \
  --p-discard-untrimmed \
  --p-minimum-length 50 \
  --p-cores 0 \
  --o-trimmed-sequences trimmed_seqs_run1.qza

# Trim Run 2
qiime cutadapt trim-paired \
  --i-demultiplexed-sequences demux_paired_run2.qza \
  --p-front-f CTTGGTCATTTAGAGGAAGTAA \
  --p-front-r GCTGCGTTCTTCATCGATGC \
  --p-adapter-f GCATCGATGAAGAACGCAGC \
  --p-adapter-r TTACTTCCTCTAAATGACCAAG \
  --p-discard-untrimmed \
  --p-minimum-length 50 \
  --p-cores 0 \
  --o-trimmed-sequences trimmed_seqs_run2.qza

# Verify primer removal
qiime demux summarize \
  --i-data trimmed_seqs_run1.qza \
  --o-visualization trimmed_summary_run1.qzv

qiime demux summarize \
  --i-data trimmed_seqs_run2.qza \
  --o-visualization trimmed_summary_run2.qzv

4. ITSxpress Attempt (Not Used)

# NOTE: ITSxpress was tested to extract ITS1 region precisely
# using HMM profiles, but removed too much sequence and results
# were unusable. Proceeded with primer-trimmed sequences instead.

# qiime q2-itsxpress trim-pair-output-unmerged \
#   --i-per-sample-sequences trimmed_seqs_run1.qza \
#   --p-region ITS1 \
#   --p-taxa ALL \
#   --o-trimmed its1_seqs_run1.qza \
#   --verbose

5. Denoise with DADA2

# DADA2 on primer-trimmed sequences
# No length truncation (trunc-len-f/r = 0) for ITS variable length
# max-ee = 2 for quality filtering
# min-overlap = 12 bp for merging paired ends

# Run 1
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs trimmed_seqs_run1.qza \
  --p-trunc-len-f 0 \
  --p-trunc-len-r 0 \
  --p-max-ee-f 2 \
  --p-max-ee-r 2 \
  --p-min-overlap 12 \
  --p-chimera-method consensus \
  --p-n-threads 0 \
  --o-table table_run1.qza \
  --o-representative-sequences rep_seqs_run1.qza \
  --o-denoising-stats denoising_stats_run1.qza \
  --verbose

# Run 2
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs trimmed_seqs_run2.qza \
  --p-trunc-len-f 0 \
  --p-trunc-len-r 0 \
  --p-max-ee-f 2 \
  --p-max-ee-r 2 \
  --p-min-overlap 12 \
  --p-chimera-method consensus \
  --p-n-threads 0 \
  --o-table table_run2.qza \
  --o-representative-sequences rep_seqs_run2.qza \
  --o-denoising-stats denoising_stats_run2.qza \
  --verbose

# Visualize denoising stats
qiime metadata tabulate \
  --m-input-file denoising_stats_run1.qza \
  --o-visualization denoising_stats_run1.qzv

qiime metadata tabulate \
  --m-input-file denoising_stats_run2.qza \
  --o-visualization denoising_stats_run2.qzv

6. Summarize Individual Runs

qiime feature-table summarize \
  --i-table table_run1.qza \
  --m-sample-metadata-file metadata_its.tsv \
  --o-visualization table_summary_run1.qzv

qiime feature-table summarize \
  --i-table table_run2.qza \
  --m-sample-metadata-file metadata_its.tsv \
  --o-visualization table_summary_run2.qzv

7. Merge Sequencing Runs

# Merge feature tables (sum counts for samples present in both runs)
qiime feature-table merge \
  --i-tables table_run1.qza \
  --i-tables table_run2.qza \
  --p-overlap-method sum \
  --o-merged-table table_merged.qza

# Merge representative sequences
qiime feature-table merge-seqs \
  --i-data rep_seqs_run1.qza \
  --i-data rep_seqs_run2.qza \
  --o-merged-data rep_seqs_merged.qza

# Summarize merged table
qiime feature-table summarize \
  --i-table table_merged.qza \
  --m-sample-metadata-file metadata_its.tsv \
  --o-visualization table_merged_summary.qzv

8. Taxonomy Assignment

# Classify ASVs using pre-trained UNITE classifier
qiime feature-classifier classify-sklearn \
  --i-classifier unite_ver2025-02-19_dynamic_s_fungi-Q2-2026.1.qza \
  --i-reads rep_seqs_merged.qza \
  --o-classification taxonomy.qza

# Visualize taxonomy assignments
qiime metadata tabulate \
  --m-input-file taxonomy.qza \
  --o-visualization taxonomy.qzv

# Create taxonomy barplots
qiime taxa barplot \
  --i-table table_merged.qza \
  --i-taxonomy taxonomy.qza \
  --m-metadata-file metadata_its.tsv \
  --o-visualization taxa_barplot.qzv

9. Filter to Fungi Only

# Remove non-fungal sequences (e.g., plant/insect host ITS)
qiime taxa filter-table \
  --i-table table_merged.qza \
  --i-taxonomy taxonomy.qza \
  --p-include k__Fungi \
  --p-exclude Unassigned \
  --o-filtered-table table_fungi_only.qza

# Final summary used for downstream analysis in R/phyloseq
qiime feature-table summarize \
  --i-table table_fungi_only.qza \
  --m-sample-metadata-file metadata_its.tsv \
  --o-visualization table_fungi_only_summary.qzv

10. Export for R/phyloseq Analysis

# Export feature table
qiime tools export \
  --input-path table_fungi_only.qza \
  --output-path exported/

# Export taxonomy
qiime tools export \
  --input-path taxonomy.qza \
  --output-path exported/

# Convert BIOM to text format (for manual inspection if needed)
biom convert \
  -i exported/feature-table.biom \
  -o exported/feature-table.tsv \
  --to-tsv

# The .qza files were then imported directly into R using qiime2R::qza_to_phyloseq()
# See Scripts/analysis.R for downstream statistical analyses

Software Versions

  • QIIME2: 2025.10 (qiime2-amplicon-2025.10)
  • DADA2: 1.32
  • Cutadapt: 4.9
  • ITSxpress: 2.1.4 (tested but not used)
  • Classifier: UNITE ver2025-02-19 (dynamic species hypothesis)
  • Rarefaction depth: 500 reads/sample

Reference: QIIME2 ITS Tutorial

Sequencing Quality & Data Processing

Before analyzing the mycobiome, I assessed sequencing quality across the two independent ITS1 sequencing runs that were merged for this analysis.

Quality Score Profiles

Figure 1: Median Phred quality scores across read positions for both sequencing runs. Forward (blue) and reverse (orange) reads shown separately. Gray lines indicate Q20 and Q30 thresholds.

Read Depth and DADA2 Processing

Figure 2: Sequencing depth and DADA2 processing, faceted by sequencing run. Top: Raw read counts per sample with rarefaction threshold (red dashed line at 500 reads). Bottom: Read retention through DADA2 pipeline showing filtering, merging, and chimera removal steps. Both plots show Run1 and Run2 side by side for comparison.

Key observations:

  • Wide variation in sequencing depth (500-9,000 reads per sample)
  • Most samples exceeded the 500-read rarefaction threshold
  • DADA2 filtering removed ~10-40% of reads based on quality scores
  • Merging step was efficient (~10% loss), reflecting good overlap between paired reads

ASV Diversity Patterns

After denoising and filtering, the final dataset contained 397 unique ASVs across all samples. ASV richness (number of unique sequences detected) varied substantially between samples and species.

ASV Richness per Sample

Figure 3: Number of unique ASVs detected per sample, colored by species. Samples ordered by total read count.

Rank-Abundance Curves

Figure 4: Rank-abundance curves showing the distribution of reads across ASVs. Each line represents one sample, colored by species. Steep curves indicate dominance by few ASVs; flatter curves indicate more even communities.

Feature Table

Click to view interactive ASV × Sample count table

Loading feature table...

Table shows read counts for each ASV (rows) across all samples (columns). Use your browser's search function (Cmd/Ctrl+F) to find specific ASVs or samples.

Alpha Diversity

Figure 5: Alpha diversity metrics by insect species (2x2 grid). Shannon entropy, observed ASV richness, Chao1 estimated richness, and Pielou's evenness all show higher diversity in beetles compared to flies.

Diversity Metric By Insect Species (df = 1) By Body Region (df = 2)
Shannon diversity χ² = 4.195, p = 0.041* χ² = 4.330, p = 0.115
Observed richness χ² = 0.138, p = 0.710 χ² = 0.728, p = 0.695
Pielou’s evenness χ² = 13.167, p < 0.001* χ² = 13.360, p = 0.001*
Chao1 richness χ² = 0.280, p = 0.597 χ² = 0.869, p = 0.648

Shannon diversity and Pielou’s evenness both differ significantly between beetle and fly species, with beetles showing higher diversity and more even communities. Observed and Chao1 richness metrics show no significant differences. Evenness also varies significantly by body region, suggesting compositional differences between abdomen and head-thorax samples in beetles.

Beta Diversity

Figure 6: Bray-Curtis PCoA colored by insect species. Circles = whole specimens (flies), triangles = abdomen, squares = head-thorax.

PERMANOVA results (999 permutations):

Comparison Distance Metric df R² F-statistic p-value
Insect Species Bray-Curtis 1 0.122 6.112 < 0.001*
Insect Species Jaccard 1 0.088 4.244 < 0.001*
Body Region Bray-Curtis 2 0.139 3.484 < 0.001*

Both Bray-Curtis and Jaccard dissimilarities show significant separation between beetle and fly mycobiomes. Body region (whole specimen vs. abdomen vs. head-thorax) also explains significant variation in fungal community composition, with the strongest effect observed for insect species identity.

Taxonomic Composition

Genus-Level Relative Abundance

Figure 7: Genus-level relative abundance across all samples. Stacked bars show fungal community composition, with samples grouped by species.

The heatmap below shows log10-transformed relative abundances for the 30 most abundant ASVs across all samples:

Figure 8: Heatmap of top 30 ASVs (log10 relative abundance). Samples ordered by species and body region. Color intensity reflects abundance from absent (yellow) to highly abundant (red).

Differential Abundance

DESeq2 identified significantly enriched fungal taxa in each species:

Figure 9: Volcano plot of differential abundance. Points above the dashed line are FDR < 0.05, beyond vertical lines are |log2FC| > 1.

Decomposer-Associated Fungi

When targeting specific genera known from decomposition and carrion ecology, prevalence varies substantially between species:

Figure 10: Prevalence of target decomposer fungi (>0.1% relative abundance).

The stacked bar charts below show relative abundances of curated decomposer-associated fungal genera across samples:

Figure 11a: Relative abundance of curated fungal genera in C. macellaria (fly) samples.

Figure 11b: Relative abundance of curated fungal genera in N. surinamensis (beetle) samples, split by body part (Abdomen vs Head-Thorax).

Discussion

  1. Species-specific mycobiomes where Beetles and flies harbor fundamentally different fungal communities, likely reflecting dietary ecology (beetles feed on carrion tissue, flies are necrophagous as larvae but not adults in our sampling).

  2. Decomposer-specialists in beetles are reflected by a high prevalence of yeasts (Rhodotorula, Yarrowia) and stinkhorn fungi (Phallus, Clathrus) in beetles aligns with their role as carrion specialists.

  3. Fungi in flies like Cladosporium are common environmental molds, possibly reflecting transient exposure rather than stable associations.

  4. Body region effects, where within-beetle differences (abdomen vs head-thorax) are modest but detectable, potentially due to microhabitat differences or organ-specific microbial niches.

References

Chakraborty, A., Modlinger, R., Ashraf, M. Z., Synek, J., Schlyter, F., & Roy, A. (2020). Core Mycobiome and Their Ecological Relevance in the Gut of Five Ips Bark Beetles (Coleoptera: Curculionidae: Scolytinae). Frontiers in Microbiology, 11. https://doi.org/10.3389/fmicb.2020.568853

Kudo, R., Masuya, H., Endoh, R., Kikuchi, T., & Ikeda, H. (2018). Gut bacterial and fungal communities in ground-dwelling beetles are associated with host food habit and habitat. The ISME Journal, 13(3), 676–685. https://doi.org/10.1038/s41396-018-0298-3

Lauber, C. L., Metcalf, J. L., Keepers, K., Ackermann, G., Carter, D. O., & Knight, R. (2014). Vertebrate Decomposition Is Accelerated by Soil Microbes. Applied and Environmental Microbiology, 80(16), 4920–4929. https://doi.org/10.1128/aem.00957-14

Shukla, S. P., Plata, C., Reichelt, M., Steiger, S., Heckel, D. G., Kaltenpoth, M., Vilcinskas, A., & Vogel, H. (2018). Microbiome-assisted carrion preservation aids larval development in a burying beetle. Proceedings of the National Academy of Sciences, 115(44), 11274–11279. https://doi.org/10.1073/pnas.1812808115

Valentini, B., Penna, M., Viazzo, M., Caprio, E., Casacci, L. P., Barbero, F., & Stefanini, I. (2024). Yeasts, arthropods, and environmental matrix: a triad to disentangle the multi-level definition of biodiversity. Scientific Reports, 14(1). https://doi.org/10.1038/s41598-024-70327-4

This is a data snapshot from 2026-02-23. Analysis methods and interpretations are subject to refinement as the project progresses.