Molecular Ecological Network Analyses: An Effective Conservation Tool for the Assessment of Biodiversity, Trophic Interactions, and Community Structure

Global biodiversity is threatened by the anthropogenic restructuring of animal communities, rewiring species interaction networks in real-time as individuals are extirpated or introduced. Conservation science and adaptive ecosystem management demand more rapid, quantitative, and non-invasive technologies for robustly capturing changing biodiversity and quantifying species interactions. Here we develop molecular ecological network analyses (MENA) as an ecosystem assessment tool to address these needs. To construct the ecological network, we used environmental DNA from feces to identify the plant and mammal diet of two carnivores: puma (Puma concolor) and bobcat (Lynx rufus); two omnivores: coyote (Canis latrans) and gray fox (Urocyon cinereoargenteus); and two herbivores: black-tailed deer (Odocoileus hemionus) and black-tailed jackrabbit (Lepus californicus) in a well-studied Californian reserve. To evaluate MENA as a comprehensive biodiversity tool, we applied our framework to identify the structure of the network, patterns of trophic interactions, key species, and to assess its utility in capturing the biodiversity of the area. The high dietary taxonomic resolution enabled the assessment of species diversity, niche breadth and overlap. The network analysis revealed a dense ecological network with a high diversity of weakly connected species and a community that is highly modular and non-nested. The significant prevalence of tri-trophic chain and exploitative competition patterns indicates (i) the removal or reintroduction of a top predator would trigger a trophic cascade within this community, directly affecting their prey and indirectly the plant communities, and (ii) the potential impact of indirect effects between two predators that consume the same prey. These results suggest that the recent resurgence of puma in the study area may impact the herbaceous and woody vegetation and the population size of other predators. This effect of fluctuating predator populations and plant communities could be predicted through MENA’s fine-scale assessment of the diet selection and the identified keystone species. Although just using a subset of species, MENA more rapidly, accurately, and effectively captured the broader biodiversity of the area in comparison to other methodologies. MENA reconstructed and unveiled the hidden complexity in trophic structure and interaction networks within the community, providing a promising toolkit for biodiversity and ecosystem management.

Global biodiversity is threatened by the anthropogenic restructuring of animal communities, rewiring species interaction networks in real-time as individuals are extirpated or introduced. Conservation science and adaptive ecosystem management demand more rapid, quantitative, and non-invasive technologies for robustly capturing changing biodiversity and quantifying species interactions. Here we develop molecular ecological network analyses (MENA) as an ecosystem assessment tool to address these needs. To construct the ecological network, we used environmental DNA from feces to identify the plant and mammal diet of two carnivores: puma (Puma concolor) and bobcat (Lynx rufus); two omnivores: coyote (Canis latrans) and gray fox (Urocyon cinereoargenteus); and two herbivores: black-tailed deer (Odocoileus hemionus) and black-tailed jackrabbit (Lepus californicus) in a well-studied Californian reserve. To evaluate MENA as a comprehensive biodiversity tool, we applied our framework to identify the structure of the network, patterns of trophic interactions, key species, and to assess its utility in capturing the biodiversity of the area. The high dietary taxonomic resolution enabled the assessment of species diversity, niche breadth and overlap. The network analysis revealed a dense ecological network with a high diversity of weakly connected species and a community that is highly modular and non-nested. The significant prevalence of tri-trophic chain and exploitative competition patterns indicates (i) the removal or reintroduction of a top predator would trigger a trophic cascade within this community, directly affecting their prey and indirectly the plant communities, and (ii) the potential impact of indirect effects between two predators that consume the same prey. These results suggest that the recent resurgence of puma in the study area may impact the herbaceous and woody vegetation and the population size of other predators. This effect of fluctuating predator populations and plant communities could be predicted through MENA's fine-scale assessment of the diet selection and the identified keystone species. Although just using a subset of species, MENA more rapidly,

INTRODUCTION
Anthropogenic impacts threaten global biodiversity, as species and their networks of interactions (ecological networks) are forced to adapt or be lost (Gonzalez et al., 2011;Dirzo et al., 2014;Munguía et al., 2016;Start et al., 2018). As a consequence of species loss and invasion, ecological networks, such as food webs and pollinator networks, are restructuring and species interactions are rewiring in real-time, fundamentally impacting whole ecosystems and their functions (Bartley et al., 2019;Daam et al., 2019). To comprehend these altered ecological networks, conservation biologists are increasingly focusing on extinction cascades caused by weakened or lost species interactions (Baiser et al., 2012;Valiente-Banuet et al., 2015;Losapio and Schöb, 2017) and are motivated to take a multitrophic perspective to address biodiversity and ecosystem functioning . However, we are in need of more rapid, non-invasive, and quantitative technologies for biodiversity assessments and monitoring for the implementation of adaptive management strategies.
Counting species (alpha diversity) rather than their interactions is often more tractable, however, it is not nearly as informative for understanding how ecosystems function (Jordan and Scheuring, 2004;McCann, 2007). Ecosystems are composed of thousands of species interactions (billions if including microbes) that directly or indirectly impact biological populations (Kéfi et al., 2016) and contribute to network robustness (Losapio and Schöb, 2017). Recent advances in network analysis enable the interpretation of large numbers of species interactions and assessment of the predictability of dynamic ecological systems (Dale and Fortin, 2010;Delmas et al., 2019). Yet, building robust ecological networks, particularly in response to perturbations and across communities that are comprised of different species, is challenging given the difficulty of observing trophic events and quantifying the strength of interactions between species. The advent of high-throughput sequencing (HTS) makes this kind of research remarkably more tractable (Roslin et al., 2019).
Innovations in HTS, and, more specifically, DNA metabarcoding, enable accurate and cost-effective biodiversity assessments at a level of taxonomic coverage and precision previously unavailable (Ji et al., 2013;Clare, 2014;Deiner et al., 2017;Pawlowski et al., 2018;Bush et al., 2019;Makiola et al., 2020). Recent studies have used metabarcoding techniques to investigate the feeding ecology of carnivores (Shehzad et al., 2012;Torre et al., 2013;Walsh, 2015;Xiong et al., 2017), herbivores (Soininen et al., 2009;Czernik et al., 2013;Kartzinel et al., 2015;Coverdale et al., 2016;Iwanowicz et al., 2016;Erickson et al., 2017;Pansu et al., 2019), and omnivores (De Barba et al., 2014;Robeson et al., 2018; for a review see Sousa et al., 2019). For example, Kartzinel et al. (2015) measured dietary niche partitioning among large herbivores and was able to address how generalist consumer species coexist on a limited range of resource types because metabarcoding allowed for unprecedented plant identification. While many of these studies targeted the diet of a single species or a trophic guild, few compared across mammal species and multiple trophic levels, a shortcoming that hinders our full knowledge of complex ecological processes . Multitrophic perspectives provide a better understanding of ecological and evolutionary processes than do typical pairwise interactions between trophic levels (Abdala-Roberts et al., 2019). Despite this, our knowledge of ecosystem function in multitrophic communities is limited to date .
Constructing multitrophic ecological networks using traditional methods is difficult and labor intensive, especially in poorly studied and highly diverse systems (Derocles et al., 2018). Non-invasively collected DNA from the environment (eDNA) [i.e., a complex mixture of genomic DNA from different organisms found in soil, water, or feces (Taberlet et al., 2018)] may help to construct complex ecological networks by rapidly detecting trophic interactions that would otherwise be impossible to observe, and to do so more quickly than traditional methods (Ruppert et al., 2019). Only recently has molecular network analysis been applied to the assessment of ecological systems. Metabarcoding studies in this context have focused on varying scales: the diet of European hake (Riccioni et al., 2018); coral reef fish (Casey et al., 2019); freshwater ecosystems ; aquatic macroinvertebrates ; and bat , forest (Evans et al., 2016), and herbivore (Kartzinel et al., 2015) communities. However, to the best of our knowledge (see Sousa et al., 2019 for a review), to date a metabarcoding study has yet to use fecal eDNA to reconstruct ecological networks and unravel trophic interactions among carnivores, omnivores, and herbivores within a terrestrial mammal community to inform biodiversity assessments.
Here we combine DNA metabarcoding and network-based approaches in a non-invasive molecular ecological network analyses (MENA) (Deng et al., 2012) and demonstrate its utility to assess biodiversity, trophic interactions, and structure within a mammalian community. We construct a food web by identifying the diets of mammals using HTS of eDNA from feces across three feeding guilds -two carnivores: puma (Puma concolor) and bobcat (Lynx rufus); two omnivores: coyote (Canis latrans) and gray fox (Urocyon cinereoargenteus); and two herbivores: black-tailed deer (Odocoileus hemionus) and black-tailed jackrabbit (Lepus californicus). To evaluate MENA as a comprehensive biodiversity assessment tool, we apply our framework to answer the following questions: (i) What is the structure of the network and who are key species? (ii) Are specific patterns of trophic interactions (i.e., trophic cascade, omnivory, and exploitative competition) occurring and prevailing within the community?; (iii) Does MENA provide an accurate portrait of biodiversity? We validate findings with long-term studies in the study area and discuss how these molecular and network metrics provide increasingly valuable information and a rapid bioassessment tool for conservation decision making.

