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.

their plant diet, but diets differ in richness

Kasper P. Hendriks‡, Karen Bisschop‡, James C. Kavanagh,

Hylke H. Kortenbosch, Anaïs E. A. Larue, Dries Bonte,

Elza J. Duijm, Joana Falcão Salles, Francisco J. Richter Mendoza,

Menno Schilthuizen, and Rampal S. Etienne

‡ Joint First Authors

Unpublished working manuscript,

ready for submission.

Abstract

Herbivore diets are often generalistic, and communities of herbivores tend to share

much of their diets. The composition of these herbivore communities is usually well

described by predictions from the Unified Neutral Theory of Biodiversity (UNTB),

which assumes functional equivalence of community members. This is in stark

contrast to niche theory, which states that species need functional differentiation to

co-exist. We tried to answer the question whether diet differentiation plays a role in

the co-existence of tens of different herbivorous land snail species that co-exist in

communities on limestone in the tropical lowlands of Malaysian Borneo. We show,

with a large metabarcoding study of the plant diet from gut contents of 658 individual

snails, that the different snail species indeed share much of their plant diet, but mean

diet richness varies strongly among species (up to 15.3×). These differences are mostly

explained by snail size, with larger snails having richer diets, and Sørensen diet

distance between species positively correlates with difference in snail species size.

Furthermore, phylogenetic analyses of the plant diet by individual snail showed

signs of clustering in ca. 28% of the individuals, suggesting diet choice, although such

clustering was weak when diets were considered by species. We discuss how the

trends in diet richness and possible diet choice we find could also be explained from

random feeding, with larger species simply eating more or less specifically, and other,

non-competitive interactions, such as snails avoiding desiccation. Our study shows

how to efficiently put the power of metabarcoding to work in unravelling complex

community processes, such as commonly encountered in tropical ecosystems, which

should serve both community ecologists and conservationists.

Keywords

Introduction

Herbivore communities often consist of many different species feeding on roughly

the same generalist diet (Belovsky 1986, Gordon and Illius 1989). Such ecological

communities have been a popular study system for a long time (Schoener 1974),

probably because some of the most important species interactions take place within

herbivore communities: competition and facilitation (Schoener 1983, Stachowicz

2001, Bruno et al. 2003). However, the very existence of such trophically similar,

horizontal communities seems paradoxical (Behmer and Joern 2008). On the one

hand, trophically similar species have a tendency to cluster together, because they

have very similar needs (Leibold and McPeek 2006); this lies at the very heart of the

definition of the community. On the other hand, it is generally accepted that too

similar species cannot co-exist indefinitely due to too strong competition, as

formalized in the ‘competitive exclusion principle’, or Gause’s law (Gause 1934, Hardin

1960).

The assembly of trophically similar communities has been described by the

Unified Neutral Theory of Biodiversity (UNTB; Hubbell 2001, Rosindell et al. 2011).

The central axiom of the UNTB is that different species in the community are

functionally equivalent (Alonso et al. 2006). This is in stark contrast with niche

theory, which highlights the importance of functional differences in co-existing

species (Grinnell 1917). Various publications described efforts to unify both theories

in the form of so-called ‘emergent neutrality’ (Holt 2006, Segura et al. 2011, Scheffer

et al. 2018, D’Andrea et al. 2019). This could arise when underlying weak stabilizing

forces (on the species) and strong equalizing forces (increasing species similarity)

work interactively. This ‘near-neutrality’ would most likely be possible in species-rich

communities (Holt 2006). Behmer and Joern (2008) showed experimentally how

dietary and nutritional differences in co-existing herbivorous species, eating the

same plant taxa, can be extremely cryptic and thus easily overlooked, again leaving a

seemingly neutral footprint. And there are other subtle ways for species to co-exist

on the same resources, as shown long ago by Hutchinson’s seminal work on the

‘paradox of the plankton’ (Hutchinson 1961), with co-occurring species differentiating

in other niche dimensions than just food or nutrients.

Here we ask how it is possible that tens of different land snail species can co-exist

within a community, and whether there is a role played by diet differentiation among

species. Communities of large mammalian herbivores have been studied thoroughly

in this respect, and from these it is known that sympatric species can co-exist due to

(subtle) differences in morphology (digestive system, body size) and spatial and

seasonal differences in their diets (Prins and Olff 1998, Kartzinel et al. 2015). Similar

studies on invertebrates are far less numerous. We studied -seemingly neutral- land

snail communities living on limestone outcrops in the tropical lowland of Sabah,

Malaysian Borneo (Schilthuizen 2011, Hendriks et al. 2019a), and applied meta

-barcoding of the gut contents of 658 individuals from 26 species to reconstruct

individual and species-level plant diets (Pompanon et al. 2012, Taberlet et al. 2012,

Kartzinel et al. 2015). Based on the observation by Schilthuizen (2011) of the seeming

neutrality of these communities, together with an apparent excess food availability,

we expected competition for food to be very low or absent. However, these communities

are characterized by two striking features that might not be in line with neutrality

requirements, and which we therefore studied in relation to the diet. First,

communities were composed of representatives from three main taxonomic clades

(Neritimorpha, Prosobranchia, and Pulmonata), with a roughly 50:50 proportion of

Pulmonata and non-Pulmonata (Schilthuizen 2011). These taxonomic groups are

considered ecologically distinct (and thus functionally different), with Pulmonata

being better dispersers and better adapted to drought and other extreme conditions

(Schilthuizen et al. 2002, 2005a), which was expected to leave its mark on the diet.

Not much is known about the ecology of the Neritimorpha in the region, but their

preference for humid limestone (Khalik et al. 2018, 2019) suggests ecological

differen-tiation from at least the Pulmonata. Second, the range of sizes of sympatric snail

species is enormous, covering five orders of magnitude (based on adult shell volume).

Furthermore, because diet is first of all dependent on the local availability of food, we

also took into account the influence of the location. Finally, we considered different

levels of diet phylogenetic clustering, corresponding to different taxonomic ranks

(but without explicit definition of the rank). Different studies have shown that

herbivore- resource interactions can occur at either a low rank (even at the level of

genotypic variation; Barbour et al. 2015) or a higher rank (Symons and Beccaloni

1999).

We found the herbivorous land snails to have a very broad diet, including 32 seed

plant families, which was shared among most snail species. Mean diet richness

differed strongly among species: up to 15.3× (when plant taxa defined at the lowest

taxonomic rank identified). Considering all 26 snail species from our study, we found

a significant positive correlation between diet richness and snail size (although only

for Prosobranchia in analyses based on higher plant taxonomic ranks). Our study of

individual snail diets showed a trend of phylogenetic clustering of the diet in ca. 28%

of individuals. This trend was not clear when diet data were considered by target

species (i.e. diet data pooled by species for the whole region), with significant

clustering only based on diet data defined at the lowest plant taxonomic rank, and for

recent plant taxa only. Mean Sørensen diet distance among species was positively

correlated with differences in species size, although when based on mean unweighted

UniFrac diet distance this trend was positive only for the Prosobranchia, and negative

when comparing species from different classification groups.

Methods

Study system and sampling

We studied species-rich communities of herbivorous land snails on six limestone

