Introduction

Understanding the ecological consequences of biodiversity loss is an increasingly important task in ecology, given the ongoing biodiversity crisis (Isbell et al. 2022). Representing the interdependencies among organisms, ecological networks reflect whether and how species interact with each other across trophic levels, playing an indispensable role in assessing ecosystem stability and integrity (De Ruiter et al. 1995, Harvey et al. 2017). Changes in network structure usually coincide with variations in diversity and community of trophic levels, which could be in turn affected by the changes in producers via trophic cascades (Barnes et al. 2018, Gonzalez et al. 2020). However, we still lack a generalizable framework for how these networks and especially their phylogenetic interdependencies respond to basal changes of biodiversity loss in ecosystems, such as changes in the tree diversity of forests (Tylianakis et al. 2008, Grossman et al. 2018). To better understand species codependence and its role for biodiversity conservation, we must further study the complex dynamics of networks from multiple angles (Tittensor et al. 2014, Brondizio et al. 2019). Interactions can be viewed from a top-down or bottom-up perspective. Higher trophic levels can impact lower trophic levels through antagonistic species interactions, while lower trophic levels can provide resources to support higher trophic levels, Previous studies have shown there might be asymmetric effects between top-down and bottom-up control by different mediators (Vidal and Murphy 2018), with both factors shaping multitrophic communities together(Hunter et al. 1992). Moreover, the diversity and community of higher trophic levels could be driven by environmental microclimate, which could be determined by plant structuring, such as canopy cover (Fornoff et al. 2021, Perlík et al. 2023). Understanding how these forces respond to changing biodiversity is therefore imperative to predict how ongoing environmental changes will impact the functioning and stability of ecosystems (Hines et al. 2019).

A crucial network that unites bottom-up and top-down processes in many ecosystems and is prone to strong alterations due to environmental change is the interaction network between parasitoids and their hosts (Tylianakis et al. 2006, Jeffs and Lewis 2013). Insect parasitoids attack and feed on and eventually kill their insect hosts (Godfray and Godfray 1994). Parasitoids are thought to be particularly sensitive to environmental changes, because species in higher trophic levels usually have smaller population sizes and the potential specialization of their hosts may cascade up to impact the parasitoids’ presence or absence. Therefore, insect host-parasitoid systems are ideal for studying the relationships between habitat-level changes and species interactions (Jeffs and Lewis 2013). Previous studies on host-parasitoid interactions mainly focused on the influence of abiotic factors, such as elevation and habitat structure (e.g. Valladares et al. 2012, Maunsell et al. 2015, Grass et al. 2018), or sometimes interaction structure (e.g. Cagnolo et al. 2011). However, the role of multiple components of plant diversity in modifying host-parasitoid interaction networks, as key biotic determinants of overall ecosystem structure, remains relatively little explored. Recent studies mainly focused on basic diversity associations (Ebeling et al. 2012, Schuldt et al. 2019, Guo et al. 2021). These studies have demonstrated both direct and indirect effects of plant diversity on both hosts and parasitoids, possibly via increased niche space and resource availability (Guo et al. 2021). Nevertheless, how these patterns propagate to their interaction networks is still unclear. Moreover, the effects of changing plant diversity are not always obvious when looking at plant species richness, and other diversity components such as plant phylogenetic diversity have been shown to better predict diversity-dependent bottom-up effects on host-parasitoid networks (e.g. Staab et al. 2016). It is especially important to take phylogenetic dependences into account when it comes to the phylogenetic structure within and across trophic levels (Webb et al. 2002, Emerson and Gillespie 2008). This makes it vital to account for multiple dimensions of biodiversity and relevant trophic interactions (e.g. taxonomic, phylogenetic, functional; Peralta et al. 2015, Volf et al. 2017, Wang et al. 2020). These components might jointly affect host-parasitoid networks in a system with high species diversity. Forests have garnered special attention lately, because they represent complex and large ecosystems susceptible to global change (add citation here). Understanding how such relationships modulate the effects of tree diversity loss on the structure and interaction strength of host-parasitoid networks clearly requires further study (Staab et al. 2016, Fornoff et al. 2019).

