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 of higher trophic levels usually coincide with variations in their diversity and community composition, 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 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 coexistence and its role for biodiversity conservation, we must further study the mechanisms on dynamics of networks from multiple perspectives (Tittensor et al. 2014, Brondizio et al. 2019). Interactions can be viewed from a top-down or bottom-up perspective. Previous studies have shown there might be asymmetric effects across trophic levels (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 also driven by microclimate, which could be determined by plant structuring, such as canopy cover (Fornoff et al. 2021, Perlík et al. 2023). Understanding how canopy cover responds 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 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 their host fluctuations may cascade up to impact the parasitoids’ presence or absence. Therefore, insect host-parasitoid systems are ideal for studying the relationships between community-level changes and species interactions (Jeffs and Lewis 2013). Previous studies mainly focused on the influence of abiotic factors on host-parasitoid interactions, 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 (i.e. taxonomic, functional and phylogenetic diversity) in modifying host-parasitoid interaction networks, as key biotic determinants of overall ecosystem structure, remains relatively little explored (but see Staab et al. 2016). Recent studies mainly focused on basic diversity associations between host and parasitoids (Ebeling et al. 2012, Schuldt et al. 2019, Guo et al. 2021). These studies have demonstrated both direct and indirect effects (i.e. one pathway and more pathways via other variables) of plant diversity on both hosts and parasitoids diversity, 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. The other diversity components (e.g. 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 (e.g. phylogenetic diversity or phylogenetic congruence) 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 (e.g. taxonomic, phylogenetic, functional) and relevant trophic interactions (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 (De Frenne et al. 2021, Popkin 2021). Understanding how multiple dimensions of biodiversity 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 discern the primary components of the diversity and composition of tree communities that affect higher trophic level interactions via quantifying the strength and complexity of associations between hosts and parasitoid. We expected that (a) multiple tree community metrics, such as species, phylogenetic and functional diversity and species community composition can structure host and parasitoid community compositions, especially via phylogenetic processes (e.g. lineages of trophic levels diverge and evolve over time), 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, mean abundance and richness are 400 and 45, respectively, for each plot) and 50 parasitoid species (38 Hymenoptera and 12 Diptera, mean abundance and richness are 14 and 9, respectively, for each plot). The full species list and their abundances are given in Table S1. Overall, our sampling was adequate for analysis (especially for hosts), as confirmed by the sampling completness evaluation (Fig. S1).

Community composition of hosts and parasitoids

Host species composition was significantly related to the species and phylogenetic composition of the trees and parasitoid communities (NMDS axis scores; i.e. preserving the rank order of pairwise dissimilarities between samples), as well as to canopy cover, tree mean phylogenetic distance (MPD), elevation, and eastness (sine-transformed radian values of aspect) (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 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). The PERMANOVA analysis also highlighted the important role of canopy cover for host and parasitoid community (Table S6-9). The Mantel test revealed a consistent pattern with the NMDS analysis, highlighting a pronounced relationship between the species composition of hosts and parasitoids (Table S10). However, the correlation between the phylogenetic composition of hosts and parasitoids was not significant.

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) for testing the hypothesis of coevolution between a clade of hosts and a clade of parasites, suggesting nonrandom associations in the phylogenetic structure of parasitoid and host communities (Fig. 2, Fig. S4).

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). Parasitoid generality was only marginally associated with canopy cover, and was not related to tree species richness or the other environmental variables. In the alternative models (tree species richness replaced by tree MPD), host vulnerability and linkage density were significantly positively related to tree MPD (Fig. 4a, 4b; Table S11), while robustness of parasitoids was negatively related to tree MPD (Fig. S5, Table S11). The results of other network metrics (parasitoid generality and interaction evenness) were consistent with those of the primary models. Tree mean nearest taxon distance (MNTD) was unrelated to any network indices.

Fig. 3.

Fig. 4.

Table 2.

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 parasitoid generality, host vulnerability, and linkage density, respectively; all values for robustness, and interaction evenness), strongly demonstrating that interactions between species were not driven by random processes.

Discussion

Our study demonstrates that tree species richness and phylogenetic diversity play key roles in modulating interacting communities of hosts and parasitoids. These interactions are further structured by the phylogenetic associations between hosts and parasitoids. Moreover, canopy cover partly determined host-parasitoid network indices patterns, including host vulnerability, linkage density, and interaction evenness. This result supports a recent finding that the structure of host-parasitoid networks is also mediated by changes in microclimate, which is directly related to canopy cover, to some extent (Fornoff et al. 2021). These patterns were highly associated with multiple tree diversity metrics (tree species, phylogenetic and functional diversity), and compositional changes which are key to understand how host-parasitoid interactions may be impacted by biodiversity loss of lower trophic levels in food webs through trait- and phylogeny-based processes.

Community composition associations