outcrops in the Lower Kinabatangan Floodplain in Sabah, Malaysian Borneo (Figure 4.1,

Table S4.1). Many community members are known to have a strong preference for

calcium carbonate, with the outcrops effectively representing habitat islands to them

(Schilthuizen et al. 2003a), meaning that communities can be defined by such outcrops

(from here on referred to as ‘locations’). Mean location community richness was 30.3

(range: 24.0-38.0; based on data from Chapter 3). Snail samples were collected from

Figure 4.1 (A) Map of limestone outcrops in the Lower Kinabatangan Floodplain, Sabah, Malaysian Borneo. The six outcrops in black were sampled for this study. (B) Batangan, one of the outcrops visited for sampling, as seen from the river. March 2015. (C) Plectostoma concinnum (Fulton, 1901), one of the target species, feeding from limestone rock. Tandu Batu, April 2016. The scale bar equals ca. 5 mm. Photographs by Kasper P. Hendriks.

5.

45



5.

50



5.

55



longitude

la

tit

ud

e

118.0

118.1

118.2

118.3

0 5 10 km

Batangan

Keruak

Tomanggong Kecil

Pangi

Tandu Batu

Tomanggong 2

SE Asia

(A)

(B)

(C)

(A)

(B)

(C)

three plots of two by two metres per location, with an inter-plot distance of at least 50

metres, along the base of the outcrop (Hendriks et al. 2019b). All samples were

collected within one week to reduce possible seasonal effects. We focussed on three

target species, each common in the region, and not closely related to one another:

Alycaeus jagori Von Martens, 1859, Georissa similis E. A. Smith, 1893 s.l. (in fact a

complex of recently described, closely related taxa; Hendriks et al. 2019b, Khalik et al.

2019), and Plectostoma concinnum (Fulton, 1901). We aimed to collect 40 individuals/

target species/plot (20 for further analyses, plus 20 more as backup and collection

voucher), searching each of four plot quadrants for 30 minutes. In addition, we

collected all other snail species encountered during this search (usually present in

small numbers only), with a maximum of 20 individuals/species/plot. Samples were

preserved in 98% ethanol and frozen directly in the field, and registered in the

BORNEENSIS collection of Universiti Malaysia Sabah, Malaysia (Table S4.2).

Trait data were collected from various literature sources on snail taxonomy

(Vermeulen 1991, Vermeulen et al. 2015, Liew 2019a, 2019b). For each species, we noted

maximum shell height and width, and calculated from these the shell volume (as ⅓ ×

height × width

2

_{) as a proxy for snail size (Table S4.3). In addition, we calculated species}

abundances from census data reported in Chapter 3 (Table S4.4).

Metabarcoding library preparation

The process of obtaining metabarcoding plant diet data from the snail gut was

described in detail in Chapter 3, and is summarized here. We extracted genomic DNA

from each separate snail’s gut contents (or from the whole snail after removing the

shell, in case the snail was too small to prepare the gut, such as in G. similis s.l.), using

Omega’s E.Z.N.A.® Mollusc DNA Kit. Library preparation involved two amplification

steps. First, the 110 bp chloroplast rbcL region was amplified in 40 PCR cycles using

the primer pair Z1f/19bR (Hofreiter et al. 2000) with Illumina adapter overhang to

the primers. Successful products were purified with 0.9× NucleoMag® NGS Clean-up

and Size Select (Macherey-Nagel). Second, we ran another 10 PCR cycles to add sample-

specific multiplex Nextera XT indexes by which samples could be demultiplexed

during subsequent bioinformatics. We measured concentrations for each sample

using the QIAxcel Advanced System with QIAxcel DNA Screening Kit (Qiagen). We

equimolarly normalized and pooled samples into three pools of ca. 280 samples each

(including three negative controls) using a QIAgility robot (Qiagen), and checked the

quality and concentration of the final pools on a Bioanalyzer using a high sensitivity

chip (Agilent). The pools were sequenced paired-end on a MiSeq PE 300bp by

BaseClear, Leiden, the Netherlands, using 50% PhiX Control v3 to reduce effects of

low library diversity.

Metabarcoding bioinformatics

We combined raw sequence data (fastq-files) from the three sequencing pools and

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) . Taxonomic assignment of ASVs (to the lowest

rank with a confidence value of > 0.70) was performed by blasting data against our

custom-built classifier, based on the seed plant (Spermatophyta) rbcL-database from

Bell et al. (2017). Non-plant ASVs and ASVs that could not be aligned were removed.

We aligned the remaining 778 true-plant ASVs using mafft v1.3.5 (Katoh and

Standley 2013), and built a maximum likelihood plant diet phylogeny using FastTree

v1.0 (Price et al. 2009) with default settings in Geneious v9.1.6 (https://www.geneious.

com).

Statistical analyses

Metabarcoding data (ASV-table, taxonomic assignments, and metadata) and trait data

were imported into r v3.6.0 (R Core Team 2018) and combined using package

‘PhyloSeq’ v1.24.2 (McMurdie and Holmes 2013). Upon inspection, 40 out of 778

true plant-ASVs were found from one or more negative controls and considered to

possibly originate from contamination, and therefore removed from all samples

before downstream analyses, with associated tip labels removed from the plant diet

phylogeny.

We first created an overview of all the plant families found from the total dataset

(n

samples

= 658). For each individual we scored the occurrence of plant diet families

and calculated percentages of occurrence by target species and location for each

plant diet family. Data for all non-target species were pooled, because sample sizes for

non-target species were low (ranging from 1 to 15; Table 4.1).

In the analyses described below all ASVs identified to originate from seed plants

were included in the analyses, including ASVs for which explicit identification down

to plant family level was not possible, because we were primarily interested in genetic

variation in the diet, not plant identification. We used the Chao1 estimator (to account

for different sample sizes, i.e. different numbers of reads per individual snail sample

from metabarcoding; Chao 1984) to calculate plant diet richness (counted as the

number ASVs) for each individual snail (i.e. not pooled by species as before), and

sorted results by species and (for target species) by location. We pairwise-tested

target species diet richness differences using a Wilcoxon signed rank test with

Bonferroni correction for pairwise testing. Subsequently, using Generalized Linear

Mixed Modelling (GLMM) with package ‘glmmTMB’ v0.2.3 (Brooks et al. 2017), we

tested for the correlation between response variable ‘mean Chao1 diet richness’ (per

species) and ‘shell volume’ (used as a proxy for snail size, see above), with ‘location’ as

random effect variable. In addition to this basic model, we ran models including also

the fixed variable ‘classification’ (with options Neritimorpha, Prosobranchia, and

Pulmonata). We modelled mean Chao1 diet richness (being count data) using a log link

function, and performed model selection (based on AIC) from the full model (which

included interactions between shell volume and classification). In an alternative

GLMM, we tested for correlation between mean Chao1 diet richness with regional

‘species abundance’ (where we used regional census data from Chapter 3; Table S4.4)

instead of shell volume, which served as a proxy for ‘regional success’ of the species.

(Shell volume and abundance were themselves not significantly correlated, with

linear model results R

2

_{= 0.052, and p = 0.282.) Both models were rerun with response}

variable data not represented by the Chao1 estimate of the original ASVs, but by Chao1