Study Site and Feces Collection
Jasper Ridge Biological Preserve (JRBP) is a small (5 km 2 ) protected area on the urban fringe of Silicon Valley, a highly urbanized region in California, United States. It is located in the eastern foothills of the Santa Cruz Mountains and surrounded by various land use types. JRBP has a Mediterranean-type climate, consisting of warm, dry summers and cool, wet winters, and averaging temperatures of 5-22 • C across the year. The mean annual precipitation is 605 mm (Zavaleta and Kettley, 2006). JRBP has been an active site of research for over a century, protected as a preserve since 1973, and the fauna and flora of the preserve have been well-documented. Dominant vegetation types include grassland, woodland, chaparral, coastal scrub, marshes, and forests, with the most common vegetation type being oak (Perea et al., 2017).
A long-term camera-trapping effort (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) at JRBP elucidated dynamic shifts within the predator and herbivore populations, triggered when an apex predator, the puma, became resident after at least a decade of very low abundance (Leempoel et al., 2019). The (re)establishment of the puma coincides with a change in coyote behavior and abundance, which in turn, opened a niche for the smaller, omnivorous gray fox. JRBP thus provides an ideal opportunity to test the MENA tool in identifying biodiversity and patterns of trophic interactions, because the flora, fauna and dynamic shifts within this system are known.
Fecal samples (scats) from puma, coyote, bobcat, gray fox, black-tailed deer, and black-tailed jackrabbit were collected across all six habitat types in JRBP within 2-week survey windows [following (Murphy et al., 2007;Ruell et al., 2009)] in the dry season (October 26 -November 12, 2017) and wet season (April 4 -17, 2018). Before the start of the survey, all scat was removed from the trails to ensure the age of the collected samples was all less than 2 days old. This kept freshness of the sample consistent and minimized fecal DNA degradation. Sampling transects were non-random and predator scat was collected along 34 paths (trails 17 km; roads 7 km) throughout JRBP (Supplementary Figure S1). Camera trapping demonstrates that predators use these pathways at a high frequency (Kohn et al., 1999;Leempoel et al., 2019). During collection, onehalf of the scat sample was left in situ to ensure that the relevant inter-and intra-species scent cues were not disturbed. Scat samples were collected into a sterile bag, using gloves to avoid contamination, and GPS coordinates and metadata were recorded at the collection site. All samples were stored at −20 • C until DNA extraction was performed. In total, + 175 km of trails were traversed through all dominate habitats over 12 collection days in each season (finding scat on 19 of the 24 traversing days).
A combination of morphological identification (i.e., size, shape, composition) and DNA sequencing (e.g., the mammal primer confirmed all predator identifications through the most abundant predator reads per sample) enabled us to determine the identity of species that deposited the scat. Multiple jackrabbit fecal samples were genotyped to confirm the identification of the collected samples, which are visually much larger, but similar in shape, to scats of the smaller brush rabbit (Sylvilagus bachmani).

Fecal DNA Extraction and Metabarcoding
We conducted all DNA extractions at Stanford University in a sterilized laminar airflow hood to avoid contamination. Before DNA extractions, fecal samples were thawed, homogenized, and processed (∼0.2 g) utilizing the Zymo Quick-DNA Fecal/Soil Microbe Miniprep Kit (Kartzinel et al., 2015). Samples were processed in small batches (∼14) with an extraction blank to monitor for potential cross-contamination in the laboratory. The eluted DNA was quantified using a Nanodrop 2000 (Thermo Fisher Scientific Inc.) and stored at −20 • C until used for PCR.
We used DNA metabarcoding on the extracted DNA to characterize carnivore, omnivore, and herbivore diets, by quantifying large numbers of short, unique regions of DNA extracted from individual fecal samples (Pompanon et al., 2012). Two different primer pairs were used to amplify mammal and plant DNA. All primers were modified with the Illumina adaptor preceding the target primers and separated by 6-N spacers as designed by Ushio et al. (2017).

Mammals
We use the MiMammal-U metabarcoding primers for the 12S mtDNA gene targeting 210 bp amplicons because it can amplify and distinguish DNA from a diverse and wide range of mammals that are well represented in public databases (MiMammal-UF: GGGTTGGTAAATTTCGTGCCAGC and MiMammal-UR: CATAGTGGGGTATCTAATCCCAGTTTG) (Ushio et al., 2017). We compared our primer reference database to the list of known JRBP mammal species (46 species) and found that these primers identify 96% of the mammals present in JRBP to the genus and species level and 4% to the family level (Leempoel et al., 2020). The PCR comprised 20 µL: 10 µl of GoTaq R Colorless Master Mix, 1 µL of each primer (5 mM), 4 µL of DNA template, and 4 µL of water. Cycling conditions used initial denaturing at 95 • C for 10 min, followed by 35 cycles of denaturing at 95 • C for 30 s, annealing at 60 • C for 30 s and extension at 72 • C for 10 s.

