Data processing

Illumina sequencing data were processed using DADA2 version 1.12.1 (ref. [58]) to produce an error-corrected table of amplicon sequence variant (ASV) abundances in each sample. We used default parameters for all DADA2 analyses except where specified below. Sequences were trimmed to remove forward and reverse primers and to a length of 201 nucleotides (forward reads) and 264 nucleotides (reverse reads) before ASV identification and merging of paired-end reads to a single ASV sequence. We assigned a taxonomy for all ASVs using the SILVA database [59] and removed all non-bacterial sequences (e.g., chloroplasts, eukaryotes; 1% of sequences) prior to further data processing. We used the decontam R package [60] to remove probable contaminant ASVs that were more abundant in negative controls than in experimental samples. We performed the decontamination on two data subsets: (1) blue tit feces samples and (2) leaf and caterpillar samples (which came from the same sequencing run) and in two steps: (1) using the PCR negative controls and (2) using the extraction and field negative controls as references to identify contaminants. The decontaminated ASV table was rarefied to 4000 reads per sample using the R package phyloseq [61]. We used a rarefaction threshold of 4000 reads to reach a plateau of ASV richness per sample while allowing the inclusion of 70% of all samples. All negative controls were excluded at the rarefaction step because they contained fewer than 4000 sequences. All the positive controls had a similar composition and contained the ASVs expected from the known composition of the mock community. The final dataset contained 37 leaf, 118 caterpillar, and 266 blue tit bacterial microbiota samples from 77 males, 64 females, and 125 nestlings. For details about rarefaction curves and control bacterial taxonomic composition, see the supplementary information (Supplementary Figs. 1 and 2).