estimates of richness from two agglomerated versions of the plant diet phylogeny,

in which all ASV-tip label taxa with a tip cophenetic distance smaller than 0.05 and

0.10, respectively, were pooled together. These datasets represent deeper plant diet

phylogenetic ranks, without simply resorting to ranks like, e.g. the plant family,

which in different clades might have a different phylogenetic meaning (Symons and

Beccaloni 1999).

To study the possibility of species and/or individual preferences for specific

plant taxa, we tested for phylogenetic non-randomness in the diet, applying theory

from community ecology (Webb et al. 2002). We calculated Standardized Effect Sizes

(SES) for Mean Taxon Distance (MTD) and Mean Nearest Taxon Distance (MNTD)

using package ‘picante’ v1.8 (Kembel et al. 2010), with ‘taxa’ defined by plant diet

phylogeny tip labels, and tip-swapping as the null model. In this sense, negative

SES-values indicate phylogenetic diet clustering (i.e. a possible sign of diet choice by

means of preference for specific taxonomic plant diet clusters), positive values

indicate overdispersion (again possible diet choice, but now for a most diverse diet,

e.g. to gather diverse nutrients, or spread the intake of plant-specific metabolites), and

values around zero indicate randomness (no diet choice). While MNTD describes

di-versification among recent phylogenetic relations, MPD does so for deeper relations

(Mazel et al. 2016). As before, we ran these analyses not only with the original plant

diet phylogeny (in which tip labels represent simply ASVs), but also with the two

agglomerated versions of the phylogeny, because we expected phylogenetic diet

clustering could be sensitive to the taxonomic rank of the diet items (Symons and

Beccaloni 1999, Barbour et al. 2015). We first calculated SES-values for data pooled by

the three target species (plus non-target species grouped) for the whole region, to

represent the ‘species diet’. Subsequently, we calculated values by individual

(including individuals from all species), and applied one-sample t-tests to confirm

mean values being different from zero. Because individual snails can only eat what is

surrounding them, SES-MPD and SES-MNTD analyses by individual were carried out

on datasets pruned to the location (i.e. all diet taxa from outside the individual’s

location were removed first). We used ANOVA to model the influence on SES-results

of explanatory variables ‘target species’ (including a group of non-target species

pooled), ‘location’, and ‘agglomeration’. We performed model selection based on AIC,

and used package ‘relaimpo’ v2.2-3 (Grömping 2006) to calculate the relative

importance of the explanatory variables, using metric LMG (partitioned R

2

_{values, as}

suggested by -and named after- Lindeman et al. 1980).

We studied target species and location diet differences as both Sørensen and

unweighted UniFrac distances (with the latter taking diet phylogenetic distances into

account), using non-metric multidimensional scaling (nMDS) and PERMutational

ANalysis Of VAriance (PERMANOVA; Anderson 2017), using package ‘vegan’ v2.5-5

(Oksanen et al. 2017). A Sørensen distance of 0 meant two samples with exactly the

same diet, while 1 meant no shared diet items. Similarly, a UniFrac distance of 0 meant

exactly the same diet and thus all branches from the diet phylogeny shared, while 1

meant no branches shared. To keep PERMANOVA possible, because it is computationally

demanding on large data ets, we pooled data by target species-plot combinations.

We studied the explanatory variables ‘target species’ and ‘location’, and used 4,999

permutations. To test for possible differences in variance among groups, we ran

BETADISPER on the exact same data, with the same number of permutations.

Finally, we calculated diet distances between all 658 individual snails (n

_sample

pairs

= n

samples

!/(n

samples

-n

combinations

)!/2 = 658!/(658-2)!/2 = 216,153) using, again,

Sørensen and unweighted UniFrac (based on the original, non-agglomerated diet

data) distances. From this, we calculated mean distance values for all species-species

combinations (n

_{species pairs}

= n

_species

!/(n

_species

-n

_combinations

)!/2 = 26!/(26-2)!/2 = 325)

(sensu Kartzinel et al. 2015). We applied GLMs using package ‘gamlls’ (Rigby and

Stasinopoulos 2017) to test whether diet distances (best fitted by a Weibull distribution)

could be explained by differences in snail size (‘shell volume ratio’, as in our models

explaining diet richness) and/or interspecific ‘classification distance’ (with categories

from within and across classification combinations from taxa Neritimorpha,

Prosobranchia, and Pulmonata). We performed model selection (based on AIC) and

used best models to predict correlations.

Results

We collected metabarcoding plant diet data for 658 individual snail samples,

distributed over 26 snail taxa, covering each of three main classification-groups

in the region: Neritimorpha (3 species), Prosobranchia (14 species), and Pulmonata

(9 species) (Table 4.1). We found the highest sample sizes, as anticipated, for the

three target species: A. jagori (160 samples), G. similis s.l. (105 samples), and P. concinnum

(283 samples).

Table 4.1 Sample sizes by higher-order classification, snail species, and location. Three target species, for which sample sizes are especially high because of targeted sampling, are highlighted. For sample sizes by plot, see Table S4.5.

Classification Taxon

Batangan Keruak Pangi Tandu Batu Tomanggong 2 Tomanggong K

ecil

Total

Neritimorpha Georissa kinabatanganensis Khalik, Hendriks,

Vermeulen & Schilthuizen, 2018 8 8

Georissa similis E. A. Smith, 1894 s.l. 3 26 25 6 20 25 105

Sulfurina martensi (Issel, 1874) 9 2 1 12

Prosobranchia Acmella cyrtoglyphe Vermeulen,

Liew & Schilthuizen, 2015 3 3

Alycaeus jagori Von Martens, 1859 1 56 39 21 43 160

Acmella striata Vermeulen,

Liew & Schilthuizen, 2015 7 7

Chamalycaeus sp. 1 1 2

Diplommatina calvula Vermeulen, 1993 3 3

Diplommatina gomantongensis

(E. A. Smith, 1894) 1 1

Diplommatina rubicunda (Von Martens, 1864) 1 1 4 6

Japonia kinabaluensis (E.A. Smith, 1895) 3 3

Japonia sp. 3 2 5

Leptopoma pellucidum (Grateloup, 1840) 1 1

Leptopoma sericatum (Pfeiffer, 1851) 12 1 1 1 15

Plectostoma concinnum (Fulton, 1901) 45 70 36 34 45 53 283

Plectostoma simplex (Fulton, 1901) 13 13

Ptychopatula orcula (Benson, 1850) 1 2 5 8

Pulmonata Everettia sp. 2 2

Kaliella accepta (Smith, 1895) 3 1 1 1 6

Kaliella barrakporensis (Pfeiffer, 1852) 2 2

Kaliella calculosa (Gould, 1852) 1 1

Kaliella scandens (Cox, 1872) 1 2 3

Microcystina appendiculata

(Von Moellendorff, 1893) 1 1

Macrochlamys tersa (Issel, 1874) 1 1 2 4

Videna metcalfei (Pfeiffer, 1845) 2 1 3

Videna sp. 1 1

A generalist diet

Of all rbcL read data (ASVs) that could be identified to originate from seed plants,

65.4% could be identified down to family rank (with > 70% confidence), representing

