Microbiome and environment explain the absence of correlations between consumers and their diet in Bornean microsnails.

Classical ecological theory posits that species partition resources such that each species occupies a unique resource niche. In general, the availability of more resources allows more species to co-occur. Thus, a strong relationship between communities of consumers and their resources is expected. However, correlations may be influenced by other layers in the food web, or by the environment. Here we show, by studying the relationship between communities of consumers (land snails) and individual diets (from seed plants), that there is in fact no direct, or at most a weak but negative, relationship. However, we found that the diversity of the individual microbiome positively correlates with both consumer community diversity and individual diet diversity in three target species. Moreover, these correlations were affected by various environmental variables, such as anthropogenic activity, habitat island size, and a possibly important nutrient source, guano runoff from nearby caves. Our results suggest that the microbiome and the environment explain the absence of correlations between diet and consumer community diversity. Hence, we advocate that microbiome inventories are routinely added to any community dietary analysis, which our study shows can be done with relatively little extra effort. Our approach presents the tools to quickly obtain an overview of the relationships between consumers and their resources. We anticipate our approach to be useful for ecologists and environmentalists studying different communities in a local food web.

Introduction

Different species within communities of ecologically similar species (guilds) can avoid competition through niche partitioning (Gause 1934, Hutchinson 1961). An important dimension of the species’ niche is formed by the resources a species can harvest (Whittaker 1972). The classical, fundamental model (i.e., each species being limited by a single resource) predicts that a greater diversity of available resources allows for more species to coexist in a local community (Tilman and Pacala 1993). Thus, a strong relationship between the diversity of a consumer community and its resources is expected (Hutchinson 1959, MacArthur 1965), and many experimental and observational studies have confirmed a positive correlation between the diversity of consumers and producers (e.g., between plant community richness and insect (Knops et al. 1999, Haddad et al. 2001), arthropod (Siemann et al. 1998, Haddad et al. 2009), butterfly (Hawkins and Porter 2003), and bird community diversities (Kissling et al. 2007). In turn, when coexisting species have a preference for the same resources, this may result in interspecific competition, that can influence the use of these resources through behavioral character displacement (Husar 1976, Moosman et al. 2012). This is commonly known from its most extreme manifestation, “ecological release,” which occurs after competitors have been removed (Kohn 1978, Kernaléguen et al. 2015). However, in many such cases it remains unclear to what extent such correlations between consumers and resources are causal (i.e., direct) or instead the result of both communities responding to yet another biotic (e.g., the microbiome) or abiotic variable the same way.

Spatial and temporal variations in abiotic factors are important drivers of diversity, with heterogeneous environments allowing more species to co‐occur than homogeneous environments (Tilman and Pacala 1993). For example, Hawkins and Porter (2003) found plant and butterfly diversity to be positively correlated, and identified variation in primary productivity and topographical variability as the most likely cause. Longmuir et al. (2007) found no correlations between communities of producers, consumers, and bacteria, but instead associations with environmental variables were important. Similarly, a third community, such as those of predators or pathogens, may influence interactions between two communities (Tilman and Pacala 1993). This is well known in tritrophic interactions of plants with herbivorous pests being mitigated by attracted predators (Heil 2008).

The microbiome has recently gained much attention as an important companion to its host, and it was even argued that the “holobiont” (i.e., the symbiotic assemblage of the host and its microbiome) is the actual unit of selection (Bordenstein and Theis 2015). Various studies have shown how the microbiome can adapt to the host’s diet and become an important codriver of diet choice (Muegge et al. 2011, Colman et al. 2012, Youngblut et al. 2019), thus directly influencing consumer–resource interactions.

Given the usually close relationship between consumers and their resources, we aimed to answer the question of how consumers’ individual diets are influenced by the community in which these individuals live, in the light of a third community, the individual’s microbiome. To this end, we studied communities of land snails (Gastropoda) living in an archipelago‐like environment of habitat islands of limestone bedrock in lowland tropical rainforest in Sabah, Malaysian Borneo (Schilthuizen et al. 2003, Schilthuizen 2011, Hendriks et al. 2019, ). We used a combination of snail consumer community census data and metabarcoding data from seed plant and bacterial genetic markers obtained from individual snail guts to reconstruct snail community, seed plant diet, and microbiome richness and diversity, and to study relationships between communities. In addition, we studied community‐level responses to the environment, because the influence of the environment on community interactions, as described above, is often important (Tilman and Pacala 1993, Hawkins and Porter 2003, Longmuir et al. 2007). Based on the MacArthur and Wilson (1963) equilibrium theory of island biogeography, we expected more dispersal and colonization (i.e., larger habitat island; less isolation; closer to a river) to result in higher community diversities. Because of the hydrophilic nature of land snails, we expected humidity during sampling to have a positive effect on community diversity (Martin and Sommer 2004). The presence of cave entrances, because of possible eutrophic conditions caused by runoff from bat and swiftlet guano, was expected to have a positive influence on plant diversity and therefore consumer communities (Sánchez‐Piñero and Polis 2000). Our expectations of the influence of anthropogenic presence were ambiguous (Luck 2007), but deemed important to include for conservation reasons.

We find no (or a weak but negative) direct correlation between consumer community and individual seed plant diet diversities. However, diversity of the third community, the individual microbiome, shows positive correlations with diversity of both the consumer community and the individual seed plant diet in three target species of consumers. Furthermore, habitat island size and distance to cave entrances have a positive influence on the diversity of all communities, whereas distance to anthropogenic activity and current humidity have a negative influence. We conclude that the microbiome and the environment are important predictors for the presence or absence of correlations between consumer and diet communities in a food web, and should be included in studies of consumers and their resources.

Methods

Study system and sampling

We studied species‐rich consumer communities of land snails, their seed plant diets, and microbiomes, in the Lower Kinabatangan Floodplain in Sabah, Malaysian Borneo (Schilthuizen et al. 2003, Schilthuizen 2011, Hendriks et al. 2019, ). Most snail community members have a preference for calcium carbonate as a substratum or depend on it. As such, they are mainly restricted to the scattered limestone outcrops within the tropical rainforest of our study region (Schilthuizen et al. 2003). In this study, we focused on six different limestone outcrops (of ~20 in the region; Fig. 1A) and collected samples from three plots per outcrop along its base (four for location Keruak). Each plot measured 2 × 2 m, and was at least 50 m from the next plot.

Fig. 1

(A) Sampling locations (i.e., limestone outcrops; in black and named) in the Lower Kinabatangan Floodplain (in blue), Sabah, Malaysian Borneo; unsampled locations in gray. Inset map © freevectormaps.com. (B) Five consumer community species. Left to right: Alycaeus jagori Von Martens, 1859, Georissa similis E. A. Smith, 1893, Plectostoma concinnum (Fulton, 1901), Diplommatina calvula Vermeulen, 1993, and Kaliella accepta (Smith, 1895). Drawings: Bas Blankevoort. Scale bars equal 1 mm.

Snail communities were censored by collecting and sorting empty snail shells from 5 L of soil debris within each plot, with the collection of shells serving as proxy for the community (below referred to as “shell consumer community”; Liew et al. 2008), in 2015 and 2016. (For the rationale of using these shells to represent contemporary communities, see Appendix S1: Methods S1, Table S1.) Shells were identified to species level based on the latest taxonomic literature (Vermeulen et al. 2015, Liew 2019, ) and counts per plot and species were collected in a community matrix (Appendix S1: Table S2A). Live snails were collected for the study of their microbiome and their seed plant diet when plots were revisited (relocated based on GPS readings and photos) in 2017 (Appendix S1: Table S3). All collections were made within a time frame of 10 d to exclude the influence of seasonal variation. We focused on three unrelated but mostly omnipresent target species: Alycaeus jagori Von Martens, 1859, Georissa similis E. A. Smith, 1893 s.l., and Plectostoma concinnum (Fulton, 1901), and aimed to collect 40 individuals for each species per plot. Additionally, we collected all live snails belonging to other species that we encountered, with a maximum of 20 individuals per species per plot (Fig. 1B). From these live snails we created a second community matrix (“live consumer community”; Appendix S1: Table S2B) with which we redid all analyses performed with the “shell consumer community” data, because these “live consumer community” data directly represented the individuals for which we also collected metabarcoding data. Also, these live individuals definitely co‐occurred and were therefore active at the same moment (daytime and season), which need not be the case for the community matrix from snail shells. Details on the target species and collection procedures can be found in Hendriks et al. (2019). We constructed a phylogeny for the snail species with five commonly used barcode markers (16S, 18S, 28S, COI, H3; Webster et al. 2012; Appendix S1: Methods S2, Fig. S1, Table S4).

To test for the influence of environmental variables associated with dispersal and colonization, we collected data from GoogleEarth on habitat island size (i.e., limestone outcrop surface as seen from above; outcrops were clearly visible on aerial photographs), isolation (shortest distance from plot to next limestone outcrop), and shortest distance to a probable vector of dispersal, the Kinabatangan River (Appendix S1: Table S3). Similarly, we collected data on variables probably associated with habitat suitability, namely anthropogenic distance (distance to closest road and plantation; also from GoogleEarth), current humidity (time since last rainfall event and the level of humidity; both scored during fieldwork), and shortest distance to the nearest cave entrance (based on data from Schilthuizen and Njunjić 2019), with bat‐ and swiftlet‐inhabited caves considered a possible heavy nutrient source (Sánchez‐Piñero and Polis 2000, Gagnon et al. 2013, Vizzini et al. 2016). No substantial changes to the outcrops (clearing of more than 1% of the surface) were recorded during fieldwork.

Metabarcoding and bioinformatics

We performed metabarcoding of seed plant rbcL and bacterial 16S rRNA genes on a maximum of 20 individuals per snail species per plot to represent the individual seed plant diet and microbiome, respectively. DNA extraction and library preparation details are described in Appendix S1: Methods S3. Raw sequence data from metabarcoding were pooled by marker, resulting in single pools for rbcL and 16S rRNA. We used QIIME2 v2017.12 (Bolyen et al. 2018) with the DADA2 (Callahan et al. 2016) philosophy and routine to denoise, apply quality control, and export representative Amplicon Sequencing Variants (ASVs; for QIIME2 scripts, see Appendix S1: Methods S4). Plant (rbcL) and bacterial (16S rRNA) origin of ASVs were confirmed by blasting our newly collected results against classifiers built from data from Bell et al. (2017) and GreenGenes v13.8 (DeSantis et al. 2006), respectively. Any ASVs found in one or more negative controls were considered to possibly originate from contamination and therefore subsequently removed from all samples. These were 40 out of 778, and 192 out of 19,542 ASVs from the rbcL and 16S rRNA data sets, respectively. Additionally, 16S rRNA data were checked not to be of host (i.e., gastropod) origin using NCBI’s nucleotide BLAST search in Geneious v9.1.6.10 For both markers, QIIME2 ASV alignments were checked by eye, and nonaligning reads removed and double‐checked to be of nontarget origin using an individual nucleotide BLAST search, after which alignments were updated using MAFFT v1.3.5 (Katoh and Standley 2013) in Geneious. Maximum‐likelihood phylogenetic trees from ASV alignments for rbcL and 16S rRNA were constructed using FastTree v1.0 (Price et al. 2009) with default settings in Geneious.

Statistical analyses

Community and phylogenetic data were imported into R v3.5.0 (R Core Team 2018) and combined into separate objects using the package “PhyloSeq” v1.24.2 (Mcmurdie and Holmes 2013): shell consumer community, live consumer community, seed plant diet, and microbiome. In the last two, the “community” is represented by the collection of ASVs from the individual diet/microbiome. For each community, we calculated the following metrics: Chao1 richness (i.e., data rarefied to account for unequal sample sizes), Shannon and Simpson diversities using “PhyloSeq”, Shannon evenness as defined by Magurran and McGill (2011), and Faith’s phylogenetic diversity (Faith 1992) using the package “picante” (Kembel et al. 2010).

To study diet and microbiome differentiation among species and locations, we performed ordination by nonmetric multidimensional scaling (nMDS) on four commonly used distance metrics (Bray‐Curtis, Jaccard, weighted UniFrac, and unweighted UniFrac). To do so, we pooled data per target species (plus a pool of all nontarget species lumped together) and by sampling plot (excluding species–plot combinations with ≤1 individual) and plotted results with 95% confidence levels by species using “ggplot2” v3.2.1 (Wickham 2016). We tested for the influence of “snail species” and “location” using PERMANOVA (Anderson 2017) through the function “adonis” from the package “vegan” (Oksanen et al. 2017). To test for correlations between diet and microbiome distances, we performed Mantel tests on all data (i.e., from all species together) using the function “mantel” from the package “vegan” with 999 permutations. We tested for correlations at both the individual level (i.e., for all individual snails for which we obtained data for diet and microbiome; n = 643) and at the plot level (i.e., pooling data by plot; n = 14). Acknowledging likely species differences, we repeated these tests specifically for each of the target species.

We obtained an initial, general idea of the correlations between Shannon diversity of consumer community (from shells), consumer diets, and consumer microbiomes using simple linear regression (using “ggpubr” v0.2.3; Kassambara 2017), for all individuals in our study, first ignoring unbalanced sampling and the effect of the variable “snail species.” Correlations were tested by plot (i.e., taking the “community diet” and “community microbiome” by pooling data by plot, using bootstrapping of the data to account for different sample sizes by plot, with sample size equal to lowest plot sample total, n = 14) and individual snail. Next, we studied the influence of the shell consumer community and “snail species” (fixed effects) on the seed plant diet and microbiome (response variables) using generalized linear mixed models (GLMM; Bolker et al. 2009), focusing only on the target species (for which sample sizes were large, i.e., >150). Because we sampled multiple individuals of the same species from each plot, we added “plot” as random effect to account for pseudoreplication. This resulted in a total of six models: three models each for diet and microbiome, using either Shannon diversity, PD, or Chao1 richness for both the response and explanatory variables. For each of these models, we determined the best‐fit distribution for the response variable using the package “fitdistrplus” v1.0.14 (Delignette‐Muller 2015) and used the function “glmmTMB” from the package “glmmTMB” v0.2.3 (Brooks et al. 2017) with default settings to fit the model. We compared each full model to simpler, nested models, and selected the best model based on the lowest AIC value from ANOVA. We performed a post hoc test on each full model to study any differences (“contrasts”) in the response to variations in the shell consumer community between the target species (adjusted for multiple comparisons based on Tukey’s method), using the function “emtrends” from the package “emmeans” v1.3.2 (Lenth 2019). We repeated the above routine for each target species separately, now with “location” as a fixed effect. In these “species models,” locations with fewer than 10 individuals for the species studied were not considered, because of convergence issues in model fitting because of too few data points. Contrasts were now studied among locations.

To study the community data in concert with the environmental variables, we applied partial‐least‐squares path modeling (Sanchez 2013) using the R package “plspm” v0.4.9 (Sanchez et al. 2017). We first studied “core models” with three latent variables (LVs) corresponding to the three communities: consumer, seed plant diet, and microbiome. LVs were constructed from combinations of summary metrics Shannon and Simpson diversity, Faith’s phylogenetic diversity (PD), and Shannon evenness (Appendix S1: Table S5). Model assessment followed Sanchez (2013). In short, we checked the “outer model” (i.e., the loading of LVs with data from measurements) for unidimensionality (which confirms that each of the LVs are loaded only with variables that describe the same phenomenon), cross loadings (to confirm that variables explain their own LV best, and not another, which would indicate an erroneous assignment of variable to LV), and positive loadings of the outer model (to confirm variables point in the same direction); negative loadings were encountered for Shannon evenness and resolved by taking negative values for this metric. The “inner model” (i.e., the relationships between the LVs) was checked for communality (to measure the part of LV’s variance explained by its variables), redundancy (how much of each endogenous LV is explained by the other LVs), and overall goodness‐of‐fit (to assess overall model quality). We performed bootstrapping to assess significance by rerunning models 999 times with data from 300 random individuals (without replacement) from any species for which we had both diet and microbiome data (so‐called “complete models”). Because unequal sample sizes per species in this model could bias the outcome of model predictions, we repeated the above routine for the three target species only (for which sample sizes were large enough) with equal sample sizes of 100 individuals/target species (so‐called “normalized models”). We calculated mean goodness‐of‐fit and path loadings and their significance (P < 0.05, two‐tailed test). In addition to these core models, we ran “full models” (again, complete and normalized) with the environmental variables mentioned above included as LVs. LVs were constructed from single (habitat island size, next habitat island distance, and cave distance), or (in three cases) multiple measurements: river distance (constructed from shortest distance to the river and altitude of the plot), anthropogenic distance (shortest distances to main road and plantation), and current humidity (time since rain and humidity level as scored during sampling on a scale from 1 to 4 with increasing humidity; not modeled against the shell consumer community, because these two cannot be correlated biologically). Because there was no a priori knowledge on what consumer community data (from shells or live data) represents the “true” community best, we ran each of the above models for both data sets.

Results

We collected and identified 11,833 empty snail shells from 19 plots, belonging to 55 species, with a mean of 657 individuals per plot (range 43–4,353; Appendix S1: Table S2A); in addition, we collected 1,494 live snails from 28 species, with a mean of 62 individuals per plot (range 33–139; Appendix S1: Table S2B). From these live snails, we randomly selected individuals for metabarcoding (aimed at 20 individuals/species/plot; total ~840 individuals), and obtained successful metabarcoding results for seed plant diet and/or microbiome for 820 individuals (Table 1). Metabarcoding was less successful for the diet (data for 645 individuals) than for the microbiome (data for 818 individuals; Appendix S1: Table S5), which was probably a result of less plant than microbial DNA from extractions. The various metrics by individual showed that the microbiome was more diverse and richer than the diet, and individual variation (also intraspecific) was substantial (Shannon diversity: diet mean 0.59 [SD 0.57] vs. microbiome 2.91 [0.85]; Simpson diversity: diet 0.32 [0.29] vs. microbiome 0.86 [0.10]; Faith’s PD: diet 0.45 [0.39] vs. microbiome 15.18 [8.02]; Shannon evenness: diet 4.96 [24.05] vs. microbiome 0.10 [0.05]; Chao1 richness: diet 3.32 [3.10] vs. microbiome 74.11 [56.53]; Appendix S1: Table S5).

Table 1

Summary of sample sizes for which successful metabarcoding data were obtained for microbiome (n = 818) and diet (n = 645), or either (as listed in the table; n = 820), sorted by species, location, and plot.

Family	Species	Batangan			Keruak				Pangi			Tandu Batu			Tomang‐gong 2			Tomang‐gong Kecil			Totals
		Plot 2	Plot 5	Plot 6	Plot 1	Plot 3	Plot 6	Plot 7	Plot 2	Plot 5	Plot 6	Plot 3	Plot 6	Plot 7	Plot 1	Plot 4	Plot 5	Plot 1	Plot 3	Plot 6
Ariophantidae	Everettia sp.											2						2			4
	Microcystina appendiculata (Von Moellendorff, 1893)										1							1			2
	Macrochlamys tersa (Issel, 1874)													1	1				2		4
Assimineidae	Acmella cyrtoglyphe Vermeulen, Liew & Schilthuizen, 2015								3												3
	Acmella striata Vermeulen, Liew & Schilthuizen, 2015											6	5	1							12
Cyclophoridae	Alycaeus jagori Von Martens, 1859				1				20	21	15		20	20		2	20	20	18	6	163
	Chamalycaeus sp.										1								1		2
	Japonia kinabaluensis (E. A. Smith, 1895)	1	2																		3
	Japonia sp.	2	1					2													5
	Leptopoma pellucidum (Grateloup, 1840)													1							1
	Leptopoma sericatum (Pfeiffer, 1851)	9		3	1						1		1								15
	Pterocyclos/Opisthoporus sp.		1										1	1				4	1		8
Diplommatinidae	Diplommatina asynaimos Vermeulen, 1993		1																		1
	Diplommatina calvula Vermeulen, 1993						2	1													3
	Diplommatina gomantongensis (E. A. Smith, 1894)					1															1
	Diplommatina rubicunda (Von Martens, 1864)							1									1	1		3	6
	Plectostoma concinnum (Fulton, 1901)	20	20	20	20	20	20	20	1	20	20	20	11	22	18	20	20	20	20	20	352
	Plectostoma simplex (Fulton, 1901)											3	11								14
Euconulidae	Kaliella accepta (Smith, 1895)					2		1		1					1			1			6
	Kaliella barrakporensis (Pfeiffer, 1852)														2						2
	Kaliella calculosa (Gould, 1852)														1						1
	Kaliella scandens (Cox, 1872)				1										2						3
Helicinidae	Sulfurina martensi (Issel, 1874)	6	2	4			2													1	15
Hydrocenidae	Georissa kinabatanganensis Khalik, Hendriks, Vermeulen & Schilthuizen, 2018				1	10															11
	Georissa similis E. A. Smith, 1894 s.l. ‡		1	9	12	20		4	20	16			16		20	13	5	20	19	2	177
Rathouisiidae	Atopos sp.										1										1
Trochomorphidae	Videna metcalfei (Pfeiffer, 1845)	1	1					1						1							4
	Videna sp.																	1			1
Totals		39	29	36	36	53	24	30	44	58	39	31	65	47	45	35	46	70	61	32	820

The three target species are printed in bold. See Appendix S1: Tables S2, S5 for sample details.

Originally described (and collected by us) as Georissa similis E. A. Smith, 1894, but recently split into a radiation of highly similar and closely related taxa (Khalik et al. 2019). With all phylogenetic relations within the radiation being much closer than those among all other taxa considered within this study, with the exception of G. nephrostoma Vermeulen, Liew and Schilthuizen, 2015, we treat G. similis s.l. as a single species in this study.

Consumer–diet correlations

Linear regression showed no significant correlation between consumers and diet diversity, at both the plot and individual levels (Fig. 2; see Appendix S1: Fig. S2 for general trends in PD and Chao1). Model selection by GLMMs confirmed the absence of a significant effect of the shell consumer community on the individual seed plant diet Shannon diversity for the three target species (Appendix S1: Table S6). Individuals of Alycaeus had the most diverse diets, with individuals from the other two target species having significantly lower diet diversities (Table 2; Appendix S1: Fig. S3A). Model selection and best models of GLMMs for Faith’s PD and Chao1 richness showed a significantly positive influence of the shell consumer community on the individual diet diversity of A. jagori only (Appendix S1: Table S7 and Fig. S3B, C). Differences among several locations were significant for A. jagori (PD and Chao1 richness) and P. concinnum (PD), with no clear geographical pattern (Appendix S1: Table S8–S10, Fig. S4A–C).

Fig. 2

Comparisons between consumer community, diet, and microbiome Shannon diversities for all samples from all species studied. (A) Diet, microbiome, and consumer community by plot, showing mean values from 1,000 bootstrapped data sets (at equal sample size by plot of n = 14), with 95% confidence intervals. (B) Diet and microbiome by individual, against consumer community by plot. Dashed lines show results from simple linear regressions. Note that an individual diet Shannon diversity of zero indicates a single diet item found from the respective individual snail.

Table 2

GLMM best model results, using function “glmmTMB” from R package “glmmTMB” v0.2.3 (Brooks et al. 2017).

Response	Metric	n	Coefficient	Estimate	SE	z‐Value	P
Diet	Diversity (Shannon)†	357	Intercept	0.98	0.04	25.93	<0.001
			Species (Georissa similis s.l.)	−0.27	0.07	−3.91	<0.001
			Species (Plectostoma concinnum)	−0.20	0.05	−4.08	<0.001
	Phylogenetic diversity (Faith’s PD)‡	357	Intercept	−1.77	0.41	−4.31	<0.001
			Shell consumer community phylogenetic diversity (PD)	0.56	0.18	3.06	0.002
			Species (G. similis s.l.)	0.74	0.81	0.91	0.361
			Species (P. concinnum)	0.84	0.45	1.86	0.063
			Shell consumer community phylogenetic diversity (PD) × Species (G. similis s.l.)	−0.85	0.41	−2.06	0.039
			Shell consumer community phylogenetic diversity (PD) × Species (P. concinnum)	−0.64	0.21	−3.06	0.002
	Richness (Chao1)‡	539	Intercept	0.75	0.29	2.57	0.010
			Shell consumer community richness (Chao1)	0.03	0.01	2.78	0.005
			Species (G. similis s.l.)	−0.31	0.47	−0.66	0.512
			Species (P. concinnum)	0.57	0.31	1.86	0.062
			Shell consumer community richness (Chao1) × Species (G. similis s.l.)	−0.02	0.02	−1.20	0.230
			Shell consumer community richness (Chao1) × Species (P. concinnum)	−0.05	0.01	−4.13	<0.001
Microbiome	Diversity (Shannon)‡	690	Intercept	0.98	0.14	7.01	<0.001
			Shell consumer community diversity (Shannon)	0.13	0.08	1.58	0.114
			Species (G. similis s.l.)	0.01	0.11	0.11	0.915
			Species (P. concinnum)	−0.21	0.10	−2.10	0.036
			Shell consumer community diversity (Shannon) × Species (G. similis s.l.)	−0.12	0.06	−2.13	0.033
			Shell consumer community diversity (Shannon) × Species (P. concinnum)	0.01	0.06	0.13	0.899
	Phylogenetic diversity (Faith's PD)§	690	Intercept	26.66	4.64	5.75	<0.001
			Shell consumer community phylogenetic diversity (PD)	−5.59	2.20	−2.54	0.011
			Species (G. similis s.l.)	6.67	3.61	1.85	0.064
			Species (P. concinnum)	−10.56	2.90	−3.64	<0.001
			Shell consumer community phylogenetic diversity (PD) × Species (G. similis s.l.)	−1.74	1.67	−1.04	0.297
			Shell consumer community phylogenetic diversity (PD) × Species (P. concinnum)	4.15	1.36	3.06	0.002
	Richness (Chao1)§	690	Intercept	126.52	29.74	4.25	<0.001
			Shell consumer community richness (Chao1)	−1.86	1.18	−1.57	0.116
			Species (G. similis s.l.)	45.04	23.69	1.90	0.057
			Species (P. concinnum)	−50.29	18.13	−2.77	0.006
			Shell consumer community richness (Chao1) × Species (G. similis s.l.)	−1.84	0.93	−1.98	0.048
			Shell consumer community richness (Chao1) × Species (P. concinnum)	1.24	0.68	1.82	0.069

Significant terms in bold. For model selection, see Appendix S1: Table S6. Note that Shannon diversity and Faith's phylogenetic diversity could not be calculated when Chao1 equals one, hence smaller sample sizes for these metrics in diet data.

Response variable modeled as a normal distribution.

Response variable modeled as a lognormal distribution.

Response variable modeled as a gamma distribution.

Core path models on shell consumer community showed no significant influence of consumer community on individual diet diversities (Fig. 3A; Appendix S1: Table S11A) for both the complete model (including all species) and the normalized models (target species only; Appendix S1: Table S11A). However, the full normalized path model based on live consumer community suggested a significant negative influence (Fig. 3B; Appendix S1: Table S11B).

Fig. 3

Significant results from partial least‐squares path modeling (PLS‐PM; Sanchez 2013) for normalized models (i.e., with equal sample sizes from target species Alycaeus jagori Von Martens, 1859, Georissa similis E. A. Smith 1893 s.l., and Plectostoma concinnum (Fulton, 1901)), based on (A) shell consumer community and (B) live consumer community. Black triangles represent the core model; black arrows represent the full model; dashed lines highlight negative path coefficients. Labels represent significant mean path coefficients (P < 0.05 from 999 bootstraps). Nonsignificant results not shown, but available in Appendix S1: Table S11.

Ordination of the diet showed much overlap between target species and different locations for all four distance metrics studied (Appendix S1: Fig. S5A–D, Table S12), indicating little diet differentiation overall. However, some differentiation was found based on unweighted UniFrac distance (Fig. 4A), with both “species” and “location” significantly explaining part of these differences (Table 3).

Fig. 4

nMDS ordination plots based on unweighted UniFrac distance from sample data pooled by species and plot for the three target species, Alycaeus jagori Von Martens, 1859, Georissa similis E. A. Smith, 1893 s.l., and Plectostoma concinnum (Fulton, 1901), and the nontarget species lumped together. (A) Diet (n = 58 from 643 snails) and (B) microbiome (n = 59 from 815 snails). Plots for which only data on one snail/species plot was available (singletons) were excluded. Numbers refer to plot identity with each location, for details of which, see Appendix S1: Table S3. Ellipses indicate 95% confidence levels. See Appendix S1: Fig. S5 for results from other metrics.

Table 3

Results from PERMANOVA and BETADISPER analyses of unweighted UniFrac data from sample data pooled by species and plot for the three target species, Alycaeus jagori Von Martens, 1859, Georissa similis E. A. Smith, 1893 s.l., and Plectostoma concinnum (Fulton, 1901), and the nontarget species lumped together.

Response variable	Explanatory variable	df	PERMANOVA				BETADISPER
			SS	pseudo‐F	R 2	Pr (>F)	SS	pseudo‐F	Pr (>F)
Diet	Species	3	1.985	2.894	0.136	<0.001	0.020	2.027	0.116
	Location	5	1.497	1.310	0.103	0.040	0.008	0.562	0.726
	Species × Location	13	2.834	0.954	0.195	0.665
	Residuals	36	8.229		0.566
	Totals	57	14.544		1.000
Microbiome	Species	3	2.466	3.700	0.156	<0.001	0.050	10.766	<0.001
	Location	5	1.782	1.604	0.113	<0.001	0.005	0.390	0.865
	Species × Location	13	3.362	1.164	0.212	0.005
	Residuals	37	8.222		0.519
	Totals	58	15.832		1.000

Statistical testing based on 4,999 permutations. Note that a significant result from PERMANOVA may be indicative of differences in dispersion (and not just centroid differences) in the distance space, for those cases where a significant result from BETADISPER for that factor also occurs. Significant PERMANOVA terms that did not also show significant differences in dispersion are directly interpretable as a shift in community structure (i.e., a centroid shift only), and shown in bold. See Appendix S1: Table S12 for results from other metrics.

Consumer–microbiome correlations

Linear regression showed a positive correlation between Shannon diversity of the microbiome and the consumer community by plot (though nonsignificant), as well as by individual (P < 0.001; Fig. 2). Model selection by GLMMs showed that the best models included interaction effects between the shell consumer community (by plot) and the species on the individual microbiome Shannon diversity, Faith’s PD, and Chao1 richness for the three target species (Table 2; Appendix S1: Table S6). Individual microbiome Shannon diversity increased with shell consumer community diversity, although not significantly (Table 2; Appendix S1: Fig. S3D). In contrast, correlations in Faith’s PD and Chao1 showed a negative trend, the former significant for A. jagori and G. similis s.l. and the latter for G. similis (Table 2; Appendix S1: Table S7, Fig. S3E–F). Significant differences among several locations were found in Shannon diversity for A. jagori (Appendix S1: Tables S8–S10, Fig. S4D–F). There was no clear pattern among the locations in differences between target species’ responses (Appendix S1: Table S10).

In contrast to the GLMMs, core path models on shell consumer community diversity showed a significantly positive relationship between consumer community and microbiome diversity, for both complete (all species; Fig. 3A) and normalized models (target species only; Appendix S1: Table S11A). In full‐path models and models based on live consumer community (Appendix S1: Table S11B), this correlation was also positive but not significant.

Ordination of the microbiome showed each target species to occupy a subset of the overall microbiome niche space for all four distance metrics studied, with some overlap between the A. jagori and P. concinnum (Fig. 4; Appendix S1: Fig. S5E, H). For all four distance metrics studied, “location” significantly explained part of the variation (Table 3; Appendix S1: Table S12).

Microbiome–diet correlations

Linear regression showed virtually no effect of Shannon diversity of the microbiome on the diet by plot, but a significant positive effect at the individual level (P < 0.001; Fig. 2). All path models confirmed a significant positive relationship between individual microbiome and diet diversities (Fig. 3; Appendix S1: Table S11). The mean path coefficient was highest in the complete core model (0.198) and lowest in the normalized full model (0.105), in both cases based on the live consumer community. Results from Mantel testing of unweighted UniFrac distances among both samples (individual snails) and plots further confirmed the positive correlation between the microbiome and the diet, whereas weighted UniFrac distances did not (Table 4). Correlations are strongest at the plot level (i.e., diet and microbiome data pooled by plot). Furthermore, we found important differences between the target species, with no correlation found for G. similis s.l.

Table 4

Results from Mantel tests of diet versus microbiome distances based on unweighted UniFrac and weighted UniFrac, respectively.

	Unweighted UniFrac		Weighted UniFrac
	r	P	r	P
By sample
All species together	0.078	0.001	−0.014	0.773
Alcaeus jagori	0.102	0.034	0.063	0.097
Georissa similis s.l.	0.000	0.477	0.048	0.154
Plectostoma concinnum	0.070	0.032	−0.010	0.621
By plot
All species together	0.547	0.001	0.112	0.200
A. jagori	0.516	0.036	0.083	0.318
G. similis s.l.	−0.119	0.710	0.131	0.257
P. concinnum	0.305	0.002	−0.008	0.501

Tests performed both by sample (i.e., by individual snail, for all samples for which both data sets were available; n = 644), and by plot (with data for all individual snail pooled; n = 19). All tests were repeated by species for the three target species, Alycaeus jagori Von Martens, 1859, Georissa similis E. A. Smith, 1893 s.l., and Plectostoma concinnum (Fulton, 1901). Significant results (P < 0.05) in bold.

Environmental effects

Full path models (i.e., models including environmental variables) showed higher goodness‐of‐fit values than core models (mean 0.389 vs. 0.152, respectively, for normalized models based on shell consumer community; Appendix S1: Table S11A), and thus explained more of the variation in the LVs (note that these goodness‐of‐fit values are only indicative and not appropriate for detailed model comparison; cf. Sanchez 2013). In the full models, correlations between core LVs (i.e., consumer, seed plant diet, and microbiome) were lower than in core models, and part of the explanatory power was “moved” to the environmental LVs (Fig. 3).

Several environmental variables were significantly correlated with community LVs (Fig. 3). We found in general that (1) diversities of consumer, seed plant diet, and microbiome communities were higher when closer to human activity (i.e., diversity decreased with anthropogenic distance), (2) larger habitat islands supported more diverse snail consumer communities with a more diverse diet and microbiome, and (3) snail consumer and diet diversity were lower near caves. Other variables (i.e., current humidity, distance to the nearest habitat island, and distance to the river) influenced one or more of the investigated communities. More humid conditions negatively influenced individual microbiome and diet diversity, more isolated locations had higher snail consumer and lower individual microbiome diversity, and more diverse diets were found further away from the river. These results are based on normalized models (i.e., used equal sample sizes of target species, for which we had most data), unless stated otherwise. Results from complete models (based on all species) were similar, but path coefficients were less often significant (Appendix S1: Table S11).

Discussion

We studied the correlations between the diversity and richness of three closely linked communities on habitat islands of limestone bedrock in lowland rainforest in Sabah, Malaysian Borneo: the communities of land snail consumers, their members’ seed plant diets, and their microbiomes. We aimed to answer the question of how consumer diets are influenced by the communities they are part of, and in addition we tested for the influence of each consumer’s microbiome and environmental variables. We found no, or at most a weak but negative, direct relationship between consumer community and individual diet (for the shell and live consumer community, respectively). However, we found both the diversity of shell consumer community and individual diet to be positively correlated with individual microbiome diversity. In some cases, responses were different for the three consumer target species studied. Moreover, environmental variables affected core community correlations, most notably anthropogenic activity, distance to the nearest cave, and habitat island size.

Our finding of no (or a weak negative) relationship between consumer community and individual seed plant diet diversity seems to contradict classical theory (Hutchinson 1959, MacArthur 1965) and general empirical findings (Siemann et al. 1998, Knops et al. 1999, Haddad et al. 2001, 2009, Kissling et al. 2007), although in our study we cannot know what portion of the resource community was not eaten. It is possible that other such cases of no or weak direct correlation have either simply rarely been published due to publication bias towards positive results (Knight 2003, Fanelli 2012), or are truly rare. However, in a study very similar to ours, dealing with three trophic communities (pelagic zooplankton, phytoplankton, and bacteria), no significant correlations between community diversities were detected either (Longmuir et al. 2007). Several large‐scale studies on aquatic ecosystems have also reported the absence of community correlations (Allen et al. 1999, Irigoien et al. 2004, Declerck et al. 2005). Declerck et al. (2005) noted that correlations are generally weaker in aquatic systems because their consumers are often filter‐feeders and as such less specific in their prey choice.

Our findings support a suggestion from Schilthuizen (2011), namely, that tropical snail communities are not strongly influenced by external biotic diversity, but instead more by “available microclimatic and microchemical gradients.” This was also found for butterflies by Hawkins and Porter (2003): the butterfly and plant diversities they studied were both influenced by the same environmental factors, namely primary productivity and topographical variability. The weak negative correlation we found shows that the more species are present in the community, the smaller the individual diet. This could result from competition for food sources, potentially leading to smaller realized niche widths. This is in fact the converse of “ecological release”: when competing species are removed, the remaining species become able to consume a wider array of resources (Kohn 1978, Kernaléguen et al. 2015).

Although we found no (or a weak) direct relationship between consumer community diversity and individual diet diversity, we did find a significant positive correlation between individual microbiome and diet diversity. The microbiome is partly composed of bacteria inherited from the parents to the offspring (vertical transfer), but also of bacteria obtained from food sources and community members (horizontal transfer), with a likely selection of (nutritionally) beneficial bacteria by the host (Watkins and Simkiss 1990, Engel and Moran 2013, Seedorf et al. 2014, Macke et al. 2017). Furthermore, fecal transplant studies (Kohl et al. 2014, Kohl and Dearing 2016) and field studies (Kohl et al. 2018) have shown a positive effect of the microbiome on host fitness, partly by allowing individuals to digest previously unpalatable food. Because different host species usually contain species‐specific microbiomes (Glasl et al. 2018, Hird et al. 2018, Sörenson et al. 2019), a more diverse consumer community could result in a more diverse individual microbiome via horizontal transmission from coexisting species. This could indirectly allow individuals to feed on a more diverse set of plants, or potentially get other benefits, such as better immunity. The latter hypothesis suggests, however, a correlation between consumer diversity and diet diversity that we did not find and is therefore less plausible. The positive correlation between microbiome and consumer diversity assumes that microbiome diversity can be used to assess generality of the diet, but no evidence was found in the literature. We suggest further experimental research in a more controlled environment (without potentially confounding effects) to explain our finding.

The influence of several environmental variables (mainly anthropogenic distance, habitat island size, and presence of cave entrances) on the core communities is often stronger than between core communities themselves. It is possible that correlations we found between communities are actually the result of similar responses to the environment (Allen et al. 1999), suggesting that resources are not a limiting factor. First, community diversities are higher in plots close to anthropogenic activity, which is likely a by‐product of horticultural and agricultural activities, which can increase tropical plant diversity (Fine 2002, Stadler et al. 2000). Similar data for snail communities are rare, but snail dispersal via the same vector is common (Dörge et al. 1999). Although richness of snail communities (in nonlimestone habitat) has been shown to drop towards agricultural activity in Nigeria (Oke and Chokor 2009), pulmonate snails (to which virtually all invasive species belong) in our study region have previously been shown to perform surprisingly well after (human) disturbance (Schilthuizen et al. 2005). Besides, anthropogenic activity may increase dispersal between limestone outcrops and in that way increase local species diversity (Cadotte 2006). Second, in agreement with the MacArthur and Wilson (1963) equilibrium theory of island biogeography, community diversities we studied increase with habitat island size. However, consumer community and individual diet diversities do not increase with decreasing distance to the next habitat island (should this be taken as a proxy for “the mainland” from the equilibrium theory). This is in line with recent work by Hendriks et al. (2019), who showed little effect of distance on the colonization opportunities in these snail consumers communities. Third, community diversities are lower close to cave entrances, which might be a result of the local eutrophic conditions caused by runoff from bat and swiftlet guano. This is in contrast with previous findings showing a positive influence of bird guano runoff on plant and consumer communities on oceanic islands (Sánchez‐Piñero and Polis 2000).

In our study we chose to investigate the influence of the consumer community on the microbiome and the diet, because the snail community is the only independent community as microbiome and diet were both reconstructed from the collected snails. We therefore oriented the arrows in the PLS path modeling in that way. However, we are aware that communities are in fact interacting with each other and that another direction in the path modeling could have been assumed. A likely alternative is that a richer plant diet allows for a richer microbiome (Reese and Dunn 2018), which in turn allows for a higher consumer diversity. For instance, it has been found that diet can rapidly alter the human (David et al. 2014) and murine gut microbiome (Carmody et al. 2019). A study in African megafauna revealed a strong correlation between the diet and microbiome compositions within and among species, but could not confirm a positive relation between dietary and microbial diversity (Kartzinel et al. 2019).

We presented the microbiome and seed plant diet diversity captured from both the entire consumer community (i.e., for all species present; complete model) and normalized data (including only the three target species; normalized model). We discussed the normalized data, as we expect unequal sample sizes (by host species) to bias the outcome of our study. We are aware that including only three species may have potential consequences. For instance, we found several different results for the different target species in the outcomes of the GLMMs, showing that different species respond differently to their surroundings. However, path model correlations between consumer, individual diet, and individual microbiome did not differ much between the complete and normalized models (with the influence of environmental variables being less often significant in the complete models), indicating that our findings are rather robust.

In conclusion, the correlation we found between the microbiome and the host snail community suggests that general ecological theories may not only hold at the host (community) level, but also at the microbiome level. Our results demonstrate that traditionally ignored food web layers (such as the microbiome) and environmental variables may explain correlations between different food web communities, which signifies the importance of the holobiont concept. With the ever‐decreasing costs of metabarcoding, we suggest the addition of microbiome inventories in any genetic dietary analysis. Our approach presents the tools needed to obtain an extended analysis of consumers and their resources in a local food web, which could greatly benefit ecologists and environmentalists.

Appendix S1

Click here for additional data file.

Data S1

Click here for additional data file.

Metadata S1

Click here for additional data file.