Fecal microbiota associated with phytohaemagglutinin‐induced immune response in nestlings of a passerine bird

Abstract The vertebrate gastrointestinal tract is inhabited by a diverse community of bacteria, the so‐called gut microbiota (GM). Research on captive mammalian models has revealed tight mutual interactions between immune functions and GM. However, our knowledge of GM versus immune system interactions in wild populations and nonmammalian species remains poor. Here, we focus on the association between GM community structure and immune response measured via the phytohaemagglutinin (PHA) skin swelling test in 12‐day‐old nestlings of a passerine bird, the barn swallow (Hirundo rustica). The PHA test, a widely used method in field ecoimmunology, assesses cell‐mediated immunity. GM structure was inferred based on high‐throughput 16S rRNA sequencing of microbial communities in fecal samples. We did not find any association between PHA response and GM diversity; however, our data revealed that the intensity of PHA response was correlated with differences in GM composition at the whole‐community level. Ten bacterial operational taxonomic units corresponding to both putative commensal and pathogens were identified as drivers of the compositional variation. In conclusion, our study suggests existence of GM versus immune system interactions in a free‐living nonmammalian species, which corresponds with previous research on captive vertebrates.

pathogens (Ivanov et al., 2009). Simultaneously, the host supports a wide range of mechanisms, usually linked with the immune genes, that regulate GM content (Benson et al., 2010;Bolnick et al., 2014).
Most current research focussed on interactions between GM, and the host immune system has used captive-bred animals as a model. However, both taxonomic and functional composition varies considerably between wild and captive populations (Kreisinger, Čížková, Vohánka, & Piálek, 2014;McKenzie et al., 2017). Similarly, both immune parameters and their interindividual variation differ between wild and captive populations due to the altered genetic background of laboratory strains, a lower prevalence of parasites and pathogens and less variation in biotic and abiotic factors involved in immune trait modulation under captive conditions (Boysen, Eide, & Storset, 2011;Flies, Mansfield, Grant, Weldele, & Holekamp, 2015). Hence, results for GM versus immune system interactions obtained from captive populations do not necessarily reflect the selective forces that shape the host's immune system over GM-associated coevolutionary history (Maizels & Nussey, 2013).
Moreover, our knowledge on host immune system versus GM interactions is largely based on mammalian species hosting different GM and having distinct immune system than other vertebrate taxa.
Specifically, bacteria from the Firmicutes and Bacteroidetes phyla typically dominate in the mammalian GM (Ley et al., 2008), whereas nonmammalian vertebrate GM may comprise taxonomically more diverse bacterial consortia (Kropáčková, Těšický, et al., 2017;Sullam et al., 2012). Consequently, further studies dealing with free-living, nonmammalian species are essential for a deeper understanding of the evolutionary forces shaping interactions between GM and the host's immune system.
Here, we study the associations between GM structure and immune response in nestlings of a free-living passerine bird, the barn swallow (Hirundo rustica). The barn swallow is a migratory, insectivorous species with complex social system that breeds in colonies (Cramp & Perrins, 1993;Petrželková et al., 2015). The GM of barn swallows and other birds differs from that of conventional mammalian models (Hird, Carstens, Cardiff, Dittmann, & Brumfield, 2014;Kropáčková, Těšický, et al., 2017;Waite & Taylor, 2014), which makes birds a valuable model group for gaining a deeper insight into GM versus immune system interactions. Various aspects of immune system function have previously been studied in barn swallows and other free-living birds, predominantly related to reproductive behavior and sexual selection (Møller, 2001;Saino, Ambrosini, Martinelli, & Møller, 2002;Saino, Ferrari, Romano, Martinelli, & Møller, 2003).
However, there have been few studies aimed at testing the association between immunity and associated microbial communities (Ruiz-Rodríguez et al., 2009).
We analyzed fecal microbiota profiles using high-throughput sequencing of 16S rRNA amplicons as a proxy for GM. Immune response was assessed via the phytohaemagglutinin (PHA) skin swelling test, which is the most widely used method for assessment of cellmediated response in field ecoimmunology (Møller, 2001;Saino et al., 2002Saino et al., , 2003Tella, Lemus, Carrete, & Blanco, 2008). The PHA assay is traditionally believed to reflect adaptive immune response mediated predominantly by T cells (Goto, Kodama, Okada, & Fujimoto, 1978;Tella et al., 2008). However, recent research suggests that immune mechanisms involved in PHA-induced swelling are more complex, comprising a strong component of innate immunity (Vinkler, Bainová, & Albrecht, 2010;Vinkler, Schnitzer, Munclinger, & Albrecht, 2012). A stronger PHA response is typically interpreted as beneficial due to its positive association with fitness-and condition-related traits (Bowers et al., 2014). Given the complex immunological background of PHA swelling, however, this may not hold universally. There are numerous examples showing no, or even a negative, relationship between PHA responsiveness and body condition, physiological stress, or health status (Møller & Petrie, 2002;Saks, Karu, Ots, & Hõrak, 2006;Vinkler et al., 2012). Despite these complexities, it is worth exploring the potential correlations between GM and PHA responsiveness as PHA-induced swelling is the most widely studied trait in ecoimmunological literature. In addition, previous research supports both positive (Saks et al., 2006) and negative (Merlo, Cutrera, & Zenuto, 2016) association between gut infection by eukaryotic parasites and PHA responsiveness, suggesting that extending such research on prokaryotic communities inhabiting the gut could be potentially fruitful. We therefore combine data on GM profiles with measures of PHA swelling in order to test whether there is any association between GM diversity and immune response in barn swallows. We also assess whether interindividual variation in PHA response is correlated with differences in GM composition and which specific bacterial taxa determine any such variation in PHA response.
Tissue thickness of the left wing web (patagium) of 12-day-old barn swallow nestlings was measured using a standard thickness gauge (Mitutoyo, Japan). Subsequently, the PHA solution (0.10 mg of PHA-P dissolved in 20 μl of DPBS) was injected and the magnitude of the swelling reaction was measured after 24 hr. Both pre-and posttreatment tissue thickness measurements were performed three times by the same person (A.P., accuracy ~0.01 mm). Repeatability of these measurements was high (intraclass correlation coefficient = 0.973 and 0.967 for pre-and posttreatment measurements, respectively). Consequently, the average tissue thickness increment between pre-and posttreatment measurements was used as an index of PHA-induced swelling in subsequent analyses.
Fecal samples of 12-day-old barn swallow nestlings were collected prior PHA injection, placed in sterile cryotubes (Simport, Canada), and stored in liquid nitrogen during field works. After the field works, samples were preserved under −80°C until DNA extractions. Further details on fecal sample collection and storage, together with a description of breeding site, are provided elsewhere .
All field procedures were conducted in accordance with the Guidelines for Animal Care and Treatment of the European Union and approved by the Animal Care and Use Committees of the Czech Academy of Sciences (041/2011) and Charles University (4789/2008-0).

