University of Groningen On the origin of species assemblages of Bornean microsnails Hendriks, Kasper

On the origin of species assemblages of Bornean microsnails

Hendriks, Kasper

DOI:

10.33612/diss.124819761

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from

it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date:

2020

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Hendriks, K. (2020). On the origin of species assemblages of Bornean microsnails. University of

Groningen. https://doi.org/10.33612/diss.124819761

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the

author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately

and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the

number of authors shown on this cover page is limited to 10 maximum.

between consumers and their diet

in Bornean snail-plant-interactions

Kasper P. Hendriks‡, Karen Bisschop‡, Hylke H. Kortenbosch,

James C. Kavanagh, Anaïs E. A. Larue, Chee-Chean Phung,

Dries Bonte, Elza J. Duijm, Joana Falcão Salles, Alex L. Pigot,

Francisco J. Richter Mendoza, Menno Schilthuizen,

Arjen G. C. L. Speksnijder, and Rampal S. Etienne

‡ Joint First Authors

Unpublished working manuscript,

under review with Ecology.

Abstract

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. Specifically, ecological communities of

consumers are expected to be strongly correlated to ecological communities of their

resources. However, these correlations may be influenced by other layers in the food

web, or by the environment. Here we show, by studying communities of consumers

(land snails) and communities of plants that make up their diet, that there is no direct,

or at most a weak but negative, correlation between their community diversities.

However, we found that microbiome diversity positively correlates with both

consumer community diversity and diet diversity in three target species. Moreover,

these correlations were affected by various environmental variables, such as

anthropogenic activity, habitat island size, and a probably important nutrient source,

guano runoff from cave entrances. Our results demonstrate that traditionally ignored

food web layers, particularly the microbiome, may mask correlations between

different communities in a food web. Hence, we advocate that microbiome inventories

are routinely added to any 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 interactions between consumers and their resources. We anticipate

our approach to be useful for ecologists and environmentalist studying different

communities in a local food web.

Keywords

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, with each species

limited by a single resource, predicts that more available resources allow for more

same-guild species in a local community (Tilman and Pacala 1993). More specifically,

the diversity of the consumer community is expected to be strongly correlated to the

diversity of the communities of their resources (Hutchinson 1959, MacArthur 1965).

Many experimental and observational studies have confirmed this correlation

indeed to be positive, 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).

However, in nature it is often not clear whether such correlations are causal (i.e.

direct), or instead the result of two trophic communities responding to an environmental

variable the same way. Spatial and temporal heterogeneous abiotic factors are

possible drivers of diversity, with variation allowing more species to co-occur than in

a homogeneous environment (Tilman and Pacala 1993). For example, Hawkins and

Porter (2003) found plant and butterfly diversities to be positively correlated, but

variation in primary productivity and topographical variability to be the likely

cause. Longmuir et al. (2007) found no correlations between the communities of

producers, consumers, and bacteria, but instead associations with environmental

variables.

Similarly, a third trophic 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). Another important candidate trophic

level of importance, the microbiome, recently gained much attention as an important

companion to its host, and it was even argued that the ‘holobiont’ is the actual unit of

selection (Bordenstein and Theis 2015). Various studies showed how the microbiome

can adapt to the host’s diet and become an important co-driver of diet choice (Muegge

et al. 2011, Colman et al. 2012, Youngblut et al. 2019), thus directly influencing

consum-er-resource interactions.

We aimed to answer the question of how seemingly neutrally assembled

communities of consumers might be influenced by their plant diet in the light of this

third trophic community, the microbiome. To this end, we studied species-rich

communities of land snails (Gastropoda) living in an archipelago-like environment of

habitat islands of limestone bedrock in lowland rainforest in Sabah, Malaysian

Borneo (Schilthuizen et al. 2003a, Hendriks et al. 2019a), focusing on three common

‘target species’. These snail communities conform to the well-known mode of few

common and many rare species, and their assembly was argued to be neutral

(Schilthuizen 2011), or possibly partly driven by rare immigration events (Hendriks et

al. 2019a). Hence, we expected no direct correlations between these communities of

consumers and their plant diet. This, together with the relative ease with which these

live snails and their empty shells (a proxy for the community) can be sampled, makes

these communities an ideal study system to answer our research question.

We applied a combination of census data and metabarcoding of multiple genetic

markers to reconstruct community richness and diversity, and to study correlations

between communities. Because these correlations can be affected by the environment

(Tilman and Pacala 1993, Hawkins and Porter 2003, Longmuir et al. 2007), we also

tested for the influence of environmental variables expected to be of importance.

Based on Island Biogeography Theory (MacArthur and Wilson 1963), we expected

more dispersal and colonization (larger habitat island; less isolation; closer to a river)

to result in higher community diversities. Due to 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, due to possible

eutrophic conditions caused by runoff from bat and swiftlet guano, was expected to

have a positive influence on plant diversity and 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 and diet

community diversities. However, diversity of the third trophic level, the microbiome,

shows positive correlations with community diversity of both the consumer and

plant diet in three target species of consumers. Furthermore, while differences in

response to community change between target consumer species are often significant,

differences among locations are not. We find several environmental variables to

correlate significantly with two or three trophic communities in the same direction,

and to affect the correlations between them. Habitat island size and distance to cave

entrances have a positive influence; distance to anthropogenic activity and current

humidity have a negative influence. We conclude that traditionally ignored food web

layers, such as the microbiome, may mask correlations between communities in a

food web, and should be included in studies of consumers and their resources, along

with environmental variables.

Methods

Study system and sampling

We studied species-rich consumer communities of land snails, their diets, and

microbiomes, in the Lower Kinabatangan Floodplain in Sabah, Malaysian Borneo

(Schilthuizen 2011). Most 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. 2003a). In this study, we focussed on six different limestone

outcrops (of ~20 in the region; Figure 3.1A) and collected samples from three plots

per outcrop along its base, each plot two by two metres and at least 50 metres from

adjacent plots.

Figure 3.1 (A) Sampling locations (i.e. limestone outcrops; in black and named) in the Lower Kinabatangan Floodplain (the river in blue), Sabah, Malaysian Borneo; unsampled locations in grey. 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.