Plants
The most widely applied DNA barcode for herbivore diet analysis is the P6 loop of the chloroplast trnL (UAA) gene (Taberlet et al., 2007). However, it has a lower taxonomic coverage than a recently designed second internal transcribed spacer of nuclear ribosomal DNA (ITS2) primer (Moorhouse-Gann et al., 2018). For this study, we targeted 200-387 bp amplicons of the ITS2 region for the amplification of plants in the diet using primers UniPlantF: TGTGAATTGCARRATYCMG and UniPlantR: CCCGHYTGAYYTGRGGTCDC (Moorhouse- Gann et al., 2018). We compared our primer reference database to the list of known JRBP plant species (762 species) and found that these primers identify 88% of the plants present at JRPB to the species and genus levels, 7% only to the family level, and 5% were not found in the database (GenBank). The first PCR comprised of 20 µL reactions using GoTaq R Colorless Master Mix, 0.6 µL (10 mM) of each primer, and 6 µL of DNA template and 2 µL of water. Thermocycling followed a program of initial denaturing at 95 • C for 10 min, followed by 35 cycles of 95 • C for 30 s, 55 • C for 30 s, and 72 • C for 30 s, with a 2-min final extension at 72 • C.
The integrity of DNA extracted by each primer was assessed by gel electrophoresis. Specifically, 3 µL of each DNA extract was analyzed in a 2% agarose gel and was visualized by UV illumination. If positive, the PCR products were cleaned with the Qiagen PCR Purification Kit (Valencia, CA). For the two-step PCR, appropriate Illumina barcodes were ligated to each sample as an index tag for each unique sample (Ushio et al., 2017). The index PCR was performed as a 20 µl reaction: 10 µl of Amplitaq Gold360 Mastermix reactions (with 2.5 mM MgCl2,200 lM each dNTP, 0.1 mg/mL BSA, 4% DMSO), 1 µl (of each primer), 3 µl of purified DNA amplicons and 6 µl of H 2 0. Cycling conditions used initial denaturing at 95 • C for 10 min, followed by 12 cycles of denaturing at 95 • C for 30 s, annealing at 60 • C for 30 s and extension at 72 • C for 10 s.
The indexed secondary PCR products were quantified using a Fragment Analyzer TM (Automated CE System from Advanced Analytical Technologies), normalized to equimolar concentrations and pooled together before purification using QIAquick PCR Purification Kit (Qiagen). Sequencing was performed on a MiSeq platform using the Reagent Kit Nano v3 for 2 × 300 bp (plant primers) and 2 × 150 bp (mammal primers) paired-end reads (Illumina, San Diego, CA, United States) and run at the Stanford University PAN Facility. PhiX DNA spike-in control (18-30%) was added to improve the data quality.

Sequence Filtering and Taxonomic Assignment
We used a series of conservative filtering steps developed by Leempoel et al. (2020), to retain as much of the "true" DNA diversity present in the scat samples as possible. We used software packages in Obitools (Boyer et al., 2016) and R version 3.6.3 (R Core Team, 2020) for demultiplexing and quality control. Each sequence was assigned to its sample of origin based on exact matches to both index identifier. Sequences were paired with Obitools illuminapairedend and aligned sequences with a score of <40 were discarded (Kartzinel et al., 2015). Quality scores of paired sequences were checked using FastQC, prior to adapter trimming (with a mismatch tolerance with primers of 10%) in Cutadapt (Martin, 2011), and low-quality sequences <Q30 were removed. After assignment of sequences to their corresponding samples, we used obiuniq to dereplicate reads into unique sequences, eliminated potential PCR and sequencing errors with obiclean, and kept only sequences occurring at least 10 times. See Boyer et al. (2016) for a more detailed explanation. Obiclean was applied sample by sample, with a maximum of one difference between two-variant sequences and a threshold ratio between counts of one, meaning that all less-abundant sequences were considered as variants. Only sequences with "head" (true sequences or chimera product and can have multiple variants) or "singleton" (either true sequences or chimeras but are not related to any other sequences) status in at least one sample were kept. Further, sequences whose status in the global dataset were more commonly "internal" (amplification/sequencing errors) than "head" or "singleton" were discarded (Giguet-Covex et al., 2014).
The remaining sequences were matched in Obitools against the reference databases built using EcoPCR to identify the molecular operational taxonomic units (mOTUs) (Blaxter et al., 2005;Ficetola et al., 2010). First, we downloaded all standard sequences for vertebrates and plants from EMBL ( 1 release 141) and converted the recovered file to EcoPCR format. EcoPCR was then used to simulate an in silico PCR, using the two primer pairs and maximum three mismatches, and a minimum and maximum length identical to the length of each metabarcode. The taxonomic assignment of mOTUs was performed using Ecotag, keeping only sequences with an identity ≥ 95%, and further inspecting and revising taxonomic assignments to ensure validity.
We removed sequences with relative read abundance (RRA) < 0.001% within samples to reduce the likelihood of including data that might be the result of cross-contamination (Kartzinel et al., 2015). Obvious contaminants (i.e., human DNA) and counter-marking species mOTUs were removede.g., fox and bobcat frequently urinate or defecate on each other's scats. Non-target mOTUs were removed because the MiMammal primer is known to identify bird sequences in the database (Ushio et al., 2017) and we identified 7 different bird mOTUs (Supplementary Table S4) (resulting removal of one coyote and bobcat scat) and the ITS2 primer amplified fungi and bacteria (Moorhouse-Gann et al., 2018). Canis lupus familiaris (domestic dog) was present in many fox samples yet the sequence was a 95% match to Vulpes vulpes (red fox) in the database, therefore we believe reads were misidentified and we removed this mOTU from diet analyses. Finally, all sequences with > 1% of the total reads in the negative controls were removed from our samples, a threshold that is 10 times more conservative than other recent studies (e.g., Siegenthaler et al., 2019).

Diet Composition and Biodiversity Analysis
We converted the filtered sequence read data into read abundance and occurrence data to examine the utility of both data types for food web construction. Three methods were used to quantify diet composition, as suggested by Deagle et al. (2019), including Frequency of Occurrence (FOO), which is calculated as the number of occurrences divided by the total number of samples; Percent of Occurrence (POO), rescaled to account for all food items; and the Relative Read Abundance (RRA) defined as the proportion of unique Illumina sequence reads in a sample divided by the final (i.e., after quality control) number of sequence reads in that sample. Fecal samples were sequenced separately for the plant and mammal DNA; therefore, we analyzed each independently for RRA and combined the FOO results in the case of omnivore diets.
Diet items were categorized into functional groups. Plants were identified as grasses, herbs, or woody vegetation (trees and shrubs), while small mammals were defined as species with an average adult body mass < 1.5 kg and large animals as those with a body mass > 1.5 kg. The mean percent of RRA was used to compare diet type (e.g., grass vs. trees vs. small mammals) among the focal species (Pansu et al., 2019).
To determine if there is an empirical relationship in the predator-prey mass ratio in JRBP, we compared predator body size (in grams) to prey body size (both natural log transformed) at the individual level as suggested by Nakazawa (2017) using Pearson's correlation coefficient. Predators can take down prey of a certain size or smaller and scavenging carcasses does not involving killing of prey. We have evidence from camera trap imagery that deer carcasses are being scavenged by bobcats and gray fox within JRBP, therefore these correlations were performed with and without the deer.
To capture important baseline information about the food web dynamics, we compared different diversity, niche breadth and diet overlap metrics. We determined diet diversity (raw richness, Shannon diversity), niche breadth (individual, Levins'), niche overlap (Pianka's niche overlap), and compared consumed prey species diversity (at mOTU level) following Razgour et al. (2011). The Shannon's diversity index (H) accounts for both abundance and evenness of the species present in the diet and was calculated for each species diet using the vegan package (Oksanen et al., 2019) in R version 3.6.3 (R Core Team, 2020). A species' niche breadth describes the suite of resources that it can use (Gaston et al., 1997). We determined how uniformly resources (mOTUs) are being utilized by each species using the standardized Levins' measure of niche breadth index (B A ) (Razgour et al., 2011;Lyngdoh et al., 2014): where p i is the proportion of fecal samples in which the mOTU i was found and n is the number of possible mOTUs in the diet. Pianka's adaptation of the niche overlap (O jk ) metric was used to determine dietary overlap among all pairs of target species (Pianka, 1973;Woodward and Hildrew, 2002;Brown et al., 2014;Pansu et al., 2019): where p ij is the proportion of prey species i in one carnivore species j diet, p ik is the proportion prey species i in another carnivore species k diet, n = total number of available prey species. Ojk = 0 represents no overlap, whereas a value of O jk = 1 represents complete overlap in prey species between all carnivore species.
To determine whether the number of collected fecal samples captured enough of the diet richness for this study, we used the function specaccum in the R package vegan to determine the mOTUs accumulation curve (without replacement) (Gotelli and Colwell, 2001). We used EcoSim null models (version 1 2 ) to test whether the extent of niche overlap is greater than expected by chance (Gotelli et al., 2015). We generated 1,000 simulated matrices of randomized mOTU diet composition, using the randomization algorithm 3, where the niche utilization values are reshuffled within each row of the matrix to detect nonrandom niche overlap patters (Winemiller and Pianka, 1990). The observed niche overlap (O jk ) is then compared to the simulated niche overlap values (Gotelli et al., 2015). For all of these pairwise comparisons, we deemed niche overlap to be significant when the observed value was greater than at least 99% of the simulated values.
We performed a non-metric multidimensional scaling (nMDS) ordination based on a Bray-Curtis dissimilarity matrix using the occurrences and the FOO of mOTUS to examine the patterns of niche space across different species (permutations = 999, trymax = 500, k = 2) (Kartzinel et al., 2015;Casey et al., 2019). We also performed principle component analysis (PCA) on the FOO of mammal and plant diet taxa to detect differences in the diet composition among species and to investigate whether certain taxa were more commonly found in certain diets.