| Microbiome profiling and bioinformatic processing of 16S rRNA data
Metagenomic DNA was isolated from fecal samples using PowerSoil Mo Bio kits (Qiagen). The V3-V4 region of 16S rRNA was amplified using S-D-Bact-0341-b-S-17 (CCTACGGGNGGCWGCAG) and S-D-Bact-0785-a-A-21 (GACTACHVGGGTATCTAATCC) primers (Klindworth et al., 2013), tagged with 10 bp oligonucleotide indices for demultiplexing. Technical PCR duplicates were prepared for all samples in order to check for microbial profile consistency.
Sequencing libraries were prepared using TruSeq Nano Kits and sequenced on Illumina Miseq using v3 chemistry.
The resulting 300 bp long paired-end reads were merged using Pear (Zhang, Kobert, Flouri, & Stamatakis, 2014) and demultiplexed using Mothur (Schloss et al., 2009). Lotus pipeline (Hildebrand, Tadeo, Voigt, Bork, & Raes, 2014) was used for quality filtering (elimination of sequences, if average Q < 30 and if average Q within 50 bp sliding dropped below 25) and elimination of chimeric sequences. Subsequently, UPARSE algorithm (Edgar, 2013) implemented in Lotus was used for clustering of resulting high-quality reads at 97% similarity threshold to operational taxonomic units (OTUs). Taxonomic assignment of representative sequences for each OTU was performed using RDP classifier and Green Genes database DeSantis et al., 2006) as a reference. Representative sequences were aligned using PyNAST (Caporaso et al., 2010) and a phylogenetic tree constructed using FastTree (Price, Dehal, & Arkin, 2010). The OTU table, sample metadata, taxonomic annotations, and phylogenetic tree were stored as a phyloseq object (McMurdie & Holmes, 2013) for further analysis. OTUs not assigned to phylum level, or those classified as chloroplasts (1% and 8.2% reads, respectively), were considered as sequencing artefacts and diet contaminants, respectively, and eliminated from all downstream analyses.
Details on laboratory procedures associated with microbiome profiling and bioinformatic processing of sequencing data were provided in a previous study on this species ).

