Data analysis
We used the proportion of time spent on each of the saplings for the analysis. We used the glmmTMB package (Brooks et al. 2017) within R 3.6.1 (R Core Team 2019) to build a generalised linear mixed model with a beta error structure and a logit link. Models that follow beta distribution only allow proportion values between 0 and 1. Therefore, 0 and 1 values in our study were converted into 0.0001 and 0.9999, respectively. We used the mixed model to determine the effect of experiment (sapling combination), bird training and their interaction on the proportion of time spent by an individual bird in proximity of each sapling. We used Anova function from ‘car’ package to calculate P-values for each variable using Wald chisquare test (Fox & Weisberg 2018). We then made pairwise comparisons between the five experiments within each of the training types using the emmeans functions in ‘emmeans’ package that adjusts P-values (Padj) following the Tukey method (Lenth 2007).
The centroided GC‐MS data in NetCDF format were further processed using XCMS online (version 2.7.2) The peaks in each sample were detected with the cent Wave algorithm (peak width 2–20 s; signal to noise threshold 10). Peak grouping across samples was restricted to peaks present in at least 50% of the samples in at least one treatment group (minfrac = 0.5). Retention time correction was accomplished with the symmetric method and nonlinear loess‐smoothing and iterated three times with decreasing bandwidth parameter for the grouping from 10 to 0.2 s. The extracted ion species were grouped according to their parent molecule into pseudospectra with the Bioconductor package camera (Kuhlet al. 2012). This resulted in a final feature table containing mass‐to‐charge ratio (m/z), retention time, peak area and the pseudospectra group for each detected ion species. Groups containing only a single ion species were considered to be artefacts and removed from the final feature table. The feature table was uploaded to MetaboAnalyst for further statistical analyses. We performed quantile normalization and Log transformation. Pareto Scaling was chosen for Principal Component Analysis. The annotation of the most significant features was done by spectral library search (NIST) and MS spectra and Kovats Index comparison with standard compounds.