5.

45



5.

50



5.

55



longitude

la

tit

ud

e

118.0

118.1

118.2

118.3

0 5 10

Batangan

Keruak

Tomanggong Kecil

Pangi

Tandu Batu

Tomanggong 2

SE Asia

(A)

(B)

(A)

(B)

Snail communities were sampled by collecting and sorting empty snail shells

from 5 l of soil debris below each plot in 2015 and 2016, with the collection of shells

serving as proxy for the community (below referred to as ‘shell consumer community’;

Liew et al. 2008). Shells were identified to species level based on the latest taxonomic

literature (Vermeulen et al. 2015, Liew 2019a, 2019b) and counts per plot and species

were collected in a community matrix (Table S3.1A). Live snails were collected for the

study of their microbiome and diet when plots were revisited (relocated based on

GPS readings and photos) in 2017 (Table 3.1A). We focussed on three unrelated 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 from other species that we

encountered, with a maximum of 20 individuals per species per plot (Figure 3.1B).

Details on the target species and collection procedures can be found in Hendriks et

al. (2019b). An additional community matrix was constructed from these live data

(‘live consumer community’; Table S3.1B). We constructed a phylogeny for the snail

species with five commonly used barcode markers (16S, 18S, 28S, COI, H3; Webster

et al. 2012). See Methods S3.1 and Table S3.2 for details.

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; 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 (Table 3.1B). 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 rain and the level of humidity; 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).

Metabarcoding and bioinformatics

We performed metabarcoding of plant rbcL and 16S rRNA genes, to represent the diet

and microbiome, respectively, on a maximum of 20 individuals per species per plot.

DNA extraction and library preparation details are included in Methods S3.2. 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

philosophy and routine (Callahan et al. 2016) to denoise, apply quality control, and

export representative Amplicon Sequencing Variants (ASVs; for qiime2 scripts, see

Methods S3.3). Plant (rbcL) and bacterial (16S rRNA) origin of ASVs were confirmed by

blasting 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 datasets, respectively. Additionally, 16S rRNA data were

checked not to be of host (gastropod) origin using NCBI’s nucleotide BLAST search

in Geneious v9.1.6 (https://www.geneious.com). For both markers qiime2 ASV-

alignments were checked by eye, and non-aligning reads removed and double-checked

to be of non-target 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, diet, and

microbiome. In the last two, the ‘community’ is represented by the collection of ASVs

from the diet/microbiome inside the individual snail gut. 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) for each

of the three target species separately and for all non-target species lumped together.

We pooled data by species and plot (excluding species-plot combinations with ≤ 1

individual) and plotted results with 95% confidence levels by species using ‘ggplot2’

(Wickham 2016).

We studied the influence of the shell consumer community and ‘snail species’

(fixed effects) on the diet and microbiome (response variables) using generalized

