DOI: 10.5061/dryad.k6djh9wff
title: Data from: Multitrophic assembly influences β-diversity across a tripartite system of flowering plants, bees, and bee-gut microbiomes

description: Theoretical frameworks of terrestrial community assembly often focus on
 single trophic levels (e.g., plants) without considering how complex
 interdependencies across different trophic levels influences assembly
 mechanisms. Yet, when multiple trophic levels are considered (e.g.,
 plant-pollinator, plant-microbe interactions) the focus is typically on
 network analyses at local spatial scales. As spatial variation in
 biodiversity (β-diversity) is increasingly being recognized for its
 relevance in understanding community assembly and conservation,
 considering how β-diversity at one trophic level may be influenced by
 assembly processes that alter abundance and composition of interacting
 communities at a different trophic level (multitrophic dependency) is
 critical. Here, we build on single trophic level community assembly
 frameworks to explore the assembly processes affecting β-diversity in
 multitrophic communities comprising flowering plants, their bee
 pollinators, and the corresponding bee-gut microbiota to better understand
 the importance of multitrophic dependency in community assembly. Using
 distance-based redundancy analysis and variation partitioning, we
 investigated community assembly processes across three interconnected
 trophic levels in two ecological regions in southern California: the Santa
 Monica Mountains and three islands of the Channel Island Archipelago. We
 found that the deterministic effects of multitrophic dependency are
 stronger on directly connected trophic levels than on indirectly connected
 trophic levels (i.e., flowers explain bee communities and bees explain
 bee-gut bacteria communities, but flowers weakly explain variation in
 bee-gut bacteria communities). We also found notable regional variation,
 where multitrophic dependency was weaker on the Channel Islands as
 ecological drift was more pronounced. Our results suggest that integrating
 the influence of multitrophic dependency on community assembly is
 important for elucidating drivers of β-diversity and that multitrophic
 dependency can be determined by the regional context in which β-diversity
 is measured. Taken together, our results highlight the importance of
 considering multiscale perspectives – both multitrophic and multiregional
 – in community assembly to fully elucidate assembly processes.