Molecular Ecological Network Analyses
Molecular ecological network analyses consisted of three main steps: (i) mapping species interactions, (ii) quantifying network structure, and (iii) quantifying specific interaction patterns.
First, we built a JRBP food web as a unipartite network, which consists of one set of nodes where two species can be connected through trophic interactions (Delmas et al., 2019). Because network metrics can be impacted by mOTU identification protocols , we were conservative with the filtering process. Thus, unique sequence reference number mOTUs were grouped at species, genus, or family level, as this approach circumvents the potential taxonomic misclassification (Lupatini et al., 2014). The network was directed from predator to prey and interactions were weighted using the FOO values. To calculate network complexity and visualize the JRBP food web, we used the software Gephi ( 3 Bastian and Heymann, 2009). The Yifan Hu layout algorithm was used for the construction of the directed network and belongs to the category of force-directed algorithms (Force Atlas and Fruchterman Reingold algorithms) (Hu, 2006). Second, we quantified overall network complexity and structure calculating the following network metrics: Connectance (proportion of realized interactions among all possible ones) to quantify the level of arrangement of a network; Average degree (average number of links per species) to understand the average level of specialization of the network; Average clustering coefficient (degree to which the neighbors of a node are connected) to estimate groupings in closely connected subsets; and Path length (average distance between any pair of nodes) to determine shortest distance between species (Delmas et al., 2017;Landi et al., 2018).
To measure the relative importance of each taxon within the network we calculated: (i) Weighted in-degree centrality, which also indicates the degree to which losing a central species from the network will trigger secondary extinctions and impact the community (Elhesha et al., 2017) and (ii) the Eigenvector centrality, which measures the node's importance within the network while accounting for the importance of its neighbors, representing potential keystone species . Keystone species play a structuring role within the food web, strongly influencing the abundances of other species (Power et al., 1996). The higher the centrality values within the network, the greater the functional role within the network .
To better understand the structure of the JRBP food web, we measured modularity and nestedness. Modularity is defined as the degree to which networks are organized in discrete groups, i.e., modules. A modular network is characterized by species belonging to the same module having high connection among themselves and few connections to species belonging to other modules (Moore et al., 2017;Landi et al., 2018;Ma et al., 2020). Modularity and module partitions were computed by simulated annealing (Doulcier et al., 2016) using rnetcarto, which allowed us also to determine which species were most connected among modules and to calculate z-values for withinmodule connectivity. Nestedness is defined as the degree to which more specialist species interact only with subsets of species that interact with more generalist species (Jonsson, 2001;Bastolla et al., 2009). The role of nestedness for biodiversity depends on the type of species interactions. In mutualistic networks, nestedness enhances the number of coexisting species, but in food webs nestedness increases niche overlap among consumers and thus may prevent their coexistence (Kondoh et al., 2010). Nestedness (η) was calculated (Bastolla et al., 2009;Losapio et al., 2019) for unipartite networks as: where n ij is the number of prey species between two predator species i and j and min (n i , n j ) is the smaller of the two values. Values of nestedness η range between 0 (perfectly non-nested, full trophic complementarity) and 100 (perfectly nested).
Thirdly, we identified specific interaction patterns using network motif analysis (Milo et al., 2002). Specifically, we focused on patterns of: (i) tri-trophic chains (e.g., predator that consumes prey who feed on plants), (ii) omnivory (species at the top of a food chain that feed on both plants and animals), and (iii) exploitative competition (two predators competing for the same shared prey) (Giling et al., 2019). These three motif patterns are represented as: We used figure the RANDESU motif finder algorithm of igraph (Csardi and Nepusz, 2006) to find and count these subnetworks. Then, we compared the JRBP food web motifs with random expectation. Notice that network links were directed from predators to prey, highlighting the structure of "who eats whom" rather the "matter and energy flow." For this reason, the motif of exploitative competition involving two predators and one prey can resemble apparent competition in theory. However, apparent competition involves two prey and one predator (Wootton, 1994). Thus, we opted for a biological definition of network motifs involving indirect interactions between two predators as exploitative competition.
To assess the significance of nestedness, modularity, trophic cascade, omnivory and exploitative competition, we used randomizations of the empirical JRBP food web. Our null model randomized the food web probabilistically maintaining total number of interactions and species-specific frequencies (Jonsson, 2001;Oksanen et al., 2019). To assess statistical significance for each of these metrics we assessed whether (i) the empirical food web metric was higher or lower than the 95% CI expected by chance, and (ii) the probability, p, that a randomization was equal or higher than observed empirically. All analyses were conducted in R version 3.6.3 (R Core Team, 2020).

Bioassessment Comparison
Mammal diversity recovered from fecal DNA was compared to an expert-curated list of known mammals in the region, and both camera trap soil eDNA surveys at JRBP. The fauna and flora of JRBP are well documented, including 762 plant species (Oakmead Herbarium, 2019) and 46 mammal species (JRBP, 2019). Camera trap and soil eDNA data sets and their analysis are described in Leempoel et al. (2020) and Leempoel et al. (2019). Briefly, during the time of this study, 18 wireless camera traps were continuously recording wildlife along trails. We used images recorded from October 15, 2017 -April 20, 2018 to overlap both fecal collection sessions for a total of 2,754 capture days (153 days × 18 cameras). To determine if some scats were over-or underrepresented, we tested for a correlation between the number of capture days per species and the number of scats. The soil eDNA study was conducted in October 2017, contemporaneously with some scat collection. It relied on 12 soil samples sampled in front of 6 of the camera traps, and used the same MiMammal-U 12 s metabarcode. We compared mammal species identified in this soil eDNA study with those identified in the scat samples.