linear mixed models (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

pseu-doreplication. This resulted in a total of six models: three models each for diet and

microbiome, using as data for the response variable and shell consumer community

either Shannon diversity, PD, or Chao1 richness. 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, selecting the best

Table 3.1 (A) Summary of sample sizes for which successful metabarcoding data were obtained for microbiome (n = 818) and diet (n = 645), or either (n = 820), sorted by species, location, and plot. The three target species are printed in bold. See Tables S3.1 and S3.3 for sample details. (B) Details on plot locations, community summary metrics, and environmental data as used in linear and path models.

(A) Location Batangan Keruak

Plot 2 5 6 1 3 6 7

Family Species

Ariophantidae Everettia sp.

Microcystina appendiculata (Von Moellendorff, 1893) Macrochlamys tersa (Issel, 1874)

Assimineidae Acmella cyrtoglyphe Vermeulen,

Liew & Schilthuizen, 2015 Acmella striata Vermeulen, Liew & Schilthuizen, 2015

Cyclophoridae Alycaeus jagori Von Martens, 1859 1

Chamalycaeus sp.

Japonia kinabaluensis (E.A. Smith, 1895) 1 2

Japonia sp. 2 1 2

Leptopoma pellucidum (Grateloup, 1840)

Leptopoma sericatum (Pfeiffer, 1851) 9 3 1

Pterocyclos/Opisthoporus sp. 1

Diplommatinidae Diplommatina asynaimos Vermeulen, 1993 1

Diplommatina calvula Vermeulen, 1993 2 1

Diplommatina gomantongensis

(E. A. Smith, 1894) 1

Diplommatina rubicunda (Von Martens, 1864) 1

Plectostoma concinnum (Fulton, 1901) 20 20 20 20 20 20 20

Plectostoma simplex (Fulton, 1901)

Euconulidae Kaliella accepta (Smith, 1895) 2 1

Kaliella barrakporensis (Pfeiffer, 1852) Kaliella calculosa (Gould, 1852)

Kaliella scandens (Cox, 1872) 1

Helicinidae Sulfurina martensi (Issel, 1874) 6 2 4 2

Hydrocenidae Georissa kinabatanganensis Khalik, Hendriks,

Vermeulen & Schilthuizen, 2018 1 10

Georissa similis E. A. Smith, 1894 s.l.1 _{1 9 12 20 4}

Rathouisiidae Atopos sp.

Trochomorphidae Videna metcalfei (Pfeiffer, 1845) 1 1 1

Videna sp.

Totals 39 29 36 36 53 24 30

1 _{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.

Pangi Tandu Batu Tomanggong 2 Tomanggong Kecil 2 5 6 3 6 7 1 4 5 1 3 6 Totals 2 2 4 1 1 2 1 1 2 4 3 3 6 5 1 12 20 21 15 20 20 2 20 20 18 6 163 1 1 2 3 5 1 1 1 1 15 1 1 4 1 8 1 3 1 1 1 3 6 1 20 20 20 11 22 18 20 20 20 20 20 352 3 11 14 1 1 1 6 2 2 1 1 2 3 1 15 11 20 16 16 20 13 5 20 19 2 177 1 1 1 4 1 1 44 58 39 31 65 47 45 35 46 70 61 32 820

Table 3.1 Continued.

(B) Location Batangan Keruak

Plot 2 51 ₆1 _{1 3 6 7}

Plot locations Latitude [°] 5.459 5.459 5.459 5.522 5.523 5.522 5.523

Longitude [°] 118.102 118.101 118.100 118.285 118.285 118.286 118.286

Consumer community summary metrics

Shell consumer community Chao1 richness 20.8 27.9 27.9 26.0 16.0 13.5 17.8

Shannon diversity 2.165 1.559 1.559 1.583 1.432 1.303 2.133

Simpson diversity 0.799 0.623 0.623 0.688 0.682 0.562 0.807

Shannon evenness 0.163 0.197 0.197 0.207 0.272 0.309 0.178

Faith's PD 1.845 2.686 2.686 2.270 1.568 1.345 1.484

Live consumer community Chao1 richness 5.0 13.0 4.0 12.0 5.0 3.0 10.0

Shannon diversity 0.909 0.855 0.810 0.626 1.183 0.356 1.101

Simpson diversity 0.448 0.347 0.418 0.307 0.663 0.163 0.492

Shannon evenness 0.683 0.563 0.89 0.892 0.525 2.559 0.467

Faith's PD 0.777 1.217 0.737 0.945 0.828 0.444 0.997

Environmental data

Latent variable (LV) Measurements

Habitat island size Limestone outcrop size [m2_{] 0.444 0.444 0.444 0.345 0.345 0.345 0.345}

Cave distance Distance to cave [m] 9 74 110 138 60 178 126

Next habitat island distance Distance to next outcrop [m] 2826 2941 3022 693 785 801 901

River distance Distance to river [m] 396 562 584 1020 981 852 798

Altitude [m] 40.17 37.76 33.92 26.55 14.57 7.93 18.57

Anthropogenic distance Distance to main road [m] 8310 8270 8310 1340 1440 1380 1500

Distance to plantation [m] 823 669 623 697 679 655 652

Current humidity Time since rain [hours] 9.0 7.0 5.0 0.0 0.0 2.5 5.0

Humidity level 1 1 1 4 4 4 2

1 _{Shell community data for Batangan plots 5 and 6 were pooled during fieldwork; here we have taken the same} data for either plot.

Pangi Tandu Batu Tomanggong 2 Tomanggong Kecil 2 5 6 3 6 7 1 4 5 1 3 6 5.532 5.533 5.538 5.604 5.605 5.604 5.524 5.531 5.525 5.507 5.509 5.510 118.306 118.304 118.306 118.339 118.343 118.342 118.302 118.304 118.307 118.297 118.299 118.300 33.3 13.0 16.5 36.3 25.0 30.2 22.5 14.0 26.0 26.6 27.0 26.0 2.365 1.563 1.332 1.588 1.295 1.253 2.301 0.964 1.670 1.522 2.484 2.059 0.832 0.708 0.507 0.652 0.548 0.525 0.864 0.366 0.689 0.619 0.876 0.745 0.133 0.278 0.285 0.191 0.25 0.245 0.16 0.417 0.203 0.202 0.13 0.153 2.517 1.104 1.467 2.431 1.966 2.383 1.665 1.374 2.042 2.475 2.375 2.479 4.0 4.0 8.0 4.0 9.0 17.0 8.0 3.0 6.0 9.5 6.5 5.0 0.906 1.061 0.709 0.591 1.516 0.946 1.014 0.724 1.002 1.368 1.248 0.671 0.553 0.625 0.386 0.280 0.715 0.511 0.541 0.442 0.557 0.700 0.665 0.307 0.796 0.68 0.876 1.22 0.317 0.543 0.507 1.258 0.62 0.351 0.447 0.927 0.718 0.761 0.667 0.560 0.943 0.913 0.905 0.568 0.866 1.136 0.915 0.843 0.748 0.748 0.748 0.129 0.129 0.129 0.347 0.347 0.347 0.025 0.025 0.025 315 476 385 8077 8078 8269 996 513 554 2167 2163 2151 265 296 817 2169 2603 2471 318 151 384 1359 1459 1571 762 561 136 222 604 562 496 888 902 59 22 66 31.25 14.49 17.56 11.04 16.22 8.07 23.59 48.72 39.63 14.49 20.31 25.08 1590 1570 2060 10300 10600 10500 716 1380 1130 264 208 211 1120 1260 954 1180 1470 1450 753 1110 489 264 232 252 6.5 5.5 3.5 19.0 16.0 12.0 2.5 8.5 5.0 28.5 30.5 33.0 2 1 2 1 1 1 3 1 1 2 1 1

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 between the

target species (adjusted for multiple comparisons based on Tukey’s method) to

variations in the shell consumer community, 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 due to 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 Modelling (Sanchez 2013) using the r package

‘plspm’ v0.4.9 (Sanchez et al. 2017). We started with a ‘core model’ with three latent

variables (LVs) corresponding to the three communities: consumer, diet, and

microbiome (Figure 3.2). LVs were constructed from combinations of summary

metrics Shannon and Simpson diversity, PD, and Shannon evenness (Table S3.3).

Model assessment followed Sanchez (2013). In short, we checked unidimensionality,

cross-loadings, and positive loadings of the outer model; negative loadings were

encountered for Shannon evenness and resolved by taking negative values for this

metric. The inner model was checked for communality, redundancy, and overall

goodness-of-fit. We performed bootstrapping to assess significance by rerunning

models 999 times with data from 300 random individuals (no replacement) from any

species for which we had both diet and microbiome data (‘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 (‘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, 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). 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 datasets.

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-4353; Table S3.1A); in addition,

we collected 1,494 live snails from 28 species, with a mean of 62 individuals per

plot (range 33-139; Table S3.1B). From these live snails, we selected individuals for

metabarcoding (aimed at 20 individuals/species/plot; total ca. 840 individuals),

and obtained successful metabarcoding results for 820 individuals (Table 3.1).

Metabarcoding was less successful for the plant diet (data for 645 individuals; for

other samples data were deficient) than for the gut microbiome (data for 818

individuals; see Table S3.3 for details), which was probably a result of less plant than

microbial DNA from extractions. Mean Chao1 richness (based on ASVs) for diet and

microbiome were 3.32 and 74.3, respectively (Table S3.3). Mean values for the other

metrics were as follows: Shannon (0.88 vs. 2.9) and Simpson diversity (0.47 vs. 0.86),

Faith’s PD (0.45 vs. 15.20), and Shannon evenness (4.85 vs. 0.10).

Consumer-diet correlations

Model selection by GLMMs showed no significant effect of the shell consumer

community on the individual diet Shannon diversity for the three target species

(Table S3.4). Alycaeus jagori individuals had the most diverse diets, with individuals

from the other two target species having significantly lower diet diversities (Table 3.2,

Figure S3.2A). Model selection and best models of GLMMs for Faith’s PD and Chao1

richness showed a significant positive influence of the shell consumer community on

the individual diet diversity of A. jagori only (Table 3.3, Figure S3.2B-C). Differences

among several locations were significant for A. jagori (Chao1 richness) and P. concinnum

(Shannon diversity), with no clear geographical pattern (Tables S3.5-S3.7, Figure

S3.3A-C). There was no clear pattern among the locations in differences between

target species’ responses (Table S3.7).

Path models on shell consumer community showed no significant correlation

between consumer community and individual diet diversities (Figure 3.2A, Table

S3.8A), for both the complete model (including all species) and the normalized models

(target species only; Table S3.8A). The full normalized path model based on live

consumer community suggested a significant negative correlation instead (Figure

3.2B, Table S3.8B). Ordination of the diet showed much overlap between target species

and different locations, indicating little diet differentiation (Figure 3.3A).

Consumer-microbiome correlations

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

Ta bl e 3 .2 G LM M b es t m od el r es ul ts , u si ng f un ct io n ‘ gl m m T M B’ f ro m R p ac ka ge ‘ gl m m T M B’ v 0. 2. 3 ( Br oo ks e t a l. 2 01 7) . S ig ni fic an t t er m s i n b ol d. Fo r m od el s el ec ti on , s ee T ab le S 3. 4. N ot e t ha t S ha nn on di ve rs it y a nd F ai th ’s p hy lo ge ne ti c d iv er si ty c ou ld n ot b e c al cu la te d w he n C ha o1 e qu al s o ne , he nc e s m al le r s am pl e s iz es f or t he se m et ri cs i n d ie t d at a. R es po ns e Me tr ic n Co effi ci ent Es ti m ate SE z-v al ue p D ie t D iv er sit y ( Sh an no n) 1 35 7 In te rc ept 0.9 8 0. 04 25 .9 3 <0 .0 01 Sp ec ie s ( G . s im ili s s .l.) -0. 27 0.0 7 -3 .9 1 <0 .0 01 Sp ec ie s ( P. c on ci nn um ) -0. 20 0. 05 -4 .0 8 <0 .0 01 Ph yl og en et ic d iv er sit y ( Fa it h' s P D ) 2 35 7 In te rc ept -1. 77 0.41 -4 .31 <0 .0 01 Sh el l c on su m er c om m un it y p hy lo ge ne ti c d iv er si ty ( PD ) 0. 56 0.18 3.0 6 0.0 02 Sp ec ie s ( G . s im ili s s .l.) 0. 74 0. 81 0.9 1 0. 361 Sp ec ie s ( P. co nc in num ) 0. 84 0.4 5 1.8 6 0. 063 Sh el l c on su m er c om m un it y p hy lo ge ne tic d iv er si ty ( PD ) * S pe ci es ( G . s im ili s s .l.) -0. 85 0. 41 -2 .0 6 0. 03 9 Sh el l c on su m er c om m un ity p hy lo ge ne tic d iv er si ty ( PD ) * S pe ci es ( P. c on ci nn um ) -0. 64 0. 21 -3 .0 6 0.0 02 Ric hn es s ( Ch ao 1) 2 53 9 In te rc ept 0. 75 0. 29 2. 57 0. 01 0 Sh el l c on su m er c om m un it y r ic hn es s ( Ch ao 1) 0. 03 0. 01 2. 78 0. 005 Sp ec ie s ( G . s im ili s s .l.) -0 .31 0.4 7 -0 .6 6 0. 51 2 Sp ec ie s ( P. co nc in num ) 0. 57 0. 31 1.8 6 0. 062 Sh el l c on su m er c om m un it y r ic hn es s ( Ch ao 1) * S pe ci es ( G . s im ili s s .l.) -0 .02 0. 02 -1. 20 0. 23 0 Sh el l c on su m er c om m un it y r ic hn es s ( Ch ao 1) * S pe ci es ( P. c on ci nn um ) -0. 05 0. 01 -4 .13 <0 .0 01 M ic ro bi om e D iv er sit y ( Sh an no n) 2 69 0 In te rc ept 0.9 8 0. 14 7.0 1 <0 .0 01 Sh el l c on su m er c om m un it y d iv er sit y ( Sh an no n) 0.1 3 0. 08 1.5 8 0. 114 Sp ec ie s ( G . s im ili s s .l.) 0. 01 0. 11 0. 11 0.9 15 Sp ec ie s ( P. c on ci nn um ) -0 .21 0. 10 -2 .10 0. 03 6 Sh el l c on su m er c om m un it y d iv er si ty ( Sh an no n) * S pe ci es ( G . s im ili s s .l.) -0 .12 0.0 6 -2 .13 0. 033 Sh el l c on su m er c om m un it y d iv er sit y ( Sh an no n) * S pe ci es ( P. co nc in num ) 0. 01 0. 06 0.1 3 0. 89 9 Ph yl og en et ic d iv er sit y ( Fa it h' s P D ) 3 69 0 In te rc ept 26 .6 6 4.6 4 5. 75 <0 .0 01 Sh el l c on su m er c om m un it y p hy lo ge ne ti c d iv er si ty ( PD ) -5 .59 2.2 0 -2 .5 4 0. 011 Sp ec ie s ( G . s im ili s s .l.) 6. 67 3. 61 1.85 0. 064 Sp ec ie s ( P. c on ci nn um ) -10 .5 6 2.9 0 -3 .6 4 <0 .0 01 Sh el l c on su m er c om m un it y p hy lo gen et ic d iv er sit y ( PD ) * S pe ci es ( G . s im ili s s .l.) -1. 74 1.67 -1. 04 0. 29 7 Sh el l c on su m er c om m un ity p hy lo ge ne tic d iv er si ty ( PD ) * S pe ci es ( P. c on ci nn um ) 4.1 5 1. 36 3.0 6 0.0 02 Ric hn es s ( Ch ao 1) 3 69 0 In te rc ept 12 6. 52 29 .74 4. 25 <0 .0 01 Sh el l c on su m er c om m un it y r ic hn es s ( Ch ao 1) -1. 86 1.18 -1. 57 0. 116 Sp ec ie s ( G . s im ili s s .l.) 45 .0 4 23 .69 1.9 0 0. 057 Sp ec ie s ( P. c on ci nn um ) -5 0. 29 18 .13 -2 .7 7 0. 00 6 Sh el l c on su m er c om m un it y r ic hn es s ( Ch ao 1) * S pe ci es ( G . s im ili s s .l.) -1 .8 4 0.9 3 -1 .9 8 0.0 48 Sh el l c on su m er c om m un it y r ic hn es s ( Ch ao 1) * S pe ci es ( P. co nc in num ) 1.2 4 0.6 8 1.82 0. 069 1 R es po ns e v ar ia bl e m od el le d a s a n or m al d is tr ib it io n. 2 R es po ns e v ar ia bl e m od el le d a s a l og no rm al d is tr ib it io n. 3 R es po ns e v ar ia bl e m od el le d a s a g am m a d is tr ib it io n.

species (Table 3.2, Table S3.4). Individual microbiome Shannon diversity increased

with shell consumer community, though not significantly (Table 3.2, Figure S3.2D).

In contrast, Faith’s PD and Chao1 showed a negative trend, the former significant for

A. jagori and G. similis s.l. and Chao1 for G. similis (Table 3.2, Table 3.3, Figure S3.2E-F).

Significant differences among several locations were found in Shannon diversity for

A. jagori (Tables S3.5-S3.7, Figure S3.3D-F). There was no clear pattern among the

locations in differences between target species’ responses (Table S3.7).

In contrast to the GLMMs, core path models on shell consumer community

diversity showed a significant positive correlation between consumer and microbiome

diversities, for both complete (all species; Figure 3.2A) and normalized models (target

species only; Table S3.8A). In full path models and models based on live consumer

community (Table S3.8B), this correlation was non-significant. Ordination showed

each target species to occupy a subset of the microbiome niche space, but with overlap

between A. jagori and P. concinnum (Figure 3.3B).

Table 3.3 Results from post hoc tests of trends in species responses to changes in shell consumer community. Results based on full models (i.e. including interactions between fixed effects ‘species’ and ‘shell consumer community’), using function ‘emtrends’ from the R package ‘emmeans’ (Lenth 2019). Trends printed in bold were significant (p < 0.05).

Response Metric Species Trend SE df p

Diet Diversity (Shannon) A. jagori 0.106 0.086 349 0.219

G. similis s.l. 0.082 0.124 349 0.510

P. concinnum -0.105 0.079 349 0.182

Phylogenetic diversity (Faith's PD) A. jagori 0.562 0.183 349 0.002

G. similis s.l. -0.287 0.391 349 0.463

P. concinnum -0.077 0.148 349 0.604

Richness (Chao1) A. jagori 0.032 0.011 531 0.006

G. similis s.l. 0.009 0.019 531 0.642

P. concinnum -0.021 0.011 531 0.059

Microbiome Diversity (Shannon) A. jagori 0.125 0.079 682 0.115

G. similis s.l. 0.006 0.079 682 0.939

P. concinnum 0.132 0.075 682 0.079

Phylogenetic diversity (Faith's PD) A. jagori -5.590 2.205 682 0.011

G. similis s.l. -7.326 2.129 682 0.001

P. concinnum -1.437 1.882 682 0.445

Richness (Chao1) A. jagori -1.860 1.185 682 0.117

G. similis s.l. -3.697 1.196 682 0.002

Microbiome-diet correlations

All path models showed a significant positive correlation between individual micro -

biome and diet diversities (Figure 3.2, Table S3.8). 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.

Environmental effects

Full path 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 the shell consumer community; Table S3.8A), and thus explained more of the

variation in the latent variables (LVs). (Note that these goodness-of-fit values are only

indicative and not appropriate for detailed model comparison, cf. Sanchez 2013.)

Figure 3.2 Significant results from Partial Least Squares Path Modelling (PLS-PM; Sanchez 2013) for normalized models (i.e. with equal sample sizes from the three 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).

MICROBIOME

DIET

habitat island size current humidity anthropogenic distance cave distance river distance

CONSUMER

next habitat island distance 1.33 0.55 0.87 0.34 0.11 -2.05 1.33 0.60 0.32 -0.63 -0.24 -0.23 -0.27

MICROBIOME

DIET

habitat island size current humidity anthropogenic distance cave distance

CONSUMER

next habitat island distance 0.52 0.39 0.11 -0.38 0.41 0.80 -0.71 -0.14 -0.22 -0.12 -0.44 0.42

(A)

(B)

0.16 0.18 0.20

(A)

(B)

In the full models, correlations between core LVs (i.e. consumer, diet, and microbiome)

were lower than in core models, and part of the explanatory power was ‘moved’ to

the environmental LVs (Figure 3.2). In the next paragraph, we discuss results for

normalized models (i.e. based on equal sample sizes of target species, for which we

had most data), unless stated otherwise. Results for complete models (based on all

species) were in the same direction, but path coefficients were less often significant

(Table S3.8).

Several environmental variables were significantly correlated with community

LVs (Figure 3.2). We found in general that (i) diversities of consumer, diet, and

microbiome communities were higher close to human activity, (ii) larger habitat

islands supported more diverse snail consumer communities with a more diverse

diet and microbiome, and (iii) snail consumer and diet diversities 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

diversities, more isolated locations had higher snail consumer and lower individual

microbiome diversity, and more diverse individual diets were found further away

from the river.

Discussion

We studied the correlations between 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’ plant diets, and gut micro -

biomes. We aimed to answer the question of how seemingly neutrally assembled

communities of consumers might be influenced by their plant diet in the light of a

third trophic community, the microbiome. We found no (shell consumer community),

or at most a weak but negative (live consumer community), direct correlation between

consumer community and individual diet. However, we found both the 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 weak negative) correlation between consumer and individual

diet diversities contradicts 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). 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. In a study very similar

to ours, also dealing with three trophic communities (pelagic zooplankton,

phyto-plankton, and bacteria), no significant correlations between community diversities

were detected either (Longmuir et al. 2007). Several other large-scale studies on

aquatic ecosystems have also reported the absence of community correlations (Allen

Figure 3.3 nMDS ordination plots based on Bray-Curtis 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 non-target 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 were available (singletons) were excluded. Numbers refer to plot identity with each location, for details of which, see Table 3.1. Ellipses indicate 95% confidence levels. 2 2 5 5 6 6 6 1 1 1 3 3 3 6 6 7 7 2 2 2 5 5 5 5 6 3 3 6 6 6 6 7 7 7 1 1 4 5 5 1 1 1 ₃ 3 3 3 6 6 6 −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0.4 −0.8 −0.6 −0.4 −0.2 0.0 0.2 nMDS1 nM DS 2 stress: 0.17 2 5 2 5 6 6 6 1 1 1 3 3 3 6 6 72 2 2 5 5 5 5 6 6 3 3 6 6 6 6 7 7 1 1 1 4 4 4 5 5 1 1 1 1 3 3 6 6 6 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 −0.8 −0.6 −0.4 −0.2 0.0 0.2 nMDS1 stress: 0.24 (A) (B) 2 2 5 5 6 6 6 1 1 1 3 3 3 6 6 7 7 2 2 2 5 5 5 5 6 3 3 6 6 6 6 7 7 7 1 1 4 5 5 1 1 1 ₃ 3 3 3 6 6 6 −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0.4 −0.8 −0.6 −0.4 −0.2 0.0 0.2 nMDS1 nM DS 2 stress: 0.17 2 5 2 5 6 6 6 1 1 1 3 3 3 6 6 72 2 2 5 5 5 5 6 6 3 3 6 6 6 6 7 7 1 1 1 4 4 4 5 5 1 1 1 1 3 3 6 6 6 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 −0.8 −0.6 −0.4 −0.2 0.0 0.2 nMDS1 stress: 0.24 (A) (B) 2 2 5 5 6 6 6 1 1 1 3 3 3 6 6 7 7 2 2 2 5 5 5 5 6 3 3 6 6 6 6 7 7 7 1 1 4 5 5 1 1 1 ₃ 3 3 3 6 6 6 −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0.4 −0.8 −0.6 −0.4 −0.2 0.0 0.2 nMDS1 stress: 0.17 2 5 2 5 6 6 6 1 1 1 3 3 3 6 6 72 2 2 5 5 5 5 6 6 3 3 6 6 6 6 7 7 1 1 14 4 5 5 1 1 1 1 3 3 6 6 6 −0.6 −0.4 −0.2 0.0 0.2 nMDS1 Location Batangan Keruak Pangi Tandu Batu Tomanggong 2 Tomanggong Kecil Species A. jagori G .similis s.l. non−target species P. concinnum (B)

(A)

(B)

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 biotic diversity, but instead more by “available microclimatic

and microchemical gradients.” This was also found by Hawkins and Porter (2003).

The butterfly and plant diversities they studied were both influenced by the same

environmental variables, namely primary productivity and topographical variability.

The weak negative correlation we found is however intriguing, because it illustrates

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

niche widths, which is the inverse of ecological release, in which case populations

living in a low competitive community are able to consume a wide array of resources

(Kernaléguen et al. 2015).

Although we found no (or weak) direct correlations between consumer

community and individual diet, we did find significant positive correlations for the

diversity of these two communities with the third community, the microbiome,

suggesting a tritrophic community system. The microbiome is partly composed of

microorganisms inherited from the parents to the offspring (vertical transfer), but

also of bacteria obtained from food sources and community members (horizontal

transfer), with 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, faecal transplant studies (Kohl et al. 2014, Kohl and Dearing 2016) and

field studies (Kohl et al. 2018) show 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 (Hird et al. 2018,

Glasl 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

co-existing species. This could allow individuals to feed on a more diverse set of

plants, or potentially get other benefits, such as better immunity. This does not

explain why we did not find a direct consumer-diet correlation, however, and we

suggest further experimental research in a more controlled environment (without

potentially confounding effects) to explain our finding.

Correlations between communities and several environmental variables (mainly

anthropogenic distance, habitat island size, and presence of cave entrances) are

stronger than between communities themselves, and it is possible that correlations

we found between communities are actually the result of similar responses to the

environment (Allen et al. 1999). 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 (Stadler et al. 2000, Fine 2002).

Similar data for snail communities are rare, but snail dispersal via crops is common

(Dörge et al. 1999). Although richness of snail communities (in non-limestone habitat)

was 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

were previously shown to perform surprisingly well after (human) disturbance

(Schilthuizen et al. 2005b). Besides, anthropogenic activity may increase dispersal

between limestone outcrops and in that way increase local species diversity (Cadotte

2006). Second, in agreement with MacArthur and Wilson’s Theory of island

biogeography (MacArthur and Wilson, 1963), community diversities we studied

increase with habitat island size. However, contrary to this theory, consumer

community and individual diet diversities do not increase with decreasing distance

to the next habitat island. This is in line with recent work by Hendriks et al. (2019a),

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 to 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), though.

We presented the diet and microbiome diversity captured from both the entire

consumer community (complete model) and from normalized data (including only

the three target species; normalized model). Above, 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,

as we found for instance different results for the different target species in the

outcome of the GLMMs. However, path model correlations between consumer,

individual diet, and individual microbiome did not differ much between the complete

and normalized models, indicating that our findings are rather robust. The influence

of the environmental variables was less significant in the complete model, however.

Our results demonstrate that traditionally ignored food web layers, such as the

microbiome, and environmental variables, may mask correlations between different

food web communities. 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

environmen-talists.

Data accessibility

All samples (both empty shells and fresh snails) were deposited in the BORNEENSIS

collection of the Institute for Tropical Biology and Conservation, Universiti Malaysia

Sabah, Kota Kinabalu, Malaysia (BORN) or the molluscan collection of the Naturalis

Biodiversity Center, Leiden, the Netherlands (RMNH). All DNA extractions were

stored at the Naturalis Biodiversity Center at -80°C for future reference. Museum IDs

for samples for which DNA extractions were performed are given in Tables S3.2 and

S3.3. Barcode marker read data for the reconstruction of the snail consumer

community phylogeny were deposited to the Barcode of Life Database (www.

boldsystems.org) as dataset ‘DS-2019PHY’ and NCBI GenBank; see Table S3.2 for

accession details. Metabarcoding data were deposited to the NCBI GenBank Sequence

Read Archive (SRA) database as project PRJNA530120 and can be retrieved from

https://www.ncbi.nlm.nih.gov/sra/PRJNA530120 upon publication; see Table S3.3 for

per sample accession details. A complete R-project zip-file with raw data and scripts

to recreate all models presented in this paper, as well as figures and tables, will be

deposited to Dataverse, Dryad, or similar upon publication.

Acknowledgements

We thank Liew Thor-Seng of ITBC, UMS, Kota Kinabalu, Malaysia, for kindly assisting

with access and export license applications, Iva Njunjić, Leonel Herrera-Alsina, and

Giacomo Alciatore for assisting with fieldwork, and UMS students Gary Kam Yin

Cheng, Foo She Fui, Choo Ming Huei, Chua Wai Min, James Yee Chun Sieng, Phung

Kin Wah, and Chee Huey Ying for identifying snail shells in the lab. Groningen

students Annabel Slettenhaar and Manon Spaans collected environmental data, and

Marti J. Anderson gave advice about statistical analyses. This research was funded

by NWO (VICI grant, number 865.13.003, awarded to R.S.E.), KNAW (Fonds Ecologie,

reference Eco/1711, awarded to K.P.H.), Leopold III-Fonds (awarded to K.B., 2017), Ghent

University (Special Research Fund, BOF, awarded to K.B.), University of Groningen

(Ubbo Emmius sandwich program, awarded to K.B.) and The Malacological Society of

London (Early Career Research Grant, awarded to K.P.H.). All samples were collected

(permit numbers JKM/MBS.1000-2/2 JLD.6 (107, 112, 114, 116, and 118)) and exported

(JKM/MBS.1000-2/3 JLD.3 (51)) under license of Sabah Biodiversity Council (SaBC).

Supplementary material

Methods S3.1 Consumer community phylogeny reconstruction.

A phylogeny for the snail community was reconstructed in order to calculate Faith’s phylogenetic diversity (Faith’s PD; Faith 1992). We used five commonly used barcode markers (16S, 18S, 18S, COI, H3; sensu Webster et al. 2012). For each species, a representative was chosen at random from our own collections made during fieldwork in the region in 2015 and 2016 (Table S3.2). When no own data were available for a species (because only shells were found, and no live animals, or laboratory procedures failed), data were retrieved from GenBank; in case data for the species were not available from GenBank, data for the most closely related species were taken (preferably same genus).

DNA extractions, primers, and amplification procedures followed Webster et al. (2012). Successful amplification products were sent to BaseClear, Leiden, the Netherlands, for Sanger sequencing in two directions. Sequence reads were checked for errors using Geneious v9.1.6 (https://www.geneious.com). Resultant sequence read data were uploaded to the Barcode of Life Database (www.boldsystems.org) and subsequently submitted to GenBank. Read data can be found from BOLD, dataset ‘DS-2019PHY’, or GenBank (accession numbers given in Table S3.2).

The phylogeny was reconstructed using beast2 v2.4.7 (Bouckaert et al. 2014), generally following the procedure described by Hendriks et al. (2019a). In short, we linked gene trees for each locus and applied a strict clock. We chose a GTR site model for each locus. We set a Yule tree prior with default clock rate, as we were only interested in relative branch lengths among species to calculate Faith’s PD. We ran the analysis for 100 million generations, sampling posterior parameter values and trees every 10,000th _{generation, after which we}

discarded a 10% burn-in. We checked convergence based on ESS values of > 200 and proper mixing of parameters over time. To obtain a final species tree (Figure S3.1) for calculation of Faith’s PD, we summarized trees with a posterior probability consensus limit of 50%.

Methods S3.2 Library preparation for metabarcoding of diet

and microbiome.

For metabarcoding of diet and microbiome, we performed whole-DNA extractions on gut contents of individual snails, or –in the case of minute species (G. similis s.l. and equally-sized species)– whole snails after removing the shell, using Omega’s E.Z.N.A.® _{Mollusc DNA Kit.}

Library preparation for the metabarcoding of microbiome and diet involved two amplification steps. A first PCR was performed with a study-related primer set and a second amplification to add multiplex Nextera XT (Illumina) indexes to each sample, so that each sample could be uniquely identified during bioinformatic analyses after sequencing.

For the diet study the first amplification was of the 110 bp chloroplast rbcL region using primers Z1af and 19bR from Hofreiter et al. (2000), with an Illumina adapter overhang to the primers. The 20 μl PCR reaction contained 1 μl 10pM each of forward primer and reverse primer, 0.5 μl 5U/μl Phire Hot Start II polymerase (Finnzymes), 5 μl PCR buffer (FinnZymes), 5 μl 2.5mM dNTPs, 1 μl 100mM BSA, and 2 μl of template. The amplification reaction was performed with the following cycle conditions: Initial denaturation 30 seconds at 98°C, followed by 40 cycles of denaturation of 10 seconds at 98°C, 10 seconds annealing at 50°C, an extension of 20 seconds at 72°C, and a final extension of 5 minutes at 72°C. We included one negative control per PCR plate.

For the microbiome study the first amplification was of the 426 bp (in E. coli) bacterial 16S rRNA V3-V4 region with primers 343F (Liu et al. 2007) and 748R (Andersson et al. 2008). These primers contained an Illumina adapter overhang. The 20 μl PCR reaction contained 10 μl 2X TaqMan™ Environmental Master Mix 2.0 (ThermoFischer Scientific), 1 μl of 10pM each of forward and reverse primer, and 2 μl of template. The amplification reaction was performed with the following cycle conditions: Initial denaturation and activation of the polymerase of 10 minutes at 95°C, followed by 35 cycles of 3 minutes of denaturation at 95°C, 30 seconds of annealing at 50°C, an extension of 40 seconds at 72°C, and a final extension of 5 minutes at 72°C. We included one negative control per PCR plate.

All generated amplicons were checked on an E-gel 96 2% agarose (Invitrogen), purified with 0.9x NucleoMag® _{NGS Clean-up and Size Select (Macherey-Nagel), and eluted in 25 μl}

Milli-Q water.

The second amplification step, to add the multiplex Nextera XT (Illumina) indexes to each individual amplified sample, was accomplished with a 20 μl PCR reaction containing 5 μl Milli-Q water, 10 μl 2X TaqMan™ Environmental Master Mix 2.0, 1 μl 10pM Nextera XT-S multiplex primer, 1 μl 10pM NexteraXT-N multiplex primer, and 3 μl amplified template. The amplification reaction was performed with the following cycle conditions: Initial denaturation and activation of the polymerase of 10 minutes at 95°C, followed by 10 cycles of 30 seconds of denaturation at 95°C, 1 minute of annealing at 55°C, an extension of 30 seconds at 72°C, and a final extension of 7 minutes at 72°C.

The concentration of the individually indexed amplicon samples was measured with a QIAxcel Advanced System (Qiagen), using a QIAxcel DNA Screening Kit (Qiagen), followed by normalization and pooling of the samples using a QIAgility robot (Qiagen). The quality and concentration of the final pool was checked on a Bioanalyzer using a High sensitivity chip (Agilent).

Paired-end sequencing of six pools from library preparation (three rbcL, three 16S rRNA; each pool with ~ 280 samples) was performed by BaseClear, Leiden, the Netherlands, on Illumina Miseq PE 300bp, using 50% PhiX Control v3 to counteract negative effects of low diversity libraries. Raw sequences (fastq-files) were deposited in the NCBI GenBank Sequence Read Archive (SRA) database as project PRJNA530120 and can be retrieved from https://www.ncbi.nlm.nih.gov/sra/PRJNA530120 upon publication (see Table S3.3 for accession details).

METHODS S3.3 Qiime2 scripts for the bioinformatics of raw read data.

We used Qiime2 v2017.12 (Bolyen et al. 2018) and raw, demultiplexed, Illumina MiSeq data, to obtain an Amplicon Sequence Variant (ASV) table, representative ASV sequences, and a taxonomic identification for the representative ASVs. Raw data can be downloaded as a set of fastq-files (two files per sample, times two markers, i.e. rbcL and 16S) from the GenBank Sequence Read Archive (SRA) as project PRJNA530120 (https://www.ncbi.nlm.nih.gov/sra/ PRJNA530120). Individual accession numbers are listed in Table S3.3.

We ran the following Qiime2 commands from a location with folder with the fastq-files, below named ‘Pool_16S_all’. The same code was repeated for the rbcL data.

Command 1: Import data into Qiime2.

qiime tools import \

--type ‘SampleData[PairedEndSequencesWithQuality]’ \ --input-path Pool_16S_all \

--source-format CasavaOneEightSingleLanePerSampleDirFmt \ --output-path demux-paired-end-snails1.qza

Command 2: Denoise data. Combine paired-end data; cut primers; truncate to marker length; remove chimeras.

qiime dada2 denoise-paired \

--i-demultiplexed-seqs demux-paired-end-snails1.qza \ --p-trim-left-f 14 \ --p-trunc-len-f 300 \ --p-trim-left-r 19 \ --p-trunc-len-r 300 \ --o-representative-sequences rep-seqs-snail1.qza \ --p-n-threads 8 \ --o-table table-snail1.qza

Command 3: Export results as biom-file.

qiime tools export \ table-snail1.qza \

--output-dir exported-feature-table-snail1

Command 4: Export a fasta-file with representative sequences for the ASVs, which allows building a phylogeny in package FastTree (performed outside of Qiime2; see main text).

qiime tools export \ rep-seqs-snail1.qza \

--output-dir rep-seq-table-snail1

Command 5: Blast representative sequences for the ASVs against a reference database, such as the GreenGenes database.

qiime feature-classifier classify-sklearn \ --i-classifier classifier_16S.qza \ --i-reads rep-seqs-snail1.qza \ --o-classification taxonomy-snail1.qza

Command 6: Export result from command 5, a taxonomy table.

qiime tools export \ taxonomy-snail1.qza \ --output-dir taxonomysnail1

Output files were loaded and combined into r v3.6.0 (R Core Team 2018) and combined using package ‘PhyloSeq’ v1.24.2 (McMurdie and Holmes 2013).

Classifiers (used in command 5) are databases used to identify ASVs and can either be downloaded online, but are better created and trained for each dataset, such that primer sequences and target regions in the database match the new data to be blasted. The code below shows how we created and trained a classifier for the diet study using a 110 bp region from the rbcL gene.

Command 6: Import a reference fasta-file and taxonomy-file, in this case from Bell et al. (2017).

qiime tools import \

--type ‘FeatureData[Sequence]’ \

--input-path rbcL_all_Jan2016.rdp_updated.fasta \ --output-path rbcL_all_Jan2016.rdp_updated.qza qiime tools import \

--type ‘FeatureData[Taxonomy]’ \

--source-format HeaderlessTSVTaxonomyFormat \ --input-path rbcL_Taxonomy_updated.txt \ --output-path ref-taxonomy-rbcL.qza

Command 7: Extract reads that contain the primer sequences we used in our study.

qiime feature-classifier extract-reads \ --i-sequences rbcL_otus.qza \

--p-f-primer ATGTCACCACCAACAGAGACTAAAGC \ --p-r-primer CTTCTTCAGGTGGAACTCCAG \ --p-trunc-len 110 \

--o-reads ref-seqs-rbcL.qza

Command 8: Train the classifier to be used in command 5.

qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads rbcL_all_Jan2016.rdp_updated.qza\ --i-reference-taxonomy ref-taxonomy-rbcL.qza\ --o-classifier classifier_rbcL_Bell2017.qza