32 plant families (Figure 4.2). We found 53.1% of plant families from each of the target

species (plus a group of non-target species), and 34.3% in all locations. The top five

(in terms of omnipresence among individuals) contains large, well-known plant families,

such as Fabaceae (bean family), Asteraceae (composite family), Brassicaceae (cabbage

Figure 4.2 Plant families found from the snail diet, sorted by (A) target species and (B) location, with sample sizes (i.e. number of individual snails) per species/location given in brackets. Dot area scales with fraction of individuals per species/location in which a plant family was found. Additional red circumferences and a number highlight the top five most eaten plant families per species/location; ties for position 5 occurred in three locations. Coloured squares highlight plant families eaten by a single species or found from a single location only. Metabarcoding reads that could be assigned to the seed plants (Spermatophyta), but not to plant family rank (i.e. ‘unassigned’), were included in further analyses, but excluded from this graph.

(A) (B) 1 1 1 1 2 2 2 3 3 2 3 4 4 5 5 3 4 4 5 5 ZingiberaceaeVitaceae Urticaceae SolanaceaeRutaceae RubiaceaeRosaceae RhizophoraceaeRhamnaceae RanunculaceaePutranjivaceae PolygonaceaePoaceae PiperaceaePinaceae OrchidaceaeMoraceae MenispermaceaeMalvaceae JuglandaceaeGeraniaceae Fagaceae Fabaceae Cyperaceae Cupressaceae CaryophyllaceaeBrassicaceae Betulaceae AsteraceaeAraceae Aquifoliaceae Anacardiaceae A. jag ori (1 60) G. sim ilis s. l. (10 5) P. co ncinn um (2 83) non− targe t speci es (11 1) Di et p la nt fa m ily 1 1 1 1 1 2 2 2 1 2 3 3 3 2 4 3 2 4 4 3 5 4 3 4 5 5 5 4 5 5 5 5 5 5 Batan gan ( 78) Kerua k (12 0) Pang i (123 ) Tand u Batu (104 ) Toma nggo ng 2 (94) Toma nggo ng Ke cil (13 9) Di et p la nt fa m ily 10 20 30 40 50 60 Percentage of individuals with plant family in diet

Plant family in diet unique to A. jagori

G. similis s.l. P. concinnum

non−target species

Top five by species/location

1 2 3 4 5

location

family), and Moraceae (fig family), each present in 15-50% of individuals. In contrast,

some plant families were just found in a single target species/location, but then

always in one or two individuals only (< 1%).

Diet richness depends on snail size

Chao1 plant diet richness varied strongly among individuals and species. For the

original dataset (plant taxa represented by ASVs), species mean richness differed by

up to 15.3×, while for the 0.05 and 0.10-agglomerated datasets this was 10.3× and 9.0×,

respectively. Closer inspection by target species showed a consistent order between

the three target species in all locations, with A. jagori having the richest diet, and G.

similis s.l. the poorest (Figure 4.3A). Most comparisons at the location scale were

significant in the Wilcoxon signed rank tests (after Bonferroni correction for pairwise

testing), and highly significant at the regional scale (Figure 4.3B). The basic GLMM

model (including all species) showed that mean Chao1 diet richness and shell volume

were significantly and positively correlated, and this was again found from the best

model including fixed variable ‘classification’ from model selection (based on AIC;

Figure 4.3C, Table 4.2). Agglomerated datasets (Figure 4.4) showed a more nuanced

trend (Tables S4.6 and S4.7), where Chao1 diet richness again correlated significantly

and positively with shell volume (but less strongly), but the best models highlighted

different trends between the three classification groups, with Prosobranchia the only

group in which this trend is present.

Alternative basic GLMM models (i.e. excluding the effect of ‘classification’) showed

that mean Chao1 diet richness and abundance were significantly and negatively

correlated in the original and most agglomerated (h = 0.10) datasets (Tables S4.6 and

S4.7). Instead, best models (based on AIC) included the effect of fixed variable

‘classification’, and trends in these best models were no longer significant. Variation

in diet richness was thus not well explained by species abundance.

Signs of diet choice

SES-MPD and SES-MNTD values in the original, and SES-MPD in the 0.05-agglomerated

dataset, showed significantly negative values for G. similis s.l., a possible sign of diet

choice in this species (Table 4.3, but see Discussion; clustering seems towards the

more derived plant taxa, see Figure 4.4). In the original dataset (no agglomeration),

SES-MNTD values were significantly negative for the other target species, too,

indicating overall clustering of recent plant taxa, while no such clustering was

detected from the agglomerated datasets.

The repeated analyses of SES-MPD and SES-MNTD by individual (instead of

species) showed for many individuals (28.1% in SES-MPD; 28.9% in SES-MNTD)

significant negative SES-values, indicating diet clustering by individual (Figure 4.5A).

Furthermore, mean values (from all individuals) were negative and significantly

different from zero (for all levels of agglomeration; based on one-sample t-test), thus

showing an overall trend towards plant diet clustering and individual diet choice

(cf. Keck and Kahlert 2019) that differed markedly from the species-level results reported

in Table 4.3. ANOVA models for response variables SES-MPD and SES-MNTD showed

that each of the explanatory variables ‘target species’, ‘location’, and ‘agglomeration’

explained part of the variance, whereas the interaction between the first two was

significant, too (Figure S4.1, Tables S4.8 and S4.9). LMG relative importance values

highlighted the agglomeration to explain most of the variation, followed by the target

species and the location (for both SES-MPD and SES-MNTD; Figure 4.5B).

Signs of diet differentiation

Target species diets could not be reliably segregated in nMDS analyses (pooled by

target species and plot) based on Sørensen distance (stress = 0.35), while this seemed

partly possible based on unweighted UniFrac distance (stress = 0.20; a stress value

larger than 0.20 is considered to indicate no safe distinction between groups).

Most notably, the UniFrac diet distances between A. jagori and G. similis s.l. showed

little overlap (Figure 4.6). PERMANOVA results with these data showed that only

‘location’ explains part of the variation, and only based on UniFrac distance (p-value

= 0.019; BETADISPER p-value = 0.535, i.e. no significant differences in dispersion among

groups; Table 4.4).

Analysis of mean species diet distances (among all species combinations) showed

that ‘classification distance’ was important in explaining unweighted UniFrac

distance, but not in explaining Sørensen distance (Table S4.10). There was a significant

positive correlation between Sørensen distance and ‘shell volume ratio’ (Figure S4.2,

Table S4.11), suggesting that an increase in size difference between two species of

snail increases the distance in their plant diets. This significant positive trend was

found only for the within-Prosobranchia classification distance group when based

on unweighted UniFrac distance, and a significant negative trend was found for the

different-classification distance group, i.e. when species compared belonged to a

different classification group (Figure 4.7, Tables S4.11 and S4.12). In this case, trends

for the other classification groups were non-significant. Mean Sørensen diet distances

within target species were large (> 0.95 for all three target species), indicating that

individuals rarely share their diets. Values were smallest for A. jagori and largest

for G. similis s.l., whereas the order was exactly the opposite in UniFrac distance.

This means that two individuals of A. jagori eat more often the same diet taxa than

two individuals of G. similis s.l., but when they do so, these are diet taxa more likely