description: <em>Study sites</em> We surveyed
 two ecological regions: Three islands within the Channel Island
 archipelago in the coast of California, namely Anacapa, Santa Cruz, and
 Santa Rosa islands, and five sites in canyon localities near Zuma Canyon
 and Stunt Ranch, belonging to the Santa Monica Mountains (hereafter SAMO)
 in California (Figure 2). The Channel Islands archipelago, which includes
 the aforementioned islands, originated approximately 20 million years ago
 through a complex interplay of plate tectonics and volcanic activity
 (Atwater, 1998; National Park Service, 2016). Over time, the islands have
 undergone dynamic changes influenced by geological processes. During the
 Pleistocene ice ages, sea levels were lower due to the accumulation of ice
 in polar regions. This lowered sea level exposed the seafloor between the
 Channel Islands and allowed them to be connected into a supper-island
 referred to as <em>Santarosae</em>. Although never connected
 to the mainland during the Quaternary, this land bridge facilitated the
 movement of terrestrial species between the islands. Conversely, during
 interglacial periods when the ice melted and sea levels rose, the islands
 became geographically distinct and independent entities, as they are today
 (Rick et al., 2005, 2012). Alongside the formation of the Channel Islands,
 the Santa Monica Mountains emerged in the mainland area as a result of
 plate tectonics (National Park Service, 2015). Consequently, both systems
 originated around the same time through similar geological processes, but
 they configure different ecological regions.
 <em>Plant community composition</em>
 To explore β-diversity across our study sites, we undertook a
 survey of floral communities within coastal shrublands. We collected data
 across two separate visits during April 2021 to May 2021, to align with
 the peak plant blooming season. Across a total of 271 2 m² plots (141
 plots in the Channel Islands and 130 plots in SAMO) we documented the
 identity of the flowering plant species present and the number of
 individual plants of each species per plot. Plant taxonomy follows the
 Jepson Manual (Baldwin et al., 2012). For our field study, we encountered
 varying permit restrictions at different locations. We were limited to
 sampling only in Cherry Canyon on Santa Rosa Island and in the East islet
 on Anacapa Island. At other study sites, our permits allowed sampling
 throughout all areas. To select the plots, we established a virtual grid
 over accessible valley shrubland areas at each site, with each grid cell
 measuring 2m<sup>2</sup> and representing a potential plot
 center. We then used an automated process to generate random numbers for
 selecting specific plots for sampling, requirement that no two plots be
 less than 10 meters apart. Once in the field, if a selected plot lacked
 flowering plants, we relocated the plot to the nearest site with flowering
 vegetation, while still respecting the minimum 10-meter spacing
 rule. <em>Bee community
 composition:</em> To comprehensively characterize
 bee visitors, two trained observers conducted 10-minute surveys at each
 plot using entomological nets when temperatures would align with bee
 activity (17-25° C). We collected all observed bee specimens, and we
 promptly placed them in a dry shipper containing liquid nitrogen to ensure
 the preservation of gut microbial communities until samples were stored at
 -80 ºC in the Entomology Department at the University of California,
 Riverside. Once in UCR, we identified the bee specimens to species or
 morphospecies at the McFrederick Lab following the keys of Michener et al.
 (1994), Michener (2007), and Ascher & Pickering (2010).
 <em>Microbial community composition:</em>
 To analyze the composition of gut-associated bacteria within the
 collected bee specimens, we carefully dissected the bee-gut within a
 sterile environment. We then applied a bead-beating procedure for 3
 minutes at 30 hertz to break the gut tissue. Following this, we extracted
 the DNA from the gut samples using the Qiagen DNeasy blood and tissue
 extraction kit and included blank samples as negative controls.
 Subsequently, we prepared amplicon libraries using 16S rRNA gene primers
 targeting the V5-V6 region, specifically primers 799F and 1115R. These
 selected primers allow paired-end sequencing with inline barcodes while
 limiting the amplification of plant chloroplasts and mitochondria, a
 characteristic corroborated by prior studies (Hanshew et al., 2013; Kembel
 et al., 2014). To standardize DNA content across libraries, we normalized
 DNA content using SequalPrep normalization plates (Invitrogen) following
 the manufacturer's recommendations. We established a library pool by
 combining 5 μl from each normalized barcoded library. To eliminate
 primer-dimers and excess master mix components, we conducted library
 purification using AMPure XP beads (Beckman Coulter). We used a
 Bioanalyzer (Agilent) to assess quality of our pooled libraries before
 sequencing on an Illumina MiSeq instrument, employing the V3 2 × 300
 reagent kit, a process carried out at the Genomics Core of UC
 Riverside. To assess quality control, eliminate
 chimers, conduct read denoising, and generate amplicon sequence variants
 (ASVs), we first used the DADA2 pipeline with default parameters (Callahan
 et al., 2016). For the assignment of ASV taxonomy, we utilized the QIIME2
 sklearn classifier, which we trained the specific region of the 16S rRNA
 gene that we amplified (799F-1115R). Our reference point was the SILVA 138
 SSURef NR99 full-length sequences and the associated taxonomy databases
 (Bokulich et al., 2018; Quast et al., 2013). Furthermore, we refined our
 dataset by eliminating ASVs classified as chloroplasts, mitochondria, or
 unidentified Eukaryotes. We assessed contamination by using the
 <em>decontam</em> tool applying a prevalence-based method
 (Davis et al., 2018), which identified contaminant taxa (Eisenhofer et
 al., 2019; Salter et al., 2014) such as: <em>Escherichia,
 Salmonella, Staphylococcus, Cutibacterium, Micrococcus</em>, and
 <em>Streptococcus, </em>which we removed from our dataset.
 After filtering and contaminant removal, we rarefied our sequence
 libraries by random subsampling to 800 reads, where our ASV accumulation
 curves flatten. This allowed us to retain most samples while still
 capturing the majority of ASV diversity. We exported this rarefied feature
 table for further microbial analyses, and we conducted all bacterial
 analyses at the ASV level except when otherwise indicated.
 <em>Environmental and spatial
 variables:</em> To assess environmental variation
 within our study plots, we obtained soil and climatic variables from
 the <em>Soil Survey Geographic Database</em> (SSURGO, 2015)
 using a raster map with 10-meters cell size. This database is the most
 detailed database for abiotic variables available for our study sites.
 Using ArcGIS PRO, we aligned these variables with our plot coordinates
 (Figure A1 – A8, Supporting Information). The variables include
 <em>Soil Erodibility Index</em> and <em>Wind Erodibility
 Index</em>, which indicates the soil's susceptibility to wind
 erosion. We also considered <em>Mean Annual Precipitation</em>
 and <em>Mean Annual Air Temperature</em>, along with
 <em>Slope</em>, <em>Elevation</em> (meters above
 sea level), and <em>Albedo</em>, which measures soil surface
 reflectivity. The selected abiotic variables have potential implications
 for the ecological dynamics in our multitrophic community. <em>Soil
 Erodibility Index</em> classifies soils according to their texture,
 mineralogy, organic matter, structure, permeability, and depth to a
 limiting layer. These characteristics can influence vegetation structure
 and plant species diversity (Rodrigues et al., 2018), which consequently,
 can affect the activity patterns of bees. Furthermore, soil features can
 affect bee nesting resources (Grundel et al., 2010), influencing bee
 species distribution. <em>Wind Erodibility</em> can alter
 floral landscapes, bee flying activity, and bee nesting sites, indirectly
 impacting bee populations and their associated gut microbiota (Hennessy et
 al., 2020, 2021).
 <em>Slope</em> and <em>Elevation</em> can
 influence plant growth and flower timing, which in turn can affect bee
 activity and species distribution (Rafferty et al., 2020).
 Lastly, <em>Albedo</em><strong> </strong>impacts
 the local temperature regime and moisture retention of the soil (Guan et
 al., 2009), further influencing the microhabitat conditions for both
 floral communities and their pollinators. To account for potential
 collinearity of our environmental variables, we conducted a Principal
 Component Analysis (PCA) on scaled environmental variables (Supporting
 figures A9-A10) separately for each area (SAMO vs. Islands) as we expected
 that each area would have different environmental variables that
 characterize the main environmental gradients. We retained the first 4 PC
 axes which combined described 97% of the variation in the environmental
 data for SAMO and 85% for the Islands. In the mainland localities (SAMO),
 PC1 represents a gradient of Wind Erodibility and Slope and describes 58%
 of the environmental variation; PC2 represents Mean Annual Precipitation
 and describes 21% of the environmental variation; PC3 represents a
 gradient of Mean Annual Air Temperature and Elevation and describes 12% of
 the environmental variation; and PC4 represents Albedo and describes 6% of
 the environmental variation. In the islands, PC1 represents a gradient of
 Mean Annual Precipitation and Mean Annual Air Temperature and describes
 34% of the environmental variation; PC2 represents a gradient of
 Elevation, Slope and Soil Erodibility; and describes 26% of the
 environmental variation; PC3 represents Wind Erodibility and describes 13%
 of the environmental variation; and PC4 represents Albedo and describes
 12% of the environmental variation.