Here, we leverage standardized trap nests for solitary cavity-nesting bees, wasps, and their parasitoids in a large-scale subtropical forest biodiversity experiment to test how multiple dimensions of tree diversity and community composition influence host-parasitoid network structure. A multi-faceted approach is particularly important when considering that associations between trophic levels might be non-random and phylogenetically structured (Volf et al. 2018, Wang et al. 2020). We aimed to quantify the strength and stability of associations between hosts and parasitoids to both discern the primary components of the diversity and composition of tree communities that affect higher trophic level interactions We expected that (a) multiple tree community metrics, such as taxonomic, phylogenetic and functional diversity and phylogenetic and species community composition can structure host and parasitoid community compositions, especially via phylogenetic processes, as species interactions often show phylogenetic conservatism (e.g. Pellissier et al. 2013, Peralta et al. 2015). Further, we hypothesized that (b) host-parasitoid networks will be more complex and stable with increasing tree species richness due to potential links from higher richness of hosts and parasitoids promoted by more tree species, and (c) both, community and interaction network changes can also be related to abiotic factors, such as canopy cover, which might play a role in structuring hymenopteran communities (Haddad et al. 2011, Fornoff et al. 2021). By better understanding these dynamics, we can begin building a generalized framework for understanding host-parasitoid interactions in forest ecosystems.

Results

Overall, 34,398 brood cells were collected from 13,267 tubes across five years of sampling (2015, 2016, 2018, 2019, and 2020). Six families of hosts and seventeen families of parasitoids were identified. Among them, we found 56 host species (12 bees and 44 wasps) and 50 parasitoid species (38 Hymenoptera and 12 Diptera). The full species list and their abundances are given in Table S1.

Community composition of hosts and parasitoids

Host species composition was significantly related to the species and phylogenetic composition of the tree and parasitoid communities (NMDS scores), as well as to canopy cover, tree MPD, elevation, and eastness (Fig. 1a, Table 1, Table S2). Parasitoid species composition was significantly associated with host phylogenetic diversity, tree functional diversity, tree MPD, eastness, and elevation, and was significantly related to tree species composition, host species composition, and canopy cover (Fig. 1b, Table 1, Table S3). Host phylogenetic composition was affected by tree species composition, tree MPD, tree functional diversity, canopy cover, eastness, elevation and was especially obviously affected by parasitoid species and phylogenetic composition (Fig. 1c, Table 1, Table S4). For parasitoid phylogenetic composition, significant relationships were found with tree species and phylogenetic composition, host species composition, tree functional diversity, canopy cover, elevation, and eastness (Fig. 1d, Table 1, Table S5).

Fig. 1.

Table 1.

An effect of host composition on the composition of the parasitoid communities was further indicated by a significant parafit test (p = 0.032), suggesting nonrandom associations in the phylogenetic structure of parasitoid and host communities (Fig. 2, Fig. S3).

Fig. 2.

Host-parasitoid network associations

The linear regression model results showed that host vulnerability, linkage density, and robustness of parasitoids were significantly negatively related to tree species richness, while remaining unaffected by other environmental covariables (Table 2, Fig. 3; except for elevation, which was marginally significantly related to robustness). Interaction evenness was significantly negatively associated with canopy cover, and interaction evenness was also negatively related to eastness (Fig. 4c; Table 2).. Generality was only marginally associated with canopy cover, and was not related to tree species richness or the other environmental factors. In the alternative models (tree species richness replaced by tree MPD), vulnerability and linkage density were significantly positively related to tree MPD (Fig. 4a, 4b; Table S6), while parasitoid robustness was negatively related to tree MPD (Fig. S4, Table S6). The results of other network metrics (generality and interaction evenness) were consistent with those of the primary models The results of the null model analysis suggested that our metrics calculated by the observed network were significantly different from a random distribution (72, 71, and 77 out of 85 values for generality, vulnerability, and linkage density, respectively; all values for robustness, and interaction evenness), strongly demonstrating that interactions between species were not driven by random processes.

Fig. 3.

Fig. 4.

Table 2.

Discussion

Our study demonstrates that tree species richness and phylogenetic diversity play key roles in modulating communities of hosts and parasitoids and their interactions with each other. These interactions are further structured by the phylogenetic associations between hosts and parasitoids. Moreover, canopy cover partly determined host-parasitoid networks, supporting a recent finding that the structure of host-parasitoid networks is also mediated by changes in microclimate to some extent (Fornoff et al. 2021). These patterns were highly associated with multiple tree diversity metrics (tree taxonomic, phylogenetic and functional diversity), and compositional changes which are key to understanding how host-parasitoid interactions may be impacted by biodiversity loss through trait- and phylogeny-based processes.