to be from phylogenetically widely different origins. P. concinnum is in between the other

two target species.

(3) (45) (26) (70) (1) * ns ns (25) (36) (56) * ** ns (6) (34) (39) ns * *** (20) (45) (21) ns *** *** (25) (53) (43) ns **** ****

Batangan Keruak Pangi Tandu Batu Tomanggong 2 Tomanggong Kecil

0 5 10 15 20 25 Ch ao 1 di et ri ch ne ss [A SV s] (A) (3) (7) (3) (8) (1) (6) (2) (160) (2) (105) (283)(13) (1) (6) (3) (2) (1) (12) (4) (3) (5) (3) (1)(15) (1) (8) ** **** **** 0 10 20 30 0 1 2 3 4 P. co ncinn um G. kin abata ngan ensis G. sim ilis s.l. K. sca nden s Japo nia sp . D. ca lvula D. rub icund a A. jag ori K. acce pta V. me tcalfe i L. se ricatum D. go manto ngen sis S. ma rtensi M. te rsa Cham alyca eus sp . Pterocycl os sp. Viden a sp. Evere ttia sp . M. ap pend iculat a J. kin abalu ensis P. sim plex L. pe llucidu m K. ca lculosa K. ba rrakp orensi s A. str iata A. cyr toglyp he Ch ao 1 di et ri ch ne ss [A SV s] Log 10_(sh ell vo lum e +1 ) [m m 3 ] (B)

Log10(shell volume + 1) [mm3]

M ea n Ch ao 1 di et ri ch ne ss [A SV s] (C) Neritimorpha Prosobranchia Pulmonata G. similis s.l. P. concinnum A. jagori Target species Classification ns Species Target species 4 8 12 16 1 2 3 4 P. co nci nn um A. ja go ri G . si m ilis s. l.

(A)

(B)

(C)

Figure 4.3 (Opposite page) (A) Chao1 plant diet richness (original data at ASV-level) by individual snail, with results sorted by target species and location, with sample sizes in brackets. Wilcoxon signed rank test results, with Bonferroni-correction for pairwise testing, are given above boxplots, with p-values defined as * < 0.05, ** < 0.01, *** < 0.001, and **** < 0.0001. (B) Same, but results from whole region combined and now for all species in our study. Additionally, grey bars indicate shell volume (as proxy for snail size), with values printed on secondary y-axis. (C) Mean Chao1 diet richness per species and location (same data as in panel B) as a function of ‘shell volume’ and snail ‘classification’. Dots represent combinations of species and location. Lines represent predictions from GLMM models with ‘location’ as random effect variable (black line: basic model with only ‘shell volume’ as fixed effect variable; coloured lines: best model with both ‘shell volume’ and ‘classification’ as fixed effect variables).

Table 4.2 Results from basic and best GLMM models explaining mean Chao1 diet richness. The basic model (top) showed that shell volume alone significantly explains mean diet richness. The best model from model selection (based on AIC), including fixed effect variable ‘classification’, showed a significant positive effect of ‘shell volume’ for each of the three classification groups, and no interaction effects. See Table S4.6 for model selection and Table S4.7 for full output from basic and best models with mean diet richness from different levels of plant diet phylogeny agglomeration. Significant terms printed in bold.

Basic model: mean Chao1 diet richness ~ log10(shell volume + 1) + (1|location)

Variable Estimate SE p

Intercept 0.968 0.257 <0.001

Log10(shell volume + 1) 0.305 0.088 <0.001 Best model: mean Chao1 diet richness ~ log10(shell volume + 1) + classification + (1|location)

Variable Estimate SE p

Intercept 0.285 0.517 0.582

Log10(shell volume + 1) 0.256 0.084 0.002 Classification (Prosobranchia) 0.755 0.485 0.120 Classification (Pulmonata) 0.987 0.491 0.044

Fi gu re 4.4 S na il p la nt d ie t p hy lo ge ny r ec on st ruc te d f ro m 1 10 b p s ec ti on o f t he c hl or op la st g en e r bc L ( se ns u H of re it er e t a l. 2 00 0) u si ng t he m ax im um li keli ho od ap pr oa ch in fa st tr ee v 1.0 ( Pr ic e e t a l. 2 00 9) w it h de fa ul t s et ti ng s. E ac h t re e t ip r ep re se nt s a p la nt d ie t t ax on r ec ov er ed f ro m t he o ve ra ll , re gi on al d at as et fr om o u r st ud y. E ac h co lo u re d ou te r ci rc le r epr es en ts th e s pe ci es d ie t o f a t ar ge t s pe ci es , w it h b la ck s qu ar es i nd ic at in g p re se nc e o f a di et ta xo n. In (A ), th e p la nt d ie t t ax a c or re sp on d t o t he A m pl ic on S eq ue nc e V ar ia nt s ( A SV s) r ec ov er ed w it h da da 2 (C al la ha n e t a l. 2 01 6) f ro m m et ab ar co di ng se que nc e da ta i n qi im e2 v 20 17 .12 (B ol ye n et a l. 20 18 ). To st ud y th e in flue nc e o f h ig he r-or de r t ax on om ic g ro up in g of t he p la nt d ie t, w e p er fo rm ed ti p ag gl om er at io n a t t ip c op he ne ti c d is ta nc es o f 0. 05 ( B) a nd 0. 10 ( C) .

(A

) n

o

ag

gl

om

er

at

io

n

(B

) h

=0

.0

5

ag

gl

om

er

at

io

n

(C

) h

=0

.1

0

ag

gl

om

er

at

io

n

A.

ja

go

ri

G

. si

m

ilis

s.

l.

P.

co

nci

nn

um

no

n−

ta

rg

et

sp

eci

es

(A)

(B)

(C

)

Ta bl e 4 .3 S pe ci es d ie ts St an da rd iz ed E ff ec t S iz es ( SE S) f or M ea n P ai rw is e D is ta nc e ( M PD ) a nd M ea n N ea re st T ax on D is ta nc e ( M N TD ) f or t he t ar ge t sp ec ie s ( pl us a p oo l o f n on -ta rg et s pe ci es ). T he t es t r es ul ts a re f ro m p er m ut at io n t es ts b as ed o n p la nt d ie t p hy lo ge ni es a nd t he o cc ur re nc es o f p la nt t ax a in s pe ci es ’ d ie ts a s s ho w n i n F ig ur e 4 .4 , a nd f or d iff er en t l ev el s o f a gg lo m er at io n ( no ne , a nd a gg lo m er at io n a t 0. 05 a nd 0. 10 t ip c op he ne ti c d is ta nc es ). A ne ga ti ve va lue su gg es ts pl an t di et phy lo ge ne ti c cl us te ri ng , z er o in di ca te s to ta l ra nd om ne ss , a nd a po si ti ve va lue ov er di sp er si on . p -v al ue s w er e ob ta in ed f ro m c om pa ri so n t o n ul l m od el r es ul ts f ro m r an do m ly s w ap pi ng t ip l ab el s ( 1,0 00 r an do m iz at io ns ). S ig ni fic an t v al ue s ( p < 0. 05 ) p ri nt ed i n b ol d. A gg lom er at io n none h = 0.05 h = 0.10 MPD MNTD MPD MNTD MPD MNTD n plant taxa SES p SES p n plant taxa SES p SES p n plant taxa SES p SES p A. jag ori 287 0.331 0.627 -3.47 5 0.001 99 3.065 1.000 4.97 5 1.000 47 3.552 1.000 5.353 1.000 G. simil is s.l. 132 -6 .885 0.001 -2.657 0.002 45 -2.216 0.015 0.605 0. 745 19 0.196 0.562 1.325 0.904 P. concinnum 329 -1.551 0.063 -1.816 0.035 94 3.594 1.000 6.124 1.000 47 4.216 1.000 6.434 1.000 non-target species 258 4.131 1.000 -3.022 0.001 95 3. 708 1.000 4.333 1.000 46 3.479 1.000 4.087 1.000