Host species composition was influenced by several factors, including the species and phylogenetic composition of trees and parasitoids, tree diversity (species richness and MPD), and other biotic factors such as elevation and slope (Fig. 1a; Table 1). 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 resource heterogeneity increasing with 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 (Fig. 1b; Table 1), 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 (Fig. 1c, d; Table 1). 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. 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 organisms highly associated with tree phylogenetic composition (e.g. caterpillars) might indirectly determine predatory wasp (host) phylogenetic composition, as recently found for interactions between plants, caterpillars and spiders (Chen et al. 2023). This could be further tested by analyzing 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 indirectly 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; Fig. 2) between the phylogenetic trees of the host and their parasitoids (see also Peralta et al. 2015). Such pattern could be further confirmed by the significant association between host phylogenetic composition and parasitoid phylogenetic composition (Fig. 1c), which suggested that their interactions are phylogenetically structured to some extent. However, this significant pattern was observed only in the NMDS analysis and not in the Mantel test, suggesting that the non-random interactions between hosts and parasitoids could not be simply predicted by their community similarity and associations between the phylogenetic composition of hosts and parasitoids are more complex and require further investigation in the future.

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 to microclimate (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 species richness did not significantly influence the diversity of hosts targeted by parasitoids (parasitoid generality), but caused a significant increase in the diversity of parasitoids per host species (host vulnerability) (Fig. 3a; Table 2). 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). Such positive relationship between host vulnerability and tree species richness suggested that host-parasitoid interactions could be driven through bottom-up effects via benefit from tree diversity. For example, parasitoid species increases more than host diversity with increasing tree species richness (Guo et al. 2021), resulting in an increase of host 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 parasitoid generality might be more complex and difficult to predict, as we found that neither tree species richness nor tree MPD were related to parasitoid generality.

Linkage density was positively related to tree species richness (Fig. 3b; Table 2), supporting the food web theory, which predicts in our case that network complexity (linkage density) depends on the number of plants (Blüthgen and Staab 2024). 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, Staab et al. 2015). For example, a network might be more sensitive to extinctions because of key species loss due to lower linkage density and lower redundancy (Blüthgen and Staab 2024). However, parasitoid robustness was unexpectedly negatively related to tree species richness (Fig. 3c; Table 2). 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 individuals per area 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 host vulnerability and linkage density (Fig. 4a, b), 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. This is in agreement with several recent studies showing that plant phylogenetic diversity not only affects herbivores but also higher tropic levels (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(Fig. 4c), further reinforcing the important role of canopy-cover modulated microclimate (likely temperature and humidity) for trophic interactions (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 experimental changes in plant communities (Fornoff et al. 2021, Guo et al. 2021). This highlights the complexities of building multi-trophic networks 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). The tree species pools of the two plots are nonoverlapping (16 species for each site). The composition of tree species within the study plots is based on a “broken-stick” design (see 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 in 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. S2).

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 (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 cmL×Ldiameter: 12.5 cm) filled with 75L±L9 (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. We evaluated our sampling completeness with r package “iNEXT” (Hsieh et al. 2016).

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 (tree MPD and MNTD) 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 used for calculation of tree functional diversity, which was calculated as the mean pairwise distance in trait values among tree species, weighted by tree wood volume, and 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, such as herbivores and predators (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 hemispherical photographs.

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 with nonmetric multidimensional scaling (NMDS) analysis based on Morisita-Horn distances. The minimum number of required dimensions in the NMDS based on the reduction in stress value was determined in the analysis (k = 2 in our case). We centered the results to acquire maximum variance on the first dimension, and used the principal components rotation in the analysis. The phylogenetic composition was calculated by mean pairwise distance among the host or parasitoid communities per plot with the R package “picante” applying the ‘mpd’ function. To test the influence of study plot heterogeneity on these relationships, we fitted their standardized values (see vectors in Table S2) to the ordination on the basis of a regression with the NMDS axis scores (Quinn and Keough 2002). NMDS was widely used to summarize the variation in species composition across plots. The two axes extracted from the NMDS represent gradients in community composition, where each axis reflects a subset of the compositional variation. These axes should not be interpreted in isolation, as the overall species composition is co-determined by their combined variation. For clarity, results were interpreted based on the relationships of variables with the compositional gradients captured by both axes together. 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). To strengthen the robustness of our findings based on NMDS, we further validated the composition results using Mantel test and PERMANOVA (with ‘adonis2’) for correlation between communities and relationships between communities and environmental variables.

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. This is widely used to assess host-parasite co-phylogeny by analyzing the congruence between host and parasite phylogenies using a distance-based matrix approach. 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 at community level with the R package “bipartite” for each plot of the two sites. Bees and wasps were considered together as hosts because there were too few abundant bee species to analyze separate interaction networks for bees and wasps. We calculated five indices to quantitatively characterize the structure of the interaction networks (Blüthgen and Staab 2024): weighted parasitoid generality (effective number of host species per parasitoid species), weighted host 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). Parasitoid generality was defined as the weighted mean number of host species per parasitoid species, with A 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. Host vulnerability was the weighted mean number of parasitoid species per host species, Vulnerability = (Bersier et al. 2002). i=l _m Robustness was defined as the area under the extinction curve, reflecting the degree of decreases of one trophic level with the random elimination species of the other trophic levels, here using the robustness index for higher trophic levels (i.e. parasitoids). The robustness calculation was also weighted based on interaction abundance. For linkage density, L_q = 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∑_jp_ij lnp_ij/ln/J, 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 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 other variables 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. S3), and variance inflation factors (VIF) of our statistical models were checked.

Data availability

Data is available on Science Data Bank at https://www.scidb.cn/en/s/Ujq22u, and will be available on the BEF-China project database at https://data.botanik.uni-halle.de/bef-china/datasets.

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

Additional information

Author contributions

C.D.Z., A.L. and M.Q.W. conceived the idea for the manuscript; C.D.Z., A.M.K., S.K.G., P.F.G. and M.Q.W. designed the study. S.K.G., P.F.G., J.J.Y., G.A.C., Z.Q.N., M.S., J.T.C., Y.L., Q.S.Z., F.F., X.S., S.L., M.M., A.M.K., A.S., X.L., K.M. and H.B. collected and/or contributed data; D.C., M.Q.W. and S.K.G. conducted the phylogenetic analysis; M.Q.W. S.KG. and P.F.G. conducted the statistical analyses and wrote the manuscript, with input from F.F., A.S., M.S. and M.M and all coauthors.

Additional files

Supplementary tables and figures