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

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.

INTRODUCTION

Vertebrates harbor diverse microbial communities in their guts (Ley et al., 2008; Qin et al., 2010) and these so‐called gut microbiota (GM) are involved in many interactions with the host. In addition to its effect on gut function (Jumpertz et al., 2011; Sekirov, Russell, Antunes, & Finlay, 2010), interactions with the host's immune system have important consequences for the host's health and fitness. The different species comprising the GM regulate the host immune system contribute to its development during early ontogenetic stages (Belkaid & Hand, 2014; Kim, Park, & Kim, 2014; Sjögren et al., 2009; Wu & Wu, 2012) and affect the host's capacity to resist invading 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 cell‐mediated response in field ecoimmunology (Møller, 2001; Saino et al., 2002, 2003; Tella, 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.

METHODS

Field data acquisition

Data on fecal microbiota and PHA response were collected during 2014 (late April – late June) from barn swallow nestlings (n = 58) distributed in 32 clutches (Czech Republic, 49° 4′ 7.762″ N, 14° 42′ 36.521″ E, Supporting Information Table S1).

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 (Kreisinger et al., 2017).

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 (v. 13_5, 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 (Kreisinger et al., 2017).

Statistical analysis

Barn swallow GM taxonomic content was visually summarized using Krona tools (Ondov, Bergman, & Phillippy, 2011). All the statistical 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 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.

RESULTS AND DISCUSSION

After all filtering steps, we obtained 947,675 high‐quality reads with a median sequencing depth per sample of 13,798 (range = 1,112–44,777) and an average number of 97% UPARSE OTUs per sample of 153.6 (range = 67–443). In line with our previous study on barn swallow nestlings from the same population (Kreisinger et al., 2017), the most abundant bacterial phyla were Proteobacteria (dominated by Serratia, Pantoea, Providencia, and Diplorickettsia) Firmicutes (dominated by genera Enterococcus, Lactococcus, Lactobacillus, and unassigned Clostridia), Bacteriodetes (dominated by the genus Dysgonomonas), and Actinobacteria. All Actinobacterial genera were represented by low percentage of reads (<1%), with Rhodococcus (0.58% of all read), Rothia (0.56% of reads), and Corynebacterium (0.55% of reads) being the most abundant (Figures 1 and 2).

Figure 1

Bar plots indicating proportions of dominant bacterial phyla and classes in barn swallow fecal microbiota samples

Figure 2

Summary of barn swallow GM taxonomic content. Rare taxa (represented by <2% reads) are labeled as “others”

We did not observe any association between PHA response and GM diversity (LME: p > 0.2 for all alpha diversity measures; Table 1). In theory, a correlation between immune functions and GM diversity might be expected as mutual interactions between immune gene allelic diversity (Bolnick et al., 2014), intensity of immune response (Hawley, Sydenstricker, Kollias, & Dhondt, 2005), parasite load (Kurtz et al., 2004; Madsen & Ujvari, 2006; Sommer, 2005), and overall GM diversity have repeatedly been reported in previous studies. However, consistent with our data, number of previous studies did not find any straightforward correlation between GM richness and immune phenotype (Chang, Hao, Offermanns, & Medzhitov, 2014; Jones et al., 2013; Vatanen et al., 2016).

Table 1

Effect of PHA response and Julian date on GM alpha diversity (assessed as number of observed OTUs, Chao1 predictions of total GM diversity and Shannon index). Shown are LME estimates (Estimate) and corresponding standard errors (SE), deviance changes due to elimination of a given term from the model (χ 2) and associated degrees of freedom (Δdf), and probability values (p)

Response	Explanatory var.	Estimate	SE	χ 2	Δdf	p
Chao1	Intercept	1.897	0.079
	Julian	−0.001	0.001	3.364	1	0.067
	PHA	0.028	0.024	1.388	1	0.239
Observed	Intercept	1.780	0.096
	Julian	−0.001	0.001	0.933	1	0.334
	PHA	0.026	0.030	1.149	1	0.284
Shannon	Intercept	0.600	0.685
	Julian	0.000	0.005	0.002	1	0.968
	PHA	0.152	0.208	0.645	1	0.422

GM: gut microbiota; LME: linear mixed‐effect model; OTUs: operational taxonomic units; PHA: phytohaemagglutinin.

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).

Figure 3

Db‐RDA ordination of GM in barn swallow nestlings. Two dissimilarity types between samples were used as a response (Bray–Curtis and weighted UniFrac), while PHA swelling and Julian date of sampling were included as explanatory variables. Variation along the first two constrained axes is shown. Strength of PHA response (in millimeters) is indicated by color intensity of plotting characters

Table 2

Effect of PHA response and Julian date on GM composition (corresponding to the first db‐RDA axis for weighted UniFrac and Bray–Curtis dissimilarity). Shown are LME estimates (Estimate) and corresponding standard errors (SE), deviance changes due to elimination of a given term from the model (χ 2) and associated degrees of freedom (Δdf), and probability values (p)