Discussion

We studied species-rich land snail communities with an abundance distribution

being well described by the UNTB (Schilthuizen 2011). Together with a seemingly

excess food availability, this led us to expect interspecific competition for food to be

very low or absent. Schilthuizen (2011), while having noted the apparent neutrality of

these snail communities, hinted at difficulties in its interpretation. Neutrality within

communities assumes functional equivalence of its members, and all its members to

belong to the same trophic guild (Rosindell et al. 2011). But these snail communities

Figure 4.5 (A) Jitter plots of Standardized Effect Sizes (SES) for Mean Pairwise Distance (MPD) and Mean Nearest Taxon Distance (MNTD) for all individual snail samples (each represented by a dot) for different levels of plant diet phylogeny agglomeration. Significant individual SES-results in black (p < 0.05), non-significant results in grey. Horizontal lines show the mean values by jitter group (with data from both significant and non-significant individual results included). Top labels show p-values for the difference of the mean from zero, based on single-sample t-test. Thus, while there is a large spread in individual SES-values, there is an overall trend of individual diet clustering. (B) Relative importance of the explanatory variables tested against SES-values using ANOVA (Table S4.9). Vertical bars indicate 95% confidence intervals from 1,000 bootstraps.

none _h=0.05 h=0.1 0 Agglomeration SE S none _h=0.05 h=0.1 0 MPD MNTD Targe t speci es Loca tion Agglo merat ion Variable LM G (% va ria tio n) Targe t speci es Loca tion Agglo merat ion 10 6 0 8 4 2 (A) (B) df=458 df=441 df=413 t=−10.3 t=−18.1 t=−21.5 p<0.001 p<0.001 p<0.001 df=458 df=441 df=413 t=−13.7 t=−20.1 t=−24.4 p<0.001 p<0.001 p<0.001 MPD MNTD −5.0 −2.5 0.0 2.5

(A)

(B)

appear not to comply with these assumptions: they contain members from several

widely different taxonomic groups, each with their own ecological characteristics,

and adult snail size covers five orders of magnitude (based on shell volume). Here we

show that (i) differences in plant diet richness among species are large, varying up to

15.3×; (ii) there is a significant positive correlation between snail size and diet richness

Figure 4.6 Non-metric Multidimensional Scaling (nMDS) of plant diet based on (A) Sørensen and (B) unweighted UniFrac distances. Each sample point represents the pool of diet taxa (ASVs) from all samples of a target species from a plot. Stress levels > 0.20 are generally considered to indicate no safe distinction between groups. Numbers refer to plots within locations (Table S4.1).

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 2 2 5 5 6 6 1 1 1 3 3 3 6 6 7 7 7 2 2 2 5 5 5 5 6 6 6 3 6 6 7 7 1 1 4 4 5 5 5 1 1 1 1 3 3 3 3 6 6 6 −0.3 −0.1 0.1 0.3 −0.2 0.0 0.2 nMDS 1 nM DS 2 Sørensen (stress: 0.35) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● _● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 2 2 5 5 6 6 6 1 1 1 3 3 3 6 6 7 7 2 ₂ 2 5 ₅ 5 5 6 6 6 3 3 6 6 6 6 7 7 7 1 1 1 4 4 5 5 5 1 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.6 −0.4 −0.2 0.0 0.2 0.4 0.6 nMDS 1 UniFrac (stress: 0.20) (A) (B) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 2 2 5 5 6 6 1 1 1 3 3 3 6 6 7 7 7 2 2 2 5 5 5 5 6 6 6 3 6 6 7 7 1 1 4 4 5 5 5 1 1 1 1 3 3 3 3 6 6 6 −0.3 −0.1 0.1 0.3 −0.2 0.0 0.2 nMDS 1 nM DS 2 Sørensen (stress: 0.35) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 2 2 5 5 6 6 6 1 1 1 3 3 3 6 6 7 7 2 ₂ 2 5 ₅ 5 5 6 6 6 3 3 6 6 6 6 7 7 7 1 1 1 4 4 5 5 5 1 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.6 −0.4 −0.2 0.0 0.2 0.4 0.6 nMDS 1 UniFrac (stress: 0.20) (A) (B) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 2 2 5 5 6 6 1 1 1 3 3 3 6 6 7 7 7 2 2 2 5 5 5 5 6 6 6 3 6 6 7 7 1 1 4 4 5 5 5 1 1 1 1 3 3 3 3 6 6 6 0.1 0.3 −0.2 0.0 0.2 nMDS 1 Sørensen (stress: 0.35) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 2 2 5 5 6 6 6 1 1 1 3 3 3 6 6 7 7 2 ₂ 2 5 ₅ 5 5 6 6 6 3 3 6 6 6 6 7 7 7 1 1 1 4 4 5 5 5 1 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.6 −0.4 −0.2 0.0 0.2 0.4 0.6 nMDS 1 UniFrac (stress: 0.20) (B) Location Batangan Keruak Pangi Tandu Batu Tomanggong 2 Tomanggong Kecil Target species A. jagori G. smilis s.l. non−target species P. concinnum

(A)

(B)

(although true only for classification group Prosobranchia when diet taxa are

agglomerated to represent higher taxonomic ranks); (iii) there are signs of significant

phylogenetic diet choice in ca. 28% of the individuals (from SES-values of MPD and

MNTD), suggesting some level of diet choice (but see below), although this trend is

weak when the diet is assessed by target species; and (iv) there is a significantly

positive correlation between difference in snail species’ size (‘shell volume ratio’)

and Sørensen diet distance, but when diet distance is calculated using unweighted

UniFrac this trend is significant only in Prosobranchia, and significantly negative

when species are from different classification groups.

We found snail plant diets to be rich, generalistic, and generally overlapping

between three target species, with all species eating roughly from the same 32 plant

families. Rich herbivore diets were explained by Freeland and Janzen (1974) possibly

to originate from inflation as a result of individuals trying to reduce the intake of

single plant-specific toxic metabolites. But this does not explain the widely varying

(up to 15.3× in plant taxa from ASVs) difference in mean diet richness among the 26

snail species in our study. We found a significant positive correlation between mean