Community composition associations

Host species composition was influenced by several factors, including the taxonomic and phylogenetic composition of trees and parasitoids, tree diversity (species richness and MPD), and other biotic factors. The effects of tree diversity and composition on host species composition agree with previous studies where solitary bee and wasp species composition were related to plant community structure (e.g. Loyola and Martins 2008). It seems likely that these results are based on bee linkages to pollen resources and predatory wasp linkages to the diverse of food sources, which may themselves be closely linked to tree species richness (Reitalu et al. 2019, Staab and Schuldt 2020).

We also found that tree MPD, FD, and species composition affect parasitoid species composition, as many studies have found significant relationships between plant and parasitoid diversity in forests, including tree phylogenetic diversity (Staab et al. 2016), functional diversity (Rodriguez et al. 2019), and structural diversity (Schuldt et al. 2019). Similar to predators, parasitoids might also be more active or efficient with increasing tree community dissimilarity due to higher prey resources or lower intraguild parasitism caused by more diverse habitats (Finke and Denno 2002). On the other hand, our results also show that host species composition and parasitoid species composition relate to each other and their phylogenetic compositions, which are structured by tree communities to some extent. This pattern could propagate to the adjacent two higher trophic interactions through both top-down and bottom-up control.

Both host and parasitoid phylogenetic composition were related to tree species composition. This pattern has important implications for cascading effects among trophic levels, in that producer communities could structure a higher trophic level community via an intermediate trophic level. As previous studies usually found weak effects of plants on higher trophic levels (e.g. Cappelli et al. 2022). However, while both host and parasitoid phylogenetic composition was related to tree MPD, only parasitoids responded to tree phylogenetic composition. This may be because there are many caterpillar-hunting wasps in our host communities, and the community composition of caterpillars were usually correlated with tree phylogenetic communities(Wang et al. 2019). Therefore, the prey highly associated with tree phylogenetic composition (e.g. caterpillars) might indirectly determine predatory wasp (host) phylogenetic composition, as recently found for plants-caterpillars-spiders (Chen et al. 2023). This could be further tested by collecting the food directly used by the wasps (caterpillars). For parasitoids, tree phylogenetic composition might drive the process of community assembly through trophic cascades (e.g. from plants to parasitoids via herbivores and host wasps) (Webb et al. 2002, Cavender-Bares et al. 2009). Additionally, parasitoid phylogenetic composition can be influenced by tree structural diversity (e.g. host availability in plots with higher heterogeneity; Schuldt et al. 2019), which can be determined by conserved traits across phylogenies (Webb et al. 2002). The phylogenetic associations between hosts and parasitoids exhibited a nonrandom structure (significant parafit correlation) between the phylogenetic trees of the host and their parasitoids (see also Peralta et al. 2015).

Interestingly, we found that only host phylogenetic composition was affected by parasitoid phylogenetic composition (and not vice-versa). This asymmetry suggests that top-down control (parasitoids to hosts) was stronger than bottom-up control (hosts to parasitoids), supporting prior hypotheses, demonstrating the strong control of parasitoids (Vidal and Murphy 2018). This pattern might reflect disproportionate interactions between hosts and parasitoids, with one trophic level more strongly influencing the other in some instances, but further study is needed.

Moreover, the species and phylogenetic composition of hosts and parasitoids was also related to abiotic factors, especially to canopy cover, which has been considered especially important (Sobek et al. 2009, Fornoff et al. 2021). In future studies, it will be useful to incorporate other, more direct metrics of microclimate, such as local temperature and humidity, to determine the proximal drivers of these microclimatic effects (Ma et al. 2010, Fornoff et al. 2021).

Community-level host-parasitoid networks

