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.

Study sites

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 Santarosae. 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.

Plant community composition

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² 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.

Bee community composition:

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).

Microbial community composition:

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 decontam tool applying a prevalence-based method (Davis et al., 2018), which identified contaminant taxa (Eisenhofer et al., 2019; Salter et al., 2014) such as: Escherichia, Salmonella, Staphylococcus, Cutibacterium, Micrococcus, and Streptococcus, 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.

Environmental and spatial variables:

To assess environmental variation within our study plots, we obtained soil and climatic variables from the Soil Survey Geographic Database (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 Soil Erodibility Index and Wind Erodibility Index, which indicates the soil's susceptibility to wind erosion. We also considered Mean Annual Precipitation and Mean Annual Air Temperature, along with Slope, Elevation (meters above sea level), and Albedo, which measures soil surface reflectivity. The selected abiotic variables have potential implications for the ecological dynamics in our multitrophic community. Soil Erodibility Index 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. Wind Erodibility 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). Slope and Elevation can influence plant growth and flower timing, which in turn can affect bee activity and species distribution (Rafferty et al., 2020). Lastly, Albedo 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.