diet richness (by species) and snail size (with ‘location’ included as random effect

Table 4.4 PERMANOVA and BETADISPER results from models explaining variation in Sørensen and unweighted UniFrac distances based on ‘target species’ and ‘location’. Sample data are based on pools of diet taxa (ASVs) from all samples of a target species from a plot. Importantly, a significant PERMANOVA result is only informative in combination with a non-significant BETADISPER result, because the variation among groups could otherwise still be explained solely by differences in group dispersions. Hence, only combinations of significant results for PERMANOVA and non-significant results for BETADISPER are in bold.

Response PERMANOVA BETADISPER

Explanatory

variable df SS pseudo-F R

2 _{Pr (>F) SS pseudo-F Pr (>F)}

Sørensen Target species 3 -0.041 -1.025 -0.085 0.995 0.018 2.295 0.086 distance Location 5 -0.132 -1.961 -0.271 1.000 0.007 0.361 0.877

Target species *

location 13 0.189 1.083 0.389 0.441

Residuals 35 0.470 0.967

Totals 56 0.486 1.000

UniFrac Target species 3 2.069 3.143 0.148 <0.001 0.032 5.170 0.003 distance Location 5 1.500 1.367 0.107 0.019 0.011 0.821 0.535

Target species *

location 13 2.727 0.956 0.195 0.664

Residuals 35 7.679 0.549

variable). Is this proof of non-neutrality? Probably not, because it is also possible that

larger snails simply eat more (because of larger mouthparts and digestive systems)

and are able to move faster (in an absolute sense), thereby covering more ground (and

plants), or that larger species unintentionally eat plants growing close to or on top of

their preferred diet plants. While such a direct size-diet richness correlation was

rejected for large mammalian herbivores (Mysterud 2000, Kartzinel et al. 2015),

we do not know of studies testing this simple correlation in snails or other invertebrates.

Alternatively, the observed correlation could be explained as a by-product of the

Jarman-Bell principle (Bell 1971, Jarman 1974, Geist 1974). This principle, based on studies of

large ungulates, but perhaps also applicable to communities of invertebrates, states

that larger herbivores need less nutrient-rich food due to their relatively large (thus

more efficient) gastrointestinal tract, and can therefore feed less specifically than

smaller species. Hence, they could feed on more different plant species, increasing

total diet richness (du Toit and Olff 2014).

Figure 4.7 Mean unweighted UniFrac plant diet distance (among all species combinations) versus shell volume ratio, for species combinations from within and between classifications, plus category ‘same species’. Dot size is scaled to the number of sample-sample comparisons available to calculate mean unweighted UniFrac distances. Model predictions from best model, which includes interaction effects of ‘shell volume ratio’ and ‘classification distance’ (Tables S4.10-S4.12). Only the model prediction trends for Prosobranchia and ‘different classification’ are significantly different from zero (non-significant trends not shown). Same species differences for the three target species are highlighted.

0.00 0.25 0.50 0.75 1.00 0 1 2 3 4

Log10(shell volume ratio)

M ea n Un iF ra c pl an t d ie t d ist an ce Classification distance

●

●

●

●

●

same species within Neritimorpha within Prosobranchia within Pulmonata different classification ● ● ● 10 100 1,000 Model predictions nsample pairs A. jagori G. similis s.l. P. concinnum

●

●

●

●

●

●

Target species

Our SES-results from phylogenetic diet clustering in individual snails (on all

three plant taxonomic levels studied, and from both MPD and MNTD) suggest

non-random diet choice in ca. 28% of individuals. The same analysis at the species

level (i.e. data pooled by species for the whole region) showed such diet clustering

only for the full dataset, and only based on the more recent plant taxa (in MNTD, not in

MPD). Thus, patterns found for individual snails do not fully transcend to the species

level. Although snails are often assumed to feed at random, several experimental

findings showed that snails can have real food preferences, both with interspecific

competition (Hatziioannou et al. 1994, Byers 2000, Riley et al. 2008) and without

(Pennings et al. 1993, Wakefield and Murray 1998, O’Rorke et al. 2016). An alternative

explanation comes from apparent diet clustering as the by-product of other,

non-dietary (selection) pressures, such as preferences for plants that allow hiding

from predators (Watanabe 1984, Alexander and Covich 1991, Nyström and Pérez 1998,

Levri 1998), or plants that offer a humid environment, helping to avoid desiccation

(Chang 1991, Lee and Silliman 2006). Hence, it is worth noting that phylogenetic diet

clustering in our study was particularly strong in G. similis s.l. Being the smallest

species in our study, it is likely more prone to desiccation than the larger species, and

thus benefits most from selecting a dark, humid place, and keep feeding there. Finally,

signs of diet differentiation were weak, with PERMANOVA of target species diets and

models of diet distances between all species showing little diet differentiation. Hence,

signs of diet differentiation, and thus possible competition, were only present at the

individual level, and not at the species level.

Snails, being semi-sessile, can eat only what is present where they live. We accounted

for this by including ‘location’ as an explanatory variable in relevant models (SES of

MPD and MNTD, and PERMANOVA of diet distances) and indeed found the location

to explain part of the variation in the response variables (roughly of the same order

as for ‘species’). An improved, but labour intensive, method would be to sample the

complete plant community from the plots, in addition to the snails, and perform plant

barcoding, which would allow the assessment of how much of the plant variation

available is truly eaten by the snails. While previous studies on snail community diet

used rather broad categories (e.g. Schamp et al. 2010), we collected data at a very high

resolution, possibly even including genetic variation within plant species, and plants

not previously recorded in the region (Azmi 1998, Boonratana 2000). As a result of

our metabarcoding marker choice of rbcL, we could not consider any non- seed plant

diet items. Therefore, our estimates of diet richness and differentiation are probably too

conserved. For example, Barker and Efford (2004) list diets for various pulmonate

families (of which several are included in our study), and these include, next to ‘live

plant material’, fungi, detritus, algae, bacteria, and (dead) animals. Ideally, our

metabarcoding for snail diet would be extended with genetic markers to pick up these

food sources, too. Furthermore, the diets we reconstructed from metabarcoding data are

snapshots in time, because snails digest their food within two or three days (Dobson

and Bailey 1982, Flari and Lazaridou-Dimitriadou 1996), for which we accounted in part by

using large sample sizes. Snails can also change their diets throughout the year

(Hatziioannou et al. 1994), and it would be worthwhile to repeat our study during

different seasons and years, because functional differences (and competition for

food) might in fact be seasonal (Hutto 1985, DuBowy 1988).

Despite these possible caveats, we believe that our approach is valuable in exploiting

the power of modern metabarcoding techniques to quickly and efficiently reconstruct

individual and species diets, study differences in diets, and search for signs of diet

choice and competition. Our approach is easily extendable to larger sample sizes,

a broader diet assessment, different trophic levels, and different regions or biomes.

Data accessibility

All sample vouchers were deposited in the BORNEENSIS collection of the Institute

for Tropical Biology and Conservation, Universiti Malaysia Sabah, Kota Kinabalu,

Malaysia (BORN), and DNA extractions at the Naturalis Biodiversity Center, Leiden,

the Netherlands, at -80°C, for future reference. Museum IDs and location details for

samples are given in Table S4.2. 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. 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. Marti Anderson commented on earlier

versions of the statistical analyses (PERMANOVA). Groningen University students

Annabel Slettenhaar and Manon Spaans collected snail trait data. 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), Malaysia.

Supplementary material