RESULTS
We collected 158 fecal samples over the two sampling periods, consisting of 15 puma, 12 coyote, 31 bobcat, 71 gray fox, 14 deer, and 15 jackrabbit samples. The majority of the samples (87%) were collected in the dry season (dry season = 137 and wet season = 12), which may be a result of seasonal usage of the JRBP by animals, seasonal differences in trail use, and/or seasonal differences in scat preservation. We combined all samples for further analysis to represent the temporal averaging of the focal species' diets and plotted sample-based speciesaccumulation curves, all of which are approaching asymptotes (Supplementary Figures S2A,D), indicating adequate sampling of guild diets for this study.

Mammal and Plant DNA Metabarcoding
We successfully sequenced either mammal DNA, plant DNA or both from 113 fecal samples post filtering (Table 1 and  Supplementary Table S1 for sample size filtering). A total of 75% of gray fox samples contained plant DNA (n = 36) and 63% contained mammal DNA (n = 30). Two of the coyote samples had both plant and mammal DNA, with one containing plant DNA only. We retained all mOTUs as representatives of the number of species that were both in the environment (biodiversity assessment) and being consumed (network assessment) at the time of collection.
From a combined 24-day sampling period, we identified a total of 212 known mOTUs with unique reference sequence numbers from 31 carnivore, 57 omnivore, and 25 herbivore fecal samples (Supplementary Table S2). These mOTUs represent 55 families and 119 genera. We combined multiple mOTUs identified to the same genus or family, making up 167 individual mOTUs (identified to each taxonomic level: 85 species/subspecies, 67 genera, 10 tribes/sub-tribes, and 5 families/subfamilies) (Supplementary Table S3). We classified 91% of the mOTUs to the genus or species/sub-species level.
The mtDNA 12S gene was sequenced to assess the mammal diet of four predator species. From the 102 samples sequenced, a total of 20,373,994 raw reads were generated with 4,459,979 reads remaining after removing host and contaminant DNA. Host DNA made up the majority of the reads for puma 63 ± 0.12%, bobcat 89 ± 0.04%, coyote 81 ± 0.09%, and gray fox 93 ± 0.02%, and therefore able to identify any unidentified scat. After taxonomic assignment, a total of three mOTUs (genera Equus, Rattus, and Felis) were removed because their abundance in the extraction-negative control was > 1% of the total reads for those mOTUs, suggesting upstream contamination. We removed samples that consisted only of host reads (i.e., samples showing no prey species other than host DNA) or samples deemed contaminated, leaving 69 samples for downstream analysis (Supplementary Table S1).
We identified a total of 16 unique mammal mOTUs from the scat samples of the predator species for downstream diet analysis ( Table 2). Predator diet diversity comprised 11 families and 13 genera. If an mOTU was identified to the genus-level and that genus was monotypic in the San Francisco Bay Area, then cf. (conferre) was indicated and the species level identification was included (i.e., Neotoma cf. fuscipes). Therefore, 75% (12/16) of the mOTUs were identified to species.
The ITS2 marker was sequenced to assess the plant diet of the two herbivore and two omnivore species (Moorhouse-Gann et al., 2018; Table 1). From the 64 samples sequenced, a total of 1,970,916 post filtering and a final 1,813,698 reads were used for downstream analysis. Plant DNA was identified as 196 unique mOTUs to the family level and below, marked by unique reference sequence numbers (Supplementary Table S5). The mOTUs represented 44 families and 106 different genera and 77 species of plants. Of the 196 unique reference sequence numbers, there were 151 individual mOTUs after combining the same family and same genus together (identified to each taxonomic level: 77 species, 61 genera, 10 tribes, and 3 families) (Supplementary Table S6). The ITS2 primers were specific in identifying plants to the species level (51% of mOTUs) and genus level (40% of mOTUs).

Diet Composition and Richness
When evaluating diet by functional groups, we found that puma consumed mainly large mammals (Figure 1 and Table 3), primarily black-tailed deer, and frequently augmented with small mammals and single occurrences of coyote, raccoon (Procyon lotor), and feral pig (Sus scrofa). The next largest predator, the coyote, consumed an even mix of large mammals, The final sample size of sequenced samples is represented in brackets ().
Frontiers in Ecology and Evolution | www.frontiersin.org   Gray fox and coyote had a combination of plant and mammal species in their diets, but RRA were not be calculated together because of differences in the primers used.
mainly black-tailed jackrabbit (Lepus cf. californicus), and small mammals, mainly broad-footed mole (Scapanus latimanus). The smaller mesopredators, the bobcat and gray fox, regularly consumed small mammals (95 and 97% of samples, respectively) including dusky-footed woodrats (Neotoma cf. fuscipes), voles (Arvicolinae), and brush rabbits (Sylvilagus bachmani macrorhinus) (Figure 1 and Tables 2, 3). All mesopredators consumed black-tailed deer, but at a very low frequency of 1-2 samples for all species. Based on PCA, it is evident that the body size of both predator and prey drive the differences in the diets of predators (Figure 2A). The two principal components explained 82.6% of the variation, with small mammals contributing 70.6% of variance to PC1 and large mammals 60.6% in the PC2. Predatorprey body mass were positively correlated when scavenged deer was included (n = 98, r = 0.55, p < 0.0001), but the correlation is stronger when deer was excluded (n = 95, r = 0.62, p < 0.0001) (Supplementary Figure S3).
Dietary prey richness was the same for the small predators [gray fox and bobcat (10 mOTUs)], 11 for the puma, but only five for the coyote, although the number of samples for the latter was low ( Table 2). The number of mammal taxa per scat ranged between 1 and 6 (individual niche breadth, mean mOTUs ± SE per species) (Supplementary Table S7). Individual niche breadth was similar among the four predator species (F 3 , 65 = 0.69, p = 0.56).
The mean plant taxa per scat sample was much larger for herbivores, with 18.63 ± 1.72 mOTUs/fecal sample, as opposed to 7.05 ± 1.06 for omnivores. The dietary dissimilarity within species was greatest for the gray fox and least for deer, while fox, deer, and jackrabbit diets were almost completely overlapping (Supplementary Figure S4). Deer are mixed feeders (i.e., grazers and browsers), and here their diet comprised of 66% woody and 33% grass and herbaceous vegetation, while jackrabbits are grazers, with 63% of their diet containing grasses and/or herbaceous plants (Figure 1 and Table 3). The plants with the highest overall FOO for the herbivores included: the trees Quercus spp. and Cedrus deodara, (FOO = 0.68, 0.40, respectively); the forbs Plantago spp. and Medicago polymorpha (0.48, 0.36); and the grasses Amelichloa spp. and Nassella spp. (0.40) (Supplementary Table S5). Note that grass was found in all deer fecal samples, however, it only represented 12% of the RRA of the diet. For the gray fox, herbs comprised just 1% of RRA, despite being frequent in the diet (FOO) ( Table 3) Table S5). Although differing in average daily intake, the gray fox diet richness was similar to that of the herbivores (Supplementary Table S7).

