Within-community variation of interspecific divergence patterns in passerine gut microbiota

Abstract

Gut microbiota (GM) often exhibit variation between different host species and co-divergence with hosts' phylogeny. Identifying these patterns is a key for understanding the mechanisms that shaped symbiosis between GM and its hosts. Therefore, both GM-host species specificity and GM-host co-divergence have been investigated by numerous studies. However, most of them neglected a possibility that different groups of bacteria within GM can vary in the tightness of their association with the host. Consequently, unlike most of these studies, we aimed to directly address how the strength of GM-host species specificity and GM-host co-divergence vary across different GM clades. We decomposed GM communities of 52 passerine species (394 individuals), characterized by 16S rRNA amplicon sequence variant (ASV) profiles, into monophyletic Binned Taxonomic units (BTUs). Subsequently, we analyzed strength of host species specificity and correlation with host phylogeny separately for resulting BTUs. We found that most BTUs exhibited significant host-species specificity in their composition. Notably, BTUs exhibiting high host-species specificity comprised bacterial taxa known to impact host's physiology and immune system. However, BTUs rarely displayed significant co-divergence with host phylogeny, suggesting that passerine GM evolution is not shaped primarily through a shared evolutionary history between the host and its gut microbes.

Keywords: 16S rRNA; co‐divergence; gut microbiota; metabarcoding; passerines.

FIGURE 1

Analytical approach used in our study. 1. GM profiling: Raw paired‐end reads were merged and sequencing errors were corrected by DADA2 to obtain amplicon sequence variants (ASVs). 2. BTUs definition: ASVs were grouped to BTUs using three methods in parallel, namely a) reference based taxonomic assignment of ASVs to genus level and reference free methods defining BTUs sharing at least b) 95% or c) 97% ASVs sequence similarity, based on distance matrix. 3. Quantification of GM‐host species specificity and co‐divergence at the BTU level: Several tests for different aspects of GM diversity were conducted on ASVs distribution within individual BTUs, namely a) GM‐host species specificity test based on permutational MANOVA of inter‐sample GM composition distances (Bray–Curtis and Jaccard), b) compositional GM‐host co‐divergence test based on PACo analyses of inter‐sample GM composition distances and host phylogeny, c) phylogenetical GM‐host co‐divergence test based on PACo analyses of ASVs phylogeny within individual BTUs and hosts phylogeny

FIGURE 2

Correlations of test effect sizes in examined BTUs. Correlation of GM‐host species specificity versus compositional and phylogenetic GM‐host co‐divergence among BTUs (represented by points). Only a subset BTUs passing filters for testing are shown. Kendall's tau was used as the correlation coefficient. The blue line and shaded area correspond to predictions of linear regression and 95% confidence intervals, respectively. Only results for Jaccard distances are shown

FIGURE 3

Heatmap of measured effect sizes in tested reference‐based BTUs. GM‐host species specificity (R ²adj) and compositional GM‐host co‐divergence (P_COMP) are based on Jaccard distances. * indicates significant result (FDR < 0.05), the dendrogram on the left represents phylogeny of the bacterial genera, that is, a phylogenetic tree based on 16S rRNA sequences extracted from GreenGenes database (97% OTUs tree version 12_10) formatted for the purpose of QIIME1)

FIGURE 4

Whole GM community classified according to outcomes of BTU‐level analyses. The figure depicts GM fractions that were not included in BTU‐level analyses, where corresponding BTU level analysis was nonsignificant, significant or significant with high effect‐size (R ²adj/P_COMP/P_PHYLO > 0.2). GM fractions are expressed as (a) total number of ASVs or (b) proportion of all reads. Note that, while modest number of ASVs were actually tested, they represented the majority of all reads in the dataset

FIGURE 5

Variation in resulting effect sizes among BTU definition methods. Box‐plots showing R ²adj values in the case of GM‐host species specificity and Procrustes correlation coefficients in the case of P_COMP or P_PHYLO using different methods of BTU definition. Outliers are marked by black points