Tree community species richness did not significantly influence the diversity of hosts targeted by parasitoids (generality), but caused a significant increase in the diversity of parasitoids per host species (vulnerability). This is likely because niche differentiation often influences network specialization via potential higher resource diversity in plots with higher tree diversity (Lopez-Carretero et al. 2014). For the significant relationship between vulnerability and tree species richness, a potential explanation is that host-parasitoid interactions could be driven through bottom-up effects. Moreover, higher trophic levels will particularly benefit from tree diversity, inducing parasitoid species increases more than host diversity with increasing tree species richness (Guo et al. 2021). These will further increase vulnerability at community level. According to the enemies hypothesis (Root 1973), which posits a positive effects of plant richness on natural enemies, the higher trophic levels in our study (e.g. predators and parasitoids) would benefit from tree diversity and regulate herbivores thereby (Staab and Schuldt 2020). Indeed, previous studies at the same site found that bee parasitoid richness and abundance were positively related to tree species richness, but not their bee hosts (Fornoff et al. 2021, Guo et al. 2021). Because our dataset considered all hosts and reflects an overall pattern of host-parasitoid interactions, the effects of tree species richness on generality might be more complex and difficult to predict, as we found that neither tree species richness nor tree MPD were related to generality. Thus, our results again indicate that top-down control imposed by higher trophic levels is more significant than bottom-up control with an increasing number of plant species.

Linkage density was positively related to tree species richness, supporting the food web theory, which predicts in our case that network complexity (linkage density) depends on the number of plants (Melián and Bascompte 2002). Although trees were not directly included as a trophic level in our networks, potential network complexity increased with tree species richness, likely enabling higher network stability/resistance (Ebeling et al. 2011). For example, a network might be more sensitive to extinctions because of key species loss due to lower linkage density and lower redundancy (e.g. Naeem and Li 1997). However, parasitoid robustness was unexpectedly negatively related to tree species richness. It was expected that higher trophic levels would be more robust, less influenced by perturbations from lower trophic levels, when plant diversity is higher, as more potential interactions at lower trophic levels should theoretically increase redundancy and resilience of connected higher levels (Blüthgen and Klein 2011, Fornoff et al. 2019). Dilution effects may explain this, as plots with higher richness held fewer individuals of a given tree species. If there are strong prey item (caterpillars, grasshoppers, etc.) preferences for one species, there may be fewer or they may be more densely aggregated and less likely to be encountered by parasitoids. This increased stochasticity in parasitoid wasps could benefit hosts by reducing parasitism pressure overall, weakening top-down controls.

Similar to tree species richness, tree MPD was also positively correlated with vulnerability and linkage density, meaning that the mean number of parasitoids per host species and number of links within the host-parasitoid system can also be promoted by tree MPD, in agreement with several recent studies (Pellissier et al. 2013, Staab et al. 2020, Wang et al. 2020). Our results suggest that the specialization and complexity of higher trophic levels can also be affected by plant phylogenetic diversity. This pattern can be traced to the effects of habitat heterogeneity caused by tree species richness and MPD on higher trophic levels via bottom-up control. The effects of tree MPD were consistent with effects of tree species richness on robustness of parasitoids to host loss. This result suggests that higher trophic levels are sensitive to changes in both plant phylogenetic relatedness and general species dissimilarity via trophic interactions, even the hosts are not all directly interacting with plants, bees excluded. Therefore, it may be that stronger linkages would be found when exclusively exploring such plant-herbivore-parasitoid systems.