Figure S4.1 Jitter plots of Standardized Effect Sizes (SES) for (A) Mean Pairwise Distance (MPD) and (B) Mean Nearest Taxon Distance (MNTD) for all snail samples (i.e. individuals) for different levels of plant diet phylogeny agglomeration (see Figure 4.4). Dots represent single samples; significant individual SES-results in black (p < 0.05). Horizontal lines show the mean values by location.

no agglomeration h = 0.05 agglomeration h = 0.10 agglomeration

A. jag ori G . si m ilis s.l. P. co nci nn um non −ta rge t sp eci es Batan gan Kerua k Pang i Tandu B atu Toma nggong 2 Toma nggo ng Ke cil Batan gan Kerua k Pang i Tandu B atu Toma nggong 2 Toma nggo ng Ke cil Batan gan Kerua k Pang i Tandu B atu Toma nggong 2 Toma nggo ng Ke cil −5.0 −2.5 0.0 2.5 −5.0 −2.5 0.0 2.5 −5.0 −2.5 0.0 2.5 −5.0 −2.5 0.0 2.5 SE S

(A) MPD

A. jag ori G . si m ilis s.l. P. co nci nn um non −ta rge t sp eci es Batan gan Kerua k Pang i Tand u Batu Toma nggo ng 2 Toma nggong Kecil Batan gan Kerua k Pang i Tand u Batu Toma nggo ng 2 Toma nggong Kecil Batan gan Kerua k Pang i Tand u Batu Toma nggo ng 2 Toma nggong Kecil −5.0 −2.5 0.0 2.5 −5.0 −2.5 0.0 2.5 −5.0 −2.5 0.0 2.5 −5.0 −2.5 0.0 2.5 SE S

(B) MNTD

Figure S4.2 Mean Sørensen diet distance (among all species combinations) versus shell volume ratio, for species combinations from within and between classifications, plus category ‘same species’. Dot size is scaled to the number of sample-sample comparisons available to calculate mean Sørensen distances. Model predictions from best model, which excludes the effect of ‘shell volume ratio’ (Tables S4.10 and S4.11). Same species differences for the three target species are highlighted. 0.85 0.90 0.95 1.00 M ea n Sø re nse n pl an t d ie t d ist an ce 0 1 2 3 4

Log10(shell volume ratio)

Classification distance

●

●

●

●

●

same species within Neritimorpha within Prosobranchia within Pulmonata different classification ● ●

●

10 100 1,000 Model prediction nsample pairs A. jagori G. similis s.l. P. concinnum

●

●

●

●

●

●

Target species

Table S4.1 Sampling location details.

Location Plot Latitude [°] Longitude [°] Sample date Batangan 2 5.459 118.102 20 November 2017 5 5.459 118.101 20 November 2017 6 5.459 118.100 20 November 2017 Keruak 1 5.522 118.285 25 November 2017 3 5.523 118.285 25 November 2017 6 5.522 118.286 25 November 2017 7 5.523 118.286 18 November 2017 Pangi 2 5.532 118.306 24 November 2017 5 5.533 118.304 24 November 2017 6 5.538 118.306 24 November 2017 Tandu Batu 3 5.604 118.339 23 November 2017 6 5.605 118.343 23 November 2017 7 5.604 118.342 23 November 2017 Tomanggong 2 1 5.524 118.302 21 November 2017 4 5.531 118.304 24 November 2017 5 5.525 118.307 21 November 2017 Tomanggong Kecil 1 5.507 118.297 19 November 2017 3 5.509 118.299 19 November 2017 6 5.510 118.300 19 November 2017

Ta bl e S 4. 2 S am pl e d et ai ls w it h G en Ba n k B io Sa m pl e a nd S R A a cc es si on n u m be rs . Gen Bank accessions Museum ID Samp le ID Species Location Plot BioSamp le SRA BO R /M O L13403.01 KPH05000.01 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295853 SRR88417 56 BO R /M O L13403.02 KPH05000.02 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295854 SRR88417 57 BO R /M O L13403.03 KPH05000.03 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295855 SRR88417 58 BO R /M O L13403.04 KPH05000.04 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295856 SRR88417 59 BO R /M O L13403.05 KPH05000.05 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295857 SRR88417 52 BO R /M O L13405.01 KPH05002.01 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295858 SRR88417 53 BO R /M O L13405.02 KPH05002.02 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295859 SRR88417 54 BO R /M O L13405.03 KPH05002.03 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295860 SRR88417 55 BO R /M O L13405.04 KPH05002.04 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295861 SRR88417 50 BO R /M O L13405.05 KPH05002.05 Plectostoma concinnum (Fu lton, 1901) K eruak 1 SAMN11295862 SRR88417 51 BO R /M O L13408.01 KPH05005.01 G eorissa kinab atang anensis Khalik, Hendriks, V ermeu len & S chilt huizen, 2018 K eruak 1 SAMN11295863 SRR8841496 BO R /M O L13410.01 KPH05007 .01 G eorissa simil is E. A. Smit h, 1894 s.l. K eruak 1 SAMN11295864 SRR8841497 BO R /M O L13410.02 KPH05007 .02 G eorissa simil is E. A. Smit h, 1894 s.l. K eruak 1 SAMN11295865 SRR8841498 BO R /M O L13415.01 KPH05012.01 K al iella scandens (Co x, 1872) K eruak 1 SAMN11295866 SRR8841499 BO R /M O L13416 .01 KPH05013.01 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN1129586 7 SRR8841500 BO R /M O L13416 .02 KPH05013.02 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295868 SRR8841501 BO R /M O L13416 .04 KPH05013.04 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295870 SRR8841503 BO R /M O L13416 .05 KPH05013.05 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295871 SRR8841493 BO R /M O L13416 .06 KPH05013.06 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295872 SRR8841494 BO R /M O L13416 .07 KPH05013.07 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295873 SRR8841623 BO R /M O L13416 .08 KPH05013.08 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN1129587 4 SRR8841622 BO R /M O L13416 .09 KPH05013.09 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN1129587 5 SRR8841625 BO R /M O L13416 .10 KPH05013.10 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN1129587 6 SRR8841624 BO R /M O L13416 .11 KPH05013.11 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295877 SRR8841619 BO R /M O L13416 .12 KPH05013.12 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295878 SRR8841618 BO R /M O L13416 .14 KPH05013.14 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295880 SRR8841620 BO R /M O L13416 .15 KPH05013.15 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295881 SRR8841627 BO R /M O L13416 .16 KPH05013.16 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295882 SRR8841626 BO R /M O L13416 .17 KPH05013.17 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295883 SRR88417 46 BO R /M O L13416 .18 KPH05013.18 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295884 SRR88417 47 BO R /M O L13416 .19 KPH05013.19 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295885 SRR88417 44 BO R /M O L13416 .20 KPH05013.20 Plectostoma concinnum (Fu lton, 1901) K eruak 7 SAMN11295886 SRR88417 45 BO R /M O L13423.01 KPH05020.01 Jap onia sp. K eruak 7 SAMN11295887 SRR88417 42 BO R /M O L13424.02 KPH05021.02 G eorissa simil is E. A. Smit h, 1894 s.l. K eruak 7 SAMN11295889 SRR88417 40