Dietary Diversity, Niche Breadth and Overlap
We identified diet diversity and niche overlap among carnivores, and across all feeding guilds. Between the two carnivores, the bobcat had a relatively wider niche breadth than the puma (Levins' measure: B A = 0.41 and 0.26, respectively) but approximately the same dietary diversity (Shannon diversity index: H = 2.13 and 2.03, respectively) (Supplementary Table S7). The coyote had the narrowest niche breadth and dietary diversity of all four predators (B A = 0.19, H = 1.47). The Shannon's Index shows that the herbivores and omnivores have a significantly more diverse diet than do the carnivores (F 2 , 110 = 29.65, p < 0.0001) (Figure 2C), along with a higher individual niche breadth (F 2 , 110 = 58.99, p < 0.0001) and wider niche breadth using Levins' Standardized Index (Supplementary Table S7).
Pianka's index was used to calculate the niche overlap of predator species, first considering only mammal species in the diet (Supplementary Table S8A). The diets of the coyote and the puma did not significantly overlap [O jk = 0.23, P(Obs > null) = 37%], as the diet of the coyote consisted mostly of small mammals and jackrabbits while the diet of the puma consisted mostly of deer. Gray fox and bobcat niches overlap significantly with each other [O jk = 0.77, P(Obs > null) = 99%], and to a lesser extent with the coyote [O jk = 0.26, 0.33, respectively; P(Obs > null) < 95%]. Gray fox and bobcat differ from the coyote mainly in the breadth of small mammals they consume. As expected, the bobcat does not significantly overlap with the puma [O jk = 0.39, P(Obs > null) = 65%]. However, the nMDS for Bray-Curtis dissimilarity of diet shows that when plant and mammal diet are considered, the carnivores group together, the herbivores group together, and the omnivores group in the first nMDS axis but split in the second nMDS axis (Figure 2D).
We expanded the dietary niche overlap analysis to include plant mOTUs for omnivores and herbivores (Supplementary Table S8B). The coyote shares the greatest niche overlap with the gray fox [Ojk = 0.32, P(Obs > null) = 99%] while sharing the least with the puma (O jk = 0.16). As expected, deer and jackrabbit diets have the greatest overlap [O jk = 0.50, P(Obs > null) = 99%]. Although there is a 50% niche overlap in the herbivores diets, the PCA visually represents how the woody and herbaceous plant taxa are driving a significant differences between the two species ( Figure 2B) Table S7).

Molecular Ecological Network Analyses
The diet data were used to construct an empirical food web composed of 151 plant (Supplementary Table S6) and 19 mammal mOTUs (Table 2), totaling 170 nodes and 310 weighted edges/links (predation interactions) by FOO ( Figure 3A). The network was directed with links going from the predator node to prey/plant node. The analysis of the topological properties showed that the JRBP network has an average clustering coefficient of 0.104 among the feeding guilds with an average path length of 1.608, average degree of 1.82, and a low connectance of 0.011.
The metrics to determine key species within the network identified similar species. The most central nodes were identified through the weighted in-degree centrality metric and the top plant mOTUS were: Quercus, Cedrus deodara, Plantago, Baccharis pilularis and mammals were: deer, vole (Arvicolinae), squirrel (Sciurus niger) and jackrabbit. Within the network, each node was given an eigenvector centrality score, the closer to 1 the greater the level of influence within the network. The most important plant nodes were: Quercus, Plantago, Cedrus deodara, Baccharis pilularis, Centaurea solstitialis, Brassica, Hedera helix; and the mammal nodes were: deer, vole, squirrel, and jackrabbit (for full list see Supplementary Table S9).
The JRBP community network was significantly more modular than expected by chance (Q = 0.36, 95% CI = 0.26-0.30, p < 0.05). The JRBP food web was comprised of four modules ( Figure 3B). The most significantly connected (Z-score > 1.96) species (hub connectors) within this community are deer (7.95), jackrabbit (7.08), gray fox (5.66), coyote (3.96) and bobcat (2.74) (D) Patterns of niche space across the six different focal species using non-metric multidimensional scaling (nMDS) for Bray-Curtis dissimilarity of frequency of occurrence data showing dissimilarity of carnivore (blue), omnivore (orange), and herbivore (yellow) diets. The distance between points represents the level of difference. The closer the species in the graph, the higher their similarity.