Interaction evenness was significantly negatively related to canopy cover, further reinforcing an important role of microclimate (likely temperature and humidity; (Sobek et al. 2009, Fornoff et al. 2021). Our results agree with a previous study on ants, where plant-insect interactions were more even with more open canopies (Dáttilo and Dyer 2014). In our case, canopy cover might change hymenopteran species evenness and then further influence interaction evenness. certain host species tended to nest in plots with higher canopy cover, which might decrease the interaction evenness by favoring parasitoids of fewer, more dominant hosts. This pattern would become more significant when more host and parasitoid species are in a plot, given the positive relationship between higher trophic level diversity and canopy cover.

Future prospects

Overall, our study enables new insights into the dynamics of host-parasitoid interactions under varying environmental conditions, an important step toward building a synthetic model for such biodiversity. A key finding was that although parasitoids and hosts respond to tree species richness, top-down control seems predominant in parasite-host interaction, though whether this holds for others antagonistic interactions requires further investigation.

Different trophic levels and functional groups of species responded differently to lower level changes (Fornoff et al. 2021, Guo et al. 2021). This highlights the complexities of building such models and calls for more studies across habitat types and taxa, to test the generality of our findings. Future studies should also consider the role of host/parasitoid functional traits, because they might play a critical role in modifying network structures and ecosystem functioning.

Materials and methods

Study sites design

This study was conducted in the BEF-China biodiversity experiment, which is the largest tree diversity experiment worldwide. The experiment is located in a subtropical forest near Xingangshan, Jiangxi province, south-east China (29°08′–29°11′N, 117°90′–117°93′E). The mean annual temperature is 16.7°C and mean annual precipitation 1821 mm (Yang et al. 2013). The experiment includes two study sites (Site A and Site B), 4 km apart from each other, that were established in 2009 (Site A) and 2010 (Site B) respectively. A total of 566 plots (25.8×25.8 m) were designed on the two sites, and per plot 400 trees were initially planted in 20 rows and 20 columns with a planting distance of 1.29 m. A tree species richness gradient (1, 2, 4, 8, 16 and 24 species) was established at each site, based on a species pool of 40 local, broadleaved tree species (Bruelheide et al. 2014).

For our study, at both sites (site A and site B) eight plots of each tree species richness level (1, 2, 4, 8) were randomly selected, as well as six and two plots of 16 and 24 mixtures. In addition, at site B eight additional monocultures were sampled (Fornoff et al. 2021), resulting 48 plots in Site B (including 16 monocultures, eight plots for each 2, 4, 8 mixtures and six and two plots of 16 and 24 mixtures. In total, 88 study plots were used (40 plots on Site A and 48 plots on Site B, see Fig. S1).

Sampling

We collected trap nests monthly to sample solitary bees and wasps (Staab et al. 2018) in the 88 plots from September to November in 2015 and April to November in 2016, 2018, 2019 and 2020. For each plot, we installed two poles with trap nests (11 m apart from each other and 9 m away from the nearest adjacent plots) along a SW–NE diagonal centrally per plot (following the design of Ebeling et al. 2012). Each pole stood 1.5 m above ground, and each trap nest consisted of two PVC tubes (length: 22 cm × diameter: 12.5 cm) filled with 75 ± 9 (SD) reed internodes of 20 cm length and diameters varying between 0.1 and 2.0 cm (Staab et al. 2014, Fornoff et al. 2021).

Every month, we sampled the reeds with nesting hymenopterans and replaced them with internodes of the same diameter. All the samples were reared in glass test tubes under ambient conditions until specimens hatched. We identified hatched hosts and parasitoids to species or morphospecies (Supplementary Table S1) based on reference specimens (vouchered at the Institute of Zoology, CAS, Beijing). We were interested in the general patterns of host-parasitoid interactions at the community level, so for the analysis we did not distinguish between the two life-history strategies of parasitoids (true parasitoids and kleptoparasitoids, including hymenopteran and dipteran parasitoids) because they both have the same ecological result, death of host brood cells.

DNA extraction and amplification

All specimens were sequenced for a region of the mitochondrial cytochrome c oxidase subunit I (COI) gene (Hebert et al. 2003). We extracted whole-genomic DNA of hosts and parasitoids using DNeasy Blood & Tissue Kits (QIAGEN GmbH, Hilden, Germany), following the manufacturer’s protocols. COI sequences of samples were amplified using universal primer pairs, LCO1490 (GGTCAACAAATCATAAAGATATTGG) as the forward primer and HCO2198 (TAAACTTCAGGGTGACCAAAAAATCA) or HCOout (CCAGGTAAAATTAAAATATAAACTTC) as the reverse primer. We carried out polymerase chain reactions (PCR) in 96-well plates with 30 μl reactions containing 10 μl ddH2O, 15 μl Premix PrimeSTAR HS (TaKaRa), 1 ul of each primer at 10 μM, and 3 μl template genomic DNA using a thermo cycling profile. The PCR procedure as follows: 94℃ for 1 min; 94℃ for 1 min, 45℃ for 1.5 min and 72℃ for 1.5 min, cycle for 5 times; 94℃ for 2 min, 58℃ for 1.5 min and 72℃ for 1 min, cycle for 36 times; 72℃ for 5 min. We performed all PCRs on an Eppendorf Mastercycler gradient, which were then visualized on a 1% agarose gel. Samples with clean single bands were sequenced after PCR purification using BigDye v3.1 on an ABI 3730xl DNA Analyser (Applied Biosystems).

Sequence alignment and phylogenetic analysis

We applied MAFFT (Misawa, Katoh, Kuma, & Miyata, 2002) to align all sequences, then translated the nucleotides into amino acids via MEGA v7.0 (Kumar, Stecher, & Tamura, 2016) to check for the presence of stop codons with manual adjustments. Host and parasitoid sequences were then aligned against the references using a Perl-based DNA barcode aligner (Chesters, 2019).

We employed two strategies to improve the phylogenetic structure of a DNA barcode phylogeny, which demonstrably improve resulting phylogeny-based diversity indices (Macías-Hernández et al. 2020). These include the integration of 1) molecular sequences of the plot data and 2) phylogenetic relationships from other molecular datasets. Integration was achieved following Wang et al. (2020) and Chesters (2020): reference DNA barcodes of Hymenoptera and Diptera were downloaded from the BOLD API (www.boldsystems.org/index.php/API_Public), which were variously processed (e.g. to retain only fully taxonomically labelled barcodes, to remove low quality or mislabeled entries, and to dereplicate to a single exemplar per species), and then aligned (Chesters 2019). A single outgroup was included for which we selected the most appropriate insect order sister to Diptera and Hymenoptera (Misof et al. 2014), a representative of the order Psocoptera (Psocidae, Psocus leidyi). We then constructed a phylogeny of the references and subjects, with references constrained according to the method described earlier (Chesters 2020). A number of backbone topologies were integrated for setting hard and soft constraints, including a transcriptomics-derived topology (Chesters 2020), a mitogenome tree of insects (Chesters 2017), Diptera-specific trees (Wiegmann et al. 2011, Cranston et al. 2012, Ament 2017) and Hymenoptera-specific trees (Peters et al. 2011, Branstetter et al. 2017, Cardinal 2018). The constrained inference was conducted with RaxML version 8 (Stamatakis 2014) under the standard GTRGAMMA DNA model with 24 rate categories. According to the backbone trees used, most taxa present were monophyletic with a notable exception of Crabronidae, for which there is emerging phylogenomic evidence of its polyphyly (Sann et al. 2018).