description: # # README This folder contains: 1\) Script1.R: Custom R script that
 allows full reproducibility of the null models for Standardize Size Effect
 calculations to test if observed β-diversity differed from stochastic
 assembly (ecological drift). 2\) Script2.R: R script that allows full
 reproducibility of variation partitioning of β-diversity using dbRDA and
 the Blanchet et al. (2008) approach. We organized this script in 7 parts
 which need to be run in order. 3\) Figure3.R: R script to create the NMDS
 representations portrayed in Figure 3. 4\) Figure4.R: R script to create
 the observed and SES deviations of β-diversity portrayed in Figure 4. 5\)
 Figure5.R: R script to create the stacked bar plots portrayed in Figure 5.
 6\) env_all_var.csv: File with abiotic environmental variables to use with
 Script1.R file. 7\) plants-site-all.csv: the plant species matrix. Species
 names are in the file "flower-names.xlsx". 8\)
 bees-site-all.csv: the bee species matrix. Species names are in the file
 "bee-names.xlsx". 9\) flower-names.xlsx: Flowering plant species
 names, which are coded with the first 3 letters of the genus and first 3
 letters of the epithet. 10\) bee-names.xlsx: Bee species names, which are
 coded with the first 3 letters of the genus and first 3 letters of the
 epithet. 11\) ASV-site.cvs: Bee-gut bacteria matrix. 12\)
 Obs-bacteria-beta.csv: Observed beta diversity in bacterial communities.
 To use with Figure 4 R script; data obtained from running Script1.R file.
 13\) Obs-bees-beta.csv: Observed beta diversity in bee communities. To use
 with Figure 4 R script; data obtained from running Script1.R file. 14\)
 Obs-plants-beta.csv: Observed beta diversity in plant communities. To use
 with Figure 4 R script; data obtained from running Script1.R file. 15\)
 bees-null-dataframe.csv: Null deviation (SES) calculated for bee
 communities. To use with Figure 4 R script; data obtained from running
 Script1.R file. 16\) plants-null-dataframe.csv: Null deviation (SES)
 calculated for plant communities. To use with Figure 4 R script; data
 obtained from running Script1.R file. 17\) bacteria-null-dataframe.csv:
 Null deviation (SES) calculated for bacteria communities. To use with
 Figure 4 R script; data obtained from running Script1.R file. 18\)
 bees-site-islands.csv: the bee species matrix for island localities.
 Species names are in the file "bee-names.xlsx". 19\)
 bees-site-mainland.csv: the bee species matrix for mainland localities.
 Species names are in the file "bee-names.xlsx". 20\)
 plants-site-islands.csv: the plant species matrix for island localities.
 Species names are in the file "flower-names.xlsx". 21\)
 plants-site-mainland.csv: the plant species matrix for mainland
 localities. Species names are in the file "flower-names.xlsx".
 22\) RDA-bacteria-abiotic.csv: To use with Figure 5 R script, data
 obtained from running Script2.R file. 23\)
 RDA-bacteria-multitrophic-bees-flowers.csv: To use with Figure 5 R script,
 data obtained from running Script2.R file. 24\)
 RDA-bacteria-multitrophic-bees.csv: To use with Figure 5 R script, data
 obtained from running Script2.R file. 25\) RDA-bees-abiotico.csv: To use
 with Figure 5 R script, data obtained from running Script2.R file. 26\)
 RDA-bees-multitrophic.csv: To use with Figure 5 R script, data obtained
 from running Script2.R file. 27\) RDA-plants-abiotic.csv: To use with
 Figure 5 R script, data obtained from running Script2.R file. 28\)
 RDA-plants-multitrophic.csv: To use with Figure 5 R script, data obtained
 from running Script2.R file.