Bioassessment Comparison
We compared our fecal eDNA data to camera trap and soil eDNA surveys from the same time period of scat collection (October 2017 -April 2018), and to an expert-curated mammal list of mammals in the study area (see Leempoel et al., 2020 for camera trap and soil eDNA methodologies). The scat sample size of predators was significantly correlated with the number of camera trap days for predator occurrences -in fact the rank order abundance in each data set was the same (Spearman's r s = 1, p(2-tailed) = 0): in order, gray fox was most abundant (scat n = 71, camera n = 186), followed by bobcat (n = 31, 136), puma (n = 15, 84) and finally coyote (n = 12, 22). This study detected 73% of known (listed) mammals in the study area (excluding bats), equal to the percent of mammals detected from soil eDNA, but greater than detected of the camera trap array (59%). Through our diet analysis, we found five species undetected by the camera traps that were small mammals (Thomomys sp., Scapanus latimanus, Peromyscus spp., and Arvicolinae), two large mammals undetected by the soil eDNA survey (Ondatra zibethicus and Procyon lotor), and an additional two species previously not known from JRBP (Mus musculus and Sus scrofa) (Figure 4). These two species, the house mouse (Mus musculus) and feral pig (Sus scrofa), are not known to inhabit the preserve, suggesting that consumption of these prey may have taken place outside the preserve.
The camera trap and soil eDNA survey both captured two species that were not found in the diets of the four predators: opossum (Didelphis virginiana) and the striped skunk (Mephitis FIGURE 3 | A subset of the molecularly constructed food web of JRBP targeting the diets of black-tailed jackrabbit, black-tailed deer, gray fox, coyote, bobcat, and puma (*) represented by 170 nodes and 310 interactions. (A) Nodes are representing the predator -prey mOTUs and are color coded according to their functional group. The size of the node is proportional to the number of links connected (in-degree), and the width of the link is proportional to the frequency of occurrence the prey item was found in the diet across scat samples. The curve of the edge is always in a clockwise direction representing a weighted link from the consumer to the food item. All of the mammal mOTUs have been labeled, along with the most consumed plant mOTUs. (B) The colors in the scaled network represent the four different modular communities detected within the JRBP community and the icon represents the key species within each module. mephitis) (Figure 4; Leempoel et al., 2019Leempoel et al., , 2020. However, one misidentified skunk scat sample was sequenced, which we included in the biodiversity assessment. Ground squirrels (Otospermophilus beecheyi) are abundant in JRBP and captured by the cameras as prey items, however, were undetected in the fecal or soil eDNA. We found that neither this species nor genus were in our 12S database relying on GenBank accessions, and thus could not be detected using our methods. The soil FIGURE 4 | Venn diagram between species recorded by camera traps, and species detected with soil eDNA and fecal eDNA both using the mtDNA 12S metabarcode. Scientific names are given at the maximum rank reached with each metabarcode. Species known to be present in the study area are in black. Species absent from the study area but detected with fecal eDNA are in gray [figure modified from Leempoel et al. (2020)].
eDNA study documented two small mammal species not found in this study: the native shrew mole (Neurotrichus gibbsii) and non-native brown rat (Rattus norvegicus).
When we include both plant genera and mammal mOTUs, we found that 32% of our mOTUs have not been recorded at JRBP previously, comprising identified plant species (41/77), genera (45/61), family/tribe (13/13) and mammal taxa (17/19) in the diet. These previously unrecord species could be new species to the reserve, or our study species may be going beyond the reserve to forage. It is also possible that detection of these new taxa is caused by poor resolution or errors in the GenBank database, or unknown sequencing errors. To be more conservative with the mOTU identification, we collapsed each plant species to the genus-level, which revealed that 71% (84/119) were on the JRBP plant list.

DISCUSSION
Studies of multitrophic interactions and ecological network analysis are increasingly being used to monitor biodiversity and trophic interactions , but the fundamental assessment and quantification of these interactions is challenging. We assembled a real-world, high-resolution multitrophic food web (Figure 3), demonstrating that molecular ecological network analyses (MENA) from readily available fecal environmental DNA (eDNA) provides a powerful assessment of an ecosystem. MENA identified the plant and mammal diet composition for six large mammals across three feeding guilds, accurately capturing the biodiversity of the area and characterizing the trophic structure and interaction patterns within the community at an exceptional taxonomic resolution, thereby improving traditional food web analysis and biomonitoring. Although previous studies have validated the use of DNA metabarcoding (Roslin et al., 2019), and created bipartite networks (Kartzinel et al., 2015;Pansu et al., 2019) or even a network of networks , few studies have yet characterized multitrophic networks with high taxonomic resolution, to characterize the extent of indirect links in a terrestrial system.

Assessing Community Structure and Keystone Species
Food web structure is one mechanism to examine ecosystem organization Thompson et al., 2012;Saint-béat et al., 2015;Monteiro and Faria, 2016) and is key in forecasting the effects of trophic degradation on ecosystem function (Duffy et al., 2007;Estes et al., 2011;Naiman et al., 2012), however, is often underutilized in bioassessments. The structure of the multitrophic network of Jasper Ridge Biological Preserve (JRBP) was found to be significantly more modular and non-nested. This highlights the presence of discrete groups of predators and preys that are strongly connected among each other and less tightly with other groups. The highly modular community indicates that disturbances (e.g., defaunation and local extinctions) would spread more slowly throughout the multiple smaller communities (Tylianakis et al., 2010). The most connected species within each hub or module were, in order of connectance, deer, jackrabbit, gray fox, and coyote. The non-nestedness of this food web indicates low niche overlap among consumers and thus a more stable coexistence (Kondoh et al., 2010).
The centrality metrics we chose identified a list of candidate keystone species that play a functional role in the food web dynamics (Power et al., 1996;Delmas et al., 2019). Studies have shown that the removal of such keystone species or nodes with a high centrality value collapse the network , thus should be prioritized for monitoring to avoid cascading extinctions within an ecosystem (Dunne et al., 2002). Within Jasper Ridge, we identified species contributing most to supporting higher trophic levels, such as the oaks, plantains, and fleshy fruit bearing shrubs/trees being fundamental to the herbivore and omnivore community and the deer, jackrabbit, vole, and squirrel are primarily supporting the predator community (Supplementary Table S9). Accurately identifying the diets of these key prey within the system will contribute to the mapping of ecological feedback loops, allowing for quantitative vegetation projections (Bowman et al., 2015). For example, the two focal herbivores of this study were the most frequently consumed prey species, and had the widest niche breadth and most diverse diets. If we focused exclusively on the frequency of functional group occurrence in the diet, we might infer that the three plant groups were consumed equally among both herbivores, however, when considering the mean abundance of reads within the diet, we see a much different conclusion (Figure 1 and Table 3). We find that both species are feeding on a broad variety of plant taxa, the deer are consuming mostly woody vegetation (trees and shrubs) and the jackrabbit more herbaceous plants. Consequently, despite the high overlap of mOTUs consumed by the herbivores, the impact on the vegetation from the removal or introduction of a prey-specific predator (e.g., puma consume mostly deer, whereas gray fox consumes more rabbits) is expected to be markedly different ( Figure 2B). With our fine-scale understanding of the diet selection by these herbivores and identified keystone species, we can now better predict the impact fluctuating predator populations may have on this ecosystem.

Patterns of Trophic Interactions
Our MENA framework provides a real-time assessment of empirical species interactions and a technique to quantify and predict the patterns of trophic interactions. Trophic interactions in food webs are, however, often modeled through species traits (McGill et al., 2006;Gravel et al., 2013), for example body mass, where predator and prey body masses scale with each other in natural food webs . This predator-prey mass ratio (PPMR) is often a driver of the patterns and strength of trophic interactions within a food web, in turn contributing to food web stability (Emmerson and Raffaelli, 2004;Brose, 2010). To be confident in the identified species interactions and trophic patterns, we confirmed that the predator and prey body masses scale with each other in JRBP, adding validity to this method and allowing for stronger predictive modeling in the future . As expected according to recent results (Allen, 2014), we found the largest predator in the region, the puma, predominantly consumed deer and other large mammals (including coyote), supplemented with small mammals. The smaller mesopredators, bobcat and gray fox, had a similar diet richness to the puma, but primarily consumed small mammals, similar to diets observed nearby in the Santa Cruz Mountains (Table 2 and Figure 2A; Smith et al., 2018).
We captured the same hierarchical patterns through network motifs as the long-term camera trap study, where the natural reoccurrence of puma to JRBP triggered a top-down effect on mesopredators and herbivores (Leempoel et al., 2019). Network motifs are the basic building blocks of communities and can be identified as overrepresented patterns throughout the network. Here we focused on the three most informative patterns: tritrophic chains, exploitative competition, and omnivory (Milo et al., 2002;Stouffer and Bascompte, 2010;Paulau et al., 2015;Delmas et al., 2019). Within the JRBP food web, tri-trophic chains and exploitative competition patterns were significantly overrepresented compared to chance expectation. These motifs may emerge through the functional importance of these direct and indirect interactions within the assembly of this ecological network. In fact, the prevalence of tri-trophic chain patterns indicate that the removal or reintroduction of top predators would trigger a trophic cascade within this community, directly affecting not only their prey but also, indirectly, the plant communities on which the prey feed. Projecting the long-term consequences of these trophic cascade motifs, we predict that the recent resurgence of puma population may positively affect the woody vegetation, favoring forest regrowth and benefiting the overall ecosystem function.
Exploitative competition predicts indirect links between predators consuming the same resources and may lead to resource exclusion by some of these competitors, results that are concordant with the decrease in coyotes following an increase in pumas in the study area (Leempoel et al., 2019). These predicted network responses are also supported by the significant niche overlap of the omnivorous coyote and gray fox diets and the subsequent predator release of gray fox that occurred at JRBP due to the exclusion of coyote by puma (Leempoel et al., 2019). This suggest the small overlap between the puma and coyote diets is driven by interspecific competition and resource partitioning, while the larger overlap of the bobcat and gray fox diet points to limited competition and a sharing of resources (Gotelli and Graves, 1996). Although omnivory patterns were not significant within this network, among the six mammals we studied, the gray fox was responsible for 58% of omnivore motifs within the food web. Since omnivory is significantly reduced without fox, it is possible that a future change in fox population would shift the prevalence and importance of omnivory within the community. Taken together, our results suggest that the JRBP food web is favoring the persistence of apex predators, despite the consequences of exploitative competition, and is enhancing the regulation of herbivores and the survival of identified key woody and herbaceous vegetation. We can then generate expectations about the impacts of shifting community dynamics, which would provide a framework for management decisions.

Capturing Biodiversity With MENA
This study showcases the value of non-invasive fecal eDNA surveys as a biodiversity assessment tool, highlighted by the depth of biodiversity captured and the greater insight into the ecology of species and their interactions. We identified a rich diversity of plant and mammal mOTUs, classifying 91% to the genus or species level (Supplementary Table S2), allowing for highresolution assessments of species dietary richness, niche breadth and overlap, and diversity of the community. To the best of our knowledge, this is the first published molecular analysis of gray fox, coyote and black-tail jackrabbit plant diets. Our approach revealed a much wider diet breadth of plants than previously recorded for gray fox (107 mOTUs), jackrabbit (97) and, to a lesser extent, coyote (24). Micro histological identification techniques used to determine fox diet found only 4-20 different plant species (Wilson and Thomas, 1999;Cunningham et al., 2006), whereas we identified 61 different genera, tripling previous estimates. We identified 65 plant genera within the black-tailed jackrabbit diet at this site, while previous morphological studies only identified 14-32 plant species across numerous sites and habitats (Fagerstone et al., 1980;Wansi et al., 1992).
We compared MENA to the plant and mammal species list of Jasper Ridge, and two concurrently run surveys, soil eDNA and camera traps (Leempoel et al., 2020), to determine its effectiveness. Although this is just a subset of the larger biodiversity network in this protected area (e.g., other mammals, microorganisms, birds, pollinators), fecal eDNA detected a greater number of mammal species over a shorter period of time than the dense array of camera traps (4 cameras/km 2 ) and a greater number of species than the soil eDNA (Figure 4). Two species, the possum and skunk, were not found in the fecal eDNA sampled diets, although we did collect a skunk scat (initially misidentified as gray fox) and so included skunk in the biodiversity assessment. Puma have previously been identified consuming both skunk and possum on camera traps, and the populations of both have declined since puma returned (Leempoel et al., 2019). The low abundances of these species could explain their absence from scat samples over the scat collection period. Fecal eDNA identified a significant proportion of the known small mammal community, which were not detected by the camera trap survey. Camera traps are set up to survey either medium to large mammals or small mammal communities, but seldom both. Accurate identification from images is also challenging (De Bondi et al., 2010;Meek and Pittet, 2012;Meek et al., 2013), and consequently, small mammals are often left out of camera trap biodiversity studies, or are poorly distinguished, yet they are important indicators of ecosystem health (Rowe and Terry, 2014). The soil eDNA detected the same small mammals as fecal eDNA, however, did not detect some of the larger mammals that occur at low occupancy in JRBP, including raccoon and muskrat. Also notable, the positive correlation between scats and camera-trap images suggests that scat collection may also be used as an indicator of relative abundance for these species, similar to previous studies on Iberian lynx (Garrote et al., 2014) and jaguar abundance (Sollmann et al., 2013).
Fecal eDNA provides a perspective of a system through a species lens, rather than that of a predefined human boundaryoften used to delimit conventional surveys or assessments. For example, the diet of plants is representative throughout a herbivores range, rather than just that recorded in JRBP (29% of mOTUs genera were previously unrecorded in JRBP). MENA can also be deployed for the early detection and monitoring of non-native species that are occurring within the animal's range (e.g., two prey items never before recoded within JRBP). For example, the presence of feral pigs (Sus scrofa) in the puma diet is an early sign that this highly invasive and destructive species [$1.5 billion in economic damage to United States agriculture and environment annually (Finzel and Baldwin, 2015]) is moving into the region, as the nearest known records are 15 km + from the study area.

Current Limitations of MENA
Molecular biodiversity assessments can only be as good as the reference database on which they rely for species identification. A good example of this is ground squirrels (Otospermophilus beecheyi), an abundant species in JRBP that we know is frequently preyed upon from photographic evidence. Ground squirrels were not detected in the diet of our four predator species and the concurrent eDNA soil study because it was not in our 12S mtDNA database. This highlights the importance of collaborative initiatives, such as The Bar Code of Life data system (Ratnasingham and Hebert, 2007), to increase reference library coverage and database capacities. Only with a comprehensive database can we thoroughly assess the biodiversity of an area using a metabarcoding approach, especially in regions that have been historically understudied .
Due to the sensitivity of eDNA, it is important to be aware of confounding factors. For example, we identified bobcat and gray fox DNA in many of each other's scat samples. We have evidence from camera traps that these species countermark one another (urinate and defecate on scats from the other species), and it is highly unlikely that they are consuming each other. Therefore, it is important to know the behavioral ecology of a study species, to rule out possible cross-contaminations, while at the same time exploring the possibility of novel and poorly recognized interactions. A second example is wind-dispersed plant spores that could seasonally contaminate scat samples. We considered possible contamination from Deodar Cedar (Cedrus deodar) pollen as it is released in the autumn (September -October) (Sharma and Khanduri, 2012). However, Deodar occurred in less than half of the samples from species consuming plants. It is possible they were consuming the seed-producing cones that were readily available. A third consideration is that the prey's diet (via the gut contents) could be represented in a predator's molecular diet analysis. It is not yet feasible to detect this, however, we did not have positive PCR results for plant DNA in the majority of the predator scats, possibly due to the larger fragment size (200-387 bp) of the ITS2 amplicon region that we used (Moorhouse-Gann et al., 2018). Regardless, as a bioassessment of the region, the identified plants were either consumed by the omnivores directly or indirectly and thus still part of the functional system. Studies reviewing these and other limitations of these techniques should be considered prior to bioassessments (Clare, 2014;Alberdi et al., 2018;Pawlowski et al., 2018;Delmas et al., 2019;McGee et al., 2019;Ruppert et al., 2019;Zinger et al., 2019).

CONCLUSION
We have demonstrated that MENA is a promising tool for monitoring biodiversity, unveiling and understanding multitrophic interactions and community structure, and identifying key and vulnerable species within a terrestrial system. Only when we know what species are present and understand their function in the ecosystem, will we be equipped to protect and holistically manage these systems. The consequences of ecological network rewiring are drastic (Olivier et al., 2019), and the current rapidity of changes to ecological networks makes it difficult to detect and respond. By repeating MENA, temporal changes can be quickly identified and tracked. For example, the impact of wildlife reintroductions or changes in wildlife abundance and occupancy can be tracked and assessed through MENA and fed back into adaptive management plans to monitor for predicted impacts (Pires, 2017). Here we have explored MENA using a few mammal species scats on three trophic levels, but future assessments could drastically broaden the picture by including more vertebrates, invertebrates, and interaction types that also include pathogens and parasites. The non-invasive and quick turn-around of MENA, along with the decreasing costs of HTS, will accelerate its implementation at local and global scales (Ruppert et al., 2019).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
JM designed the study, performed the field and lab work, and wrote the first draft. JM and KL performed the bioinformatics related to DNA metabarcoding. GL and JM performed the network analysis. EH designed and supervised the research. GL, KL, and EH contributed to writing and editing the manuscript. All authors contributed to the article and approved the submitted version. FUNDING JM acknowledges support from the Philippe S. Cohen Graduate Fellowship and JRBP Mellon grant. KL acknowledges support from the Swiss National Science Foundation Early Post doc Mobility grant (P2ELP3_175075), JRBP Mellon grant. GL acknowledges support from the Swiss National Science Foundation (P2ZHP3_187938). The work was made possible by funding from the HHMI Professorship (EH).