Tree phylogenetic diversity, functional diversity, and environmental covariates

The phylogenetic diversity of the tree communities was quantified by calculating wood volume-weighted phylogenetic Mean Pairwise Distance (MPD) (Tucker et al. 2017). Tree wood volume was estimated from data on basal area and tree height (Bongers et al. 2021) measured in the center of each plot. Moreover, to represent variations towards the tips of the phylogeny beyond MPD, we additionally calculated Mean Nearest Taxon Distance (MNTD), which is a measure that quantifies the distance between each species and its nearest neighbor on the phylogenetic tree (Webb 2000). Phylogenetic metrics of trees were calculated based on a maximum likelihood phylogenetic tree available for the tree species in our study area (Purschke et al. 2017). Considering that predatory wasps mainly feed on herbivorous caterpillars, we calculated tree functional diversity to test the indirect effects on hymenopteran communities and relevant network indices. Specifically, seven leaf traits were expressed as Rao’s Q (Ricotta and Moretti 2011), including specific leaf area, leaf toughness, leaf dry matter content, leaf carbon content, ratio of leaf carbon to nitrogen, leaf magnesium content, and leaf calcium content. These functional traits were commonly related to higher trophic levels in our study area (Wang et al. 2020, Chen et al. 2023), which are the main food resources of our predatory wasps. All of the traits were measured on pooled samples of sun-exposed leaves of a minimum of five tree individuals per species following standard protocols (Pérez-Harguindeguy et al. 2003).

As our analyses mainly compare community patterns among study plots, we additionally considered potential effects of environmental variation by using plot means of slope, elevation, “eastness” (sine-transformed radian values of aspect), and “northness” (cosine-transformed radian values of aspect) as environmental covariates that characterize the heterogeneity of the study plots. We also accounted for the potential effects of canopy cover at plot level for host-parasitoid interactions, as it can structure hymenopteran communities (Perlík et al. 2023). Canopy cover was calculated as in Fornoff et al. (2021) based on above-nest images.

Statistical analysis