Response	Explanatory var.	Estimate	SE	χ 2	Δdf	p
Bray–Curtis	Intercept	−5.053	0.953
	Julian	0.024	0.006	13.479	1	<0.001
	PHA	0.911	0.294	8.724	1	0.003
UniFrac	Intercept	−2.074	1.262
	Julian	0.007	0.008	0.871	1	0.351
	PHA	0.769	0.366	4.523	1	0.033

GM: gut microbiota; LME: linear mixed‐effect model; PHA: phytohaemagglutinin.

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). There are several plausible explanations for the negative correlation between abundances of Enterococcus OTU and PHA response. 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). Unfortunately, we are not able to distinguish between these two alternatives as 16S rRNA region used in our study does not allow reliable species‐level assignment of this particular OTU. Interestingly, association of Enterococcus loads and phenotype was observed in another study on passerine juveniles (González‐Braojos, Vela, Ruiz‐de‐Castañeda, Briones, & Moreno, 2012). In particular, Enterococcus loads were negatively correlated with growth rates, that is, the phenotype trait that can covary with PHA response (Lifjeld, Dunn, & Whittingham, 2002) and other immune parameters as well (van der Most, de Jong, Parmentier, & Verhulst, 2011). Reduced PHA response was also associated higher abundances of OTU from genus Rickettsia, an insect‐borne intracellular pathogen (Parola & Didier, 2001) commonly detected bird GM (Kropáčková, Těšický, et al., 2017).

Table 3

OTUs associated with PHA response according to generalized LMEs for data with negative binomial errors. Shown are taxonomic assignation of individual OTUs, LME‐based estimates, and corresponding standard errors (SE) for OTU abundance versus PHA response association, together with their probability (p) and q values

OTU name	Estimate	SE	p Value	q Value	Class	Family	Genus
OTU_105	−5.7648	1.0746	<0.001	0.0005	Lactobacillales	Streptococcaceae	Lactococcus
OTU_446	−5.5817	1.2421	0.0001	0.0035	Clostridiales	Unassigned	Unassigned
OTU_77	−3.774	0.9505	0.0005	0.0111	Lactobacillales	Enterococcaceae	Enterococcus
OTU_155	−3.5573	0.989	0.0013	0.0207	Pseudomonadales	Pseudomonadaceae	Unassigned
OTU_50	−2.6995	0.7293	0.001	0.0186	Rickettsiales	Rickettsiaceae	Rickettsia
OTU_280	2.7659	0.8043	0.002	0.0218	Lactobacillales	Streptococcaceae	Streptococcus
OTU_45	2.9869	0.8524	0.0017	0.0207	Bacillales	Staphylococcaceae	Staphylococcus
OTU_235	5.7617	1.6456	0.0017	0.0207	Bacillales	Bacillaceae	Bacillus
OTU_104	7.9446	1.4446	<0.001	0.0005	Flavobacteriales	Flavobacteriaceae	Unassigned
OTU_119	24.4801	0.2463	<0.001	<0.001	Bacteroidales	Porphyromonadaceae	Dysgonomonas

LME: linear mixed‐effect model; OTUs: operational taxonomic units; PHA: phytohaemagglutinin.

On the contrary, intensity of PHA response tended to increase with increasing abundance of OTUs from genus Bacillus, Staphylococcus (putatively Staphylococcus saprophyticus; 100% identity according to blastn searches), Dysgonomonas and Streptococcus. Bacteria from genus Bacillus include many common gut symbionts of vertebrates. On the other hand, S. saprophyticus is opportunistic pathogen causing inflammatory diseases of urinary tract in humans (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). However, their effect on physiology and health of avian are still unknown. Many Streptococcus species are vertebrate commensal, but some represent opportunistic pathogens of various vertebrate taxa including birds (Benskin, Wilson, Jones, & Hartley, 2009). Unfortunately, our data did not allow reliable species‐level assignation of this particular OTU.

The contrasting effect of Rickettsia and Staphylococcus OTUs on PHA response suggests that putative bacterial pathogens can be associated both with attenuation and enhancement of PHA response. A similar contradictory pattern has been observed in the case of putative commensal or beneficial bacteria, with Lactococcus OTU abundance, in particular, being negatively related to PHA swelling and Bacillus sp. being positively related. We propose that these seemingly contrasting results are related to both interaction complexity between bacteria and the vertebrate immune system and to the complex immunological background of the PHA swelling response. Further studies targeting specific components of the bird immune system are required, therefore, in order to obtain a better understanding of how (a) host immune system interacts with GM and (b) how the overall pattern of such interactions differ from well‐established mammalian models.

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTIONS

TA and JK designed the study. AP, MA, TA, RM, and OT carried out field sampling. LS, J‐FM, and RM performed laboratory analysis. JK and LS performed data analysis. TA, LS, and JK secured funding. JK drafted the manuscript. All authors provided helpful comments and recommendations and approved the final version of the manuscript.

DATA ACCESSIBILITY

FASTQ files for individual samples are available at European Nucleotide Archive (project accession number: PRJEB27618, http://www.ebi.ac.uk/ena/data/view/PRJEB27618).

Click here for additional data file.

Click here for additional data file.