| Statistical analysis
Barn swallow GM taxonomic content was visually summarized using Krona tools (Ondov, Bergman, & Phillippy, 2011). All the statistical F I G U R E 1 Bar plots indicating proportions of dominant bacterial phyla and classes in barn swallow fecal microbiota samples analyses were conducted using packages running under R 3.4.3 software (R Core Team, 2016). As we detected significant correlation between PHA response and Julian date of fecal sample collection (Pearson r = 0.324, p < 0.01), we controlled all subsequent statistical analyses for effect of sampling date. Association between microbiota diversity (i.e., number of observed OTUs, Chao1 diversity estimates and Shannon diversities) and PHA response or Julian date of sample collection was tested using linear mixed-effect models (LME, R package lme4; Bates, Mächler, Bolker, & Walker, 2015) with Gaussian distribution of errors. Nest identity was included as a random intercept.
Next, weighted UniFrac (Lozupone & Knight, 2005) and Bray-Curtis community dissimilarity between samples were calculated based on sample-specific OTU proportions. The effect of PHA-induced response and Julian date was assessed using distance-based redundancy analysis (db-RDA, Legendre & Anderson, 1999), with the matrix of between-sample dissimilarities included as a response.
Permutation-based ANOVA (anova.cca function from R package vegan, Oksanen et al., 2013) was then used to test for significance of the constrained db-RDA axes. According to these analyses, only the first constrained db-RDA axis was significant (p < 0.001 for both UniFrac and Bray-Curtis dissimilarity), and the effect of the second constrained db-RDA axis was nonsignificant (p > 0.3 in both cases).
We then extracted the scores for the first db-RDA axis and tested whether they were significantly associated with PHA response and/ or Julian date using LME. We argue that this analysis method is preferable to the default anova.cca, as this function cannot effectively F I G U R E 2 Summary of barn swallow GM taxonomic content. Rare taxa (represented by <2% reads) are labeled as "others" account for pseudoreplications induced by sampling multiple individuals from the same nest.
Associations between abundance of OTUs and PHA response were tested using generalized LMEs from R package BhGLM for data with negative binomial distribution of errors (Zhang et al., 2017).
OTU-specific read counts within individual samples were included as a response, while Julian date of sample collection and PHA response were included as explanatory variables. Log-transformed total number of reads per sample was specified as model offset (i.e., assuming number of reads per given OTU to be proportional to total number of reads per individual sample) and clutch identity as random effect.
The qvalue method (Storey & Tibshirani, 2003) was subsequently used to account for false discoveries due to multiple testing. To optimize sensitivity of OTU-level analyses, we applied "independent filtering" procedure (Bourgon, Gentleman, & Huber, 2010) using DESeq2 R package (Love, Huber, & Anders, 2014) and considered only OTUs that passed this step (n = 196 OTUs, representing 96% of all high-quality reads). Procrustean analysis revealed high congruence between original and subsetted microbial profiles (Procrustean r = 0.9974, Procrustean sum of squares = 0.0052, p < 0.0001), suggesting that resulting OTU subset covered representative variation in the GM content.  (Figures 1 and 2).

| RE SULTS AND D ISCUSS I ON
We did not observe any association between PHA response and GM diversity (LME: p > 0.2 for all alpha diversity measures; Table 1).
Despite the lack of any relationship between GM diversity and PHA response, we observed significant correlation between magnitude of PHA swelling and variation in GM composition at the whole-community level, suggesting that individuals with similar GM composition had a similar PHA response. This association was specifically implied based on db-RDA ordination (Figure 3). In addition, LMEs running on scores for the first db-RDA axis (using both Bray-Curtis and weighted UniFrac distance) revealed a significant effect of PHA response after statistical control for Julian date of sample collection, while the effect of Julian date itself was significant only in the case of db-RDA for Bray-Curtis dissimilarity (Table 2).
According to OTU-centered negative binomial LMEs, ten OTUs represented by relatively low number of reads (~2.2% in total) exhibited significant association with the intensity of PHA response (Table 3 and Supporting Information File S1). Two of these OTUs, belonging to Lactic Acid Bacteria from genus Enterococcus and Lactococcus, were negatively related to PHA swelling. Intensity of PHA-induced swelling seems to strongly reflect general proinflammatory potential of given individual (Vinkler et al., 2010).
Consequently, observed negative correlation between Lactococcus abundances and PHA response can be related to anti-inflammatory effect that was previously described for some probiotic species from this genus (Han, Lee, Park, & Paik, 2015;Luerce et al., 2014). Similarly, as in the case of Lactococcus, some Enterococcus species exhibit probiotic properties. However, Enterococcus genus includes also several pathogenic strains, whose infection can directly affect host's immunity (Fisher & Phillips, 2009 (Hovelius & Mårdh, 1984) and Dysgonomonas can cause gut inflammation in immunocompromised human subjects (Bangsborg, Frederiksen, & Bruun, 1990). Both these OTUs were previously detected in bird GM (Kropáčková, Pechmanová, et al., 2017;Kropáčková, Těšický, et al., 2017;Xenoulis et al., 2010). We propose that these seemingly contrasting results are related to both interaction complexity between bacteria and the verte-

CO N FLI C T O F I NTE R E S T
None declared.