All analyses were conducted in R 4.1.2 with the packages ape, vegan, picante, bipartite, and caper (http://www.R-project.org). Prior to analysis, samples from the five years (2015, 2016, 2018, 2019, and 2020) were pooled at the plot level to discern overall and generalizable effects permeating this system. We excluded three plots with no living trees because of high mortality, resulting in 85 plots in the final analysis.

Composition of trees, hosts and parasitoids

The species and phylogenetic composition of trees, hosts, and parasitoids were quantified at each plot with nonmetric multidimensional scaling (NMDS) analysis based on Morisita-Horn distances. The phylogenetic composition was calculated by mean pairwise distance among the host or parasitoid communities per plot with the R package “picante.” To test the influence of study plot heterogeneity on these relationships, we fitted their standardized values to the ordination on the basis of a regression with the NMDS axis scores (Quinn and Keough 2002). For the analysis, we considered tree species richness, tree functional and phylogenetic diversity, canopy cover, and environmental covariates (elevation, eastness, northness, and slope) as plot characteristics. We assessed the significance of correlations with permutation tests (permutation: n= 999).

Phylogenetic match of hosts and parasitoids

In addition, we used a parafit test (9,999 permutations) with the R package “ape” to test whether the associations were non-random between hosts and parasitoids. The species that were not attacked by parasitoids or failed to generate sequences were excluded from the analyses. For species abundance and composition, see Table S1.

Host-parasitoid interactions

We constructed quantitative host-parasitoid networks with the R package “bipartite” for the level of host and parasitoid communities across for each plot of the two sites. All hosts were grouped because there were too few abundant bee species to separate them out from other hosts. We calculated six indices to quantitatively characterize the structure of the interaction networks: weighted generality (effective number of host species per parasitoid species), weighted vulnerability (Effective number of parasitoid species attacking a host species), robustness (degree of network stability), linkage density (degree of network specialization), and interaction evenness (degree of network evenness). Generality was defined as the weighted mean number of host species per parasitoid species, , with A_j being the number of interactions of parasitoid species j, m the total number of interactions of all species, and H_j the Shannon diversity of interactions of species j. Vulnerability was the weighted mean number of parasitoid species per host species, (Bersier et al. 2002). Robustness was defined as the area under the extinction curve, reflecting the degree of decreases of one trophic level with the elimination of the other trophic levels, here using the robustness index for higher trophic levels (i.e. parasitoids). For linkage density, , we used the realized proportion of possible links between the two trophic levels as the mean number of interactions per species across the entire network (Tylianakis et al. 2007). Interaction evenness was defined as E_S = − ∑_i ∑_j p_ijlnp_ij /lnIJ, which is used to describe Shannon’s evenness of network interactions (Dormann et al. 2009). To check whether all network indices significantly differ from chance across all study plots, we used Patefield null models (Dormann et al. 2009) to compare observed indices with simulated values (10,000 times).

Linear mix-effects models

To test the effects of tree species richness, tree phylogenetic, and functional diversity, as well as canopy cover and the other environmental covariates (including slope, elevation, eastness, and northness) on the six network indices, we used linear mixed-effect models. For our analyses, the study sites were considered as random effects, and the others were treated as fixed effects (see above). Given the strong correlation between tree species richness and tree MPD (Pearson’s r = 0.74, p<0.01), we excluded tree MPD in the models where tree species richness was a predictor. To evaluate the potential effects caused by tree MPD, we also ran alternative models where tree species richness was replaced with tree MPD. We simplified all models by gradually removing non-significant factors to obtain the most parsimonious model with the lowest AICc. To ensure that the analyses were not strongly affected by multicollinearity, the correlations among all predictors were tested (Fig. S2), and variance inflation factors (VIF) of our statistical models were checked.

Acknowledgements

We thank the BEF-China consortium for support (especially, Bo Yang). We thank Mr. Yinquan Qi for his help with the collection. This work was supported by the National Key Research Development Program of China (2022YFF0802300), the National Science Fund for Excellent Young Scholars (32122016), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB310304), the National Natural Science Foundation, China (32100343, 32070465), and the National Science Fund for Distinguished Young Scholars (31625024). MQW was supported by the Alexander von Humboldt research fellowships. CDZ’s lab is funded by the Key Program of the National Natural Science Foundation of China (No. 32330013) and also has been continuously supported by grants from the Key Laboratory of the Zoological Systematics and Evolution of the Chinese Academy of Sciences (grant number 2008DP173354).