Study organism and bacterial sampling
In order to investigate the skin microbiota of wild zebra finches (Taeniopygia guttata), individuals from a breeding colony around the UNSW Arid Zone Research Station at Fowlers Gap, Australia (31°05’13.1’’ S 141°42’17.4’’ E), details, see: [25, 36], were sampled during September and October 2016. Study animals comprised twelve zebra finch families, consisting of three to four individuals (n = 40; 9 females, 8 males, 23 nestlings). Families were defined as two parents and two nestlings (family of four individuals, n = 4) or two parents and one nestling (family of three individuals, n = 1) or one parent with two nestlings (family of three individuals, n = 7). Samples that did not fit in one of these groups were excluded. Parents were always a male and a female of unknown age, the offspring was between six and ten days old at the day of sampling. Offspring sex was not identifiable at this developmental stage.
For the sampling, adults were caught with nest box traps. When a bird was caught, it was retrieved from the nest box and sampled immediately, wearing a fresh pair of sterile gloves for every bird. Adults were then held in bird bags until catching was completed and afterwards the chicks were sampled, taking each out of the nest and sampling it separately, wearing a fresh pair of gloves. All birds were sampled with a nylon flocked swab (ESwab, Copan Italia, Italy) on the bare skin parts around the preen gland. In order to do so, a swab was stroked on the skin adjacent to the preen gland, sparing the skin of the gland itself to avoid taking up preen gland secretions (for a detailed description see Engel et al. 2018). As our aim was to sample skin bacteria, the feathers of adult birds were parted and the swab rubbed on the skin, trying to avoid feather contact. Nevertheless, touching the down feathers could not be avoided completely. Samples from every family were taken at the same day, whenever possible. Swabs were transferred to liquid Amies medium immediately and transported on ice in a cooler box until storage at -20 °C within 5 hours.
Shipping of the frozen samples occurred in February 2017 on dry ice from Australia. In Germany, DNA of all samples was extracted within one week after arrival using the BiOstic Bacteremia DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA, USA) according to the manufacturer's instructions with minor adjustments. The entire volume of 800 µl liquid Amies medium containing floating bacteria was used (instead of a sample of 1800 µl blood as instructed in the kit’s manual). Furthermore, samples in bead tubes were shaken vertically for 10 min at 50 Hz in a TissueLyser instrument (TissueLyser LT, Qiagen, Venlo, Netherlands).
Library preparation and DNA sequencing
Library preparation was performed following the Illumina 16s Metagenomic Library Prep Guide, 15044223-b and the detailed description provided by Jervis-Bardy et al. . We used a slightly modified protocol to accommodate our need for low DNA concentrations. The primer pair we used in our study is derived from Klindworth et al.  and targeted the hypervariable regions V3 and V4 of the 16S rRNA gene (S-D-Bact-0341-b-S-17: 5’-CCTACGGGNGGCWGCAG-3’ and S-D-Bact-0785-a-A-21: 5’-GACTACHVGGGTATCTAATCC-3’). Primers were connected to Illumina overhang adapter sequences as follows:
forward: 5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG -3’,
For a detailed description see Supporting information S1 in .
The libraries were sequenced on an Illumina MiSeq system (Illumina, Inc., San Diego, CA, USA) in paired-end mode (2 x 300 sequencing cycles) at the CeBiTec, Bielefeld University, using the MiSeq Reagent Kit v3 (Illumina, Inc., San Diego, CA, USA) following the manufacturer’s manual.
Bioinformatic data processing
Data was processed analogous to our previous study , with some minor adjustments, as described briefly in the following. MiSeq paired-end reads were assembled with Flash v1.2.11  in an iterative manner coupled with read clipping based on average quality score values using sickle v1.33  in order to achieve a maximum assembly rate. Reads initially failing the read assembly were clipped to a q20 quality threshold and tried to assemble again. This process was repeated with a repeatedly increased quality threshold by 3 up to the limit of q35. Assembled PE reads were then screened for matching forward and reverse amplification primers, in both directions consecutively, using cutadapt v2.1, and matched primers were trimmed during this procedure. Reads matching the primers in reverse direction were reverse complemented hereby. Then, mothur v.1.41.3  was used to de-replicate, align, filter, and de-noise reads in accordance to the mothur MiSeq standard operating procedure. Read sets were consecutively screened for chimeras and clustered into OTUs with a percent identity threshold of 97% using the USEARCH v8.0.147 tool suite . Finally, each OTU was taxonomically labelled based on naïve Bayesian read classification in mothur using the full SILVA database v132  as a reference and a confidence cut-off of 80%. Please see Engel et al. 2008 for more details on the bioinformatics processing, invoked commands and used parameters.
Data filtering and statistical analyses
All statistical analyses were performed using R version 3.3.3.  and Primer-e software version 7 . After bioinformatic processing, we applied additional filtering steps to the OTU table as described in Pankoke et al. . In short, only bacterial OTUs that could be classified at phylum level were used for further filtering steps. OTUs with less than 10 read counts per cluster as well as OTUs that did not occur in at least two samples were removed from the dataset. After filtering, 1378 OTUs with a mean read count of 24194 ± 10885 were used for a rarefaction analysis. Samples were rarefied to 5536 read counts, using rrarefy() from the vegan R package  resulted in 1233 OTUs.
Using the respective function implemented in vegan, we estimated bacterial diversity by calculating the Shannon diversity index (alpha diversity). To overcome limitations of unbalanced sampling, i.e. not always two parents and two siblings per family, we compared diversity scores with linear models on reduced datasets. To estimate the impact of bird age on microbial diversity, we selected the fully replicated families consisting of two adults (female, male) and two nestlings (offspring) (n = 4). To estimate the influence of sex on microbial diversity, we analysed the paired adults from those n = 5 families that consisted of two adults and at least one nestling. And to estimate the influence of family on diversity scores, we analysed all n = 12 families that consisted of three to four individuals (one or two adults and one or two nestlings). Residuals of the models were visually inspected for normality and deviation from homoscedasticity .
The rarefied OTU abundance data were transformed (log10(x + 1)) prior to statistical analysis. For all further analyses investigating beta diversity, we computed a pairwise similarity matrix based on the Bray-Curtis similarity index and analysed potential differences between a priori defined groups (here: age, sex and families) using a non-parametric analysis of similarities (ANOSIM). We performed the ANOSIM using Primer 7 . To visualise the similarities or differences in microbial communities, we used a non-metric multidimensional scaling plot (nMDS).
Moreover, we explored whether the amount of shared skin bacterial communities differed between members within a family. Therefore, we calculated pairwise similarity indices (based on Bray-Curtis similarity) for the following pairs: parent-parent, parent-offspring and offspring-offspring using the rarefied OTU table (using Primer 7). To test whether similarities differed among these groups we used a Kruskal-Wallis test, followed by post-hoc Wilcoxon rank sum tests for comparisons between groups (using R).
To test the hypothesis that spatial or environmental factors influenced the skin microbiome of the zebra finch, the rarefied OTU table was first converted to a distance matrix using the Bray-Curtis distance. Next, we calculated a distance matrix of the spatial coordinates of the nest boxes using distm(coordinates (), fun = distHaversine) from the packages geosphere  and sp . Finally, we performed a Mantel test as implemented in vegan to analyse whether the distance matrix of the rarefied OTU table and the distance matrix of the spatial coordinates were correlated. The matrices used for the analyses are attached as additional file 1.
Taxonomic profiling of the skin microbiota of the zebra finch was done according to Pankoke et al. . In short, read counts of all OTUs belonging to one specific taxon were summed up and normalized by the total number of read counts, yielding proportions. Each taxonomic level was then deduced from the normalized OTU counts by summing up the taxonomic classifications on different levels. The respective proportions of the different taxa were ranked by decreasing value and plotted as stack bars. All taxa with less than 0.5% of prevalence were summed up and were being classified as of “low prevalence”.