Tree diversity and functional leaf traits drive herbivore‐associated microbiomes in subtropical China

Abstract Herbivorous insects acquire microorganisms from host plants or soil, but it remains unclear how the diversity and functional composition of host plants contribute to structuring herbivore microbiomes. Within a controlled tree diversity setting, we used DNA metabarcoding of 16S rRNA to assess the contribution of Lepidoptera species and their local environment (particularly, tree diversity, host tree species, and leaf traits) to the composition of associated bacterial communities. In total, we obtained 7,909 bacterial OTUs from 634 caterpillar individuals comprising 146 species. Tree diversity was found to drive the diversity of caterpillar‐associated bacteria both directly and indirectly via effects on caterpillar communities, and tree diversity was a stronger predictor of bacterial diversity than diversity of caterpillars. Leaf toughness and dry matter content were important traits of the host plant determining bacterial species composition, while leaf calcium and potassium concentration influenced bacterial richness. Our study reveals previously unknown linkages between trees and their characteristics, herbivore insects, and their associated microbes, which contributes to developing a more nuanced understanding of functional dependencies between herbivores and their environment, and has implications for the consequences of plant diversity loss for trophic interactions.


| INTRODUC TI ON
Insect symbionts comprise of bacteria, fungi, and viruses, and persist both in the insects and on the cuticle of their exoskeleton and engage in a variety of interactions (Frago et al., 2012;Klepzig & Six, 2004).
Herbivorous insects can acquire specific symbiont bacterial species from host plants or environment, such as the soil in which the host plants grow (Kikuchi et al., 2012;Sugio et al., 2015). Many symbionts found in the gut of insects are allied to microorganisms of the immediate environment (Frago et al., 2012). In addition to gut symbionts, herbivorous insects are known to acquire intra-and extracellular microbes also from host plants (Caspi et al., 2011;Li et al., 2016).
Studies of bacterial communities of Lepidoptera have increased rapidly in past several years, mostly because Lepidoptera possess extraordinary species richness and are important forest and agricultural pests (Belda et al., 2011;Broderick et al., 2004;Gayatri et al., 2012;Xia et al., 2013;Xiang et al., 2006). However, some recent studies reported that caterpillars lack resident gut microbes (Hammer et al., 2017;Hammer et al., 2019). Most of the microbes found in the Lepidoptera gut were found to be shared with the surrounding leaf surface or with the soil in which the host plant grew, and the bacteria seem to have no effects on growth and survival of the caterpillar (Hammer et al., 2017;Whitaker et al., 2016). If indeed lacking a persistent microbiome, then any microbial-driven processes in caterpillars might be more susceptible to environmental influences than would otherwise be.
What kind of environmental factors influence the microbial community of herbivorous insects? And to what extent did these factors influence the herbivore microbes? To answer these questions, we analyze the relationships between the bacterial composition and diversity of host plant and herbivores. This is conducted in the "BEF-China" experiment, a large-scale forest biodiversity experiment incorporating random extinction scenarios of tree species, and used to estimate the ecological effects of biodiversity loss (Bruelheide et al., 2014). We target microbial symbionts on and in the caterpillar body, using 16S rRNA gene sequencing ( Figure 1). We hypothesize that the community composition of herbivore-associated microbes is driven by specific aspects of their surroundings (i.e., leaf traits, host tree species identity, and diversity of tree species in the host tree community) and by the composition and diversity of the host herbivore species. We aim to assess the relative importance of direct host-mediated effects versus environmentally mediated effects on the herbivore microbiome.

K E Y W O R D S
16S rRNA, Bacteria, BEF-China, herbivore-associated microbiome, leaf characteristics, Lepidoptera, Microbial ecology F I G U R E 1 Overview of the study; (a) location of the study site; Xin-Gang mountain, Jiangxi Province (29°08′-29°11′N, 117°90′-117°93′E), southeast China, with a typical subtropical climate; (b) two example plots in the study site, with tree species richness of 2 & 4; (c) presence / absence of three Lepidoptera species in the two plots; (d) the relationships between trees, lepidopteran samples and their associated bacterial OTUs 2 | MATERIAL S AND ME THODS

| Experimental design for study sites
The study was conducted in the BEF-China forest biodiversity experiment, which was established in the southeast of China (Xingangshan, Jiangxi Province, 29°08′-29°11′N, 117°90′-117°93′E) in 2009. The study area has a subtropical monsoon climate with an annual mean temperature of 16.7°C and precipitation of 1,821 mm (Yang et al., 2013). The 38.4 ha study area consists of two sites including a total of 566 plots (25.8 *25.8 m/ plot, 271 plots in site "A" and 295 plots in site "B"). In each plot, 400 trees were planted in 20 rows and 20 columns. The species pool includes 40 species of trees. Species were selected for each plot according to a random broken stick design for extinction scenarios of

| Sample collection of individual caterpillars
We focus on the larval stage, being the primary feeding stage of Lepidoptera and typically the focus for studies of herbivory. The collection of lepidopteran larvae used herein has been previously reported, with focus on tree diversity effects on herbivores themselves in Wang et al. (2019Wang et al. ( , 2020; thus, here we focus on the microbiomes associated with these herbivores. In October 2018, we selected individuals from the caterpillar samples for extraction of bacterial DNA. We recognize that DNA contamination could be an issue. We therefore took steps to avoid this. Firstly, the caterpillar samples were cleaned with sterile water and 75% ethanol and then kept in separate tubes with 99.5% ethanol. Secondly, we extracted whole lepidopteran larvae, but we broke their bodies sufficiently that we can obtain their gut microbiota. All individuals were stored in a − 20 ℃ freezer prior to DNA extraction. To ensure the comparability of caterpillar-associated bacteria across the tree diversity levels, we selected the individual caterpillars randomly from different plots based on the BEF-China design (118 caterpillars came from monocultures, and 110, 104, 97, 98, and 107 came from the mixtures of 2, 4, 8, 16, and 24 species, respectively; details see from Table 1). Those Lepidoptera species that had the highest abundances in each plot were selected as representative for the plot. As the number of caterpillars among plots and that of different Lepidoptera species varied greatly, we selected the caterpillars to sample as many different tree species per plot as possible. Because of tree mortality, the whole sample included 37 tree species from 54 plots (Table 1). For each plot, all bacterial sequences obtained from the selected caterpillars served as the bacterial community of the Lepidoptera larvae from the given plot. Totally, this resulted in the selection of 444 trees covering the full tree diversity gradient and a total of 634 caterpillar individuals.
As the current strategy for sampling caterpillar-associated bacteria at the plot level has some limitations (particularly, the bacterial composition of some rare Lepidoptera species was not considered), we tested whether the results were affected by sample size. To this end, we used additional linear models to check the relationships between bacterial and tree species richness at the tree richness level (Figure 3a and 3b). We found that the different Lepidoptera species feeding on the same tree species showed similar bacterial communities (relative abundance of bacterial phyla; Figure S1), and there was no significant difference in bacterial diversity (results not shown).

| DNA extraction, amplification, quantitation, and sequencing
Total DNA was extracted from the individual caterpillars (because many larvae were very small or less than 5 mm) using Qiagen DNeasy Tissue Kit (QIAGEN GmbH, Hilden, Germany), following the manufacturer's protocol. Samples were processed using sterile tools and conditions. The DNA extracts were quantified using the Qubit 4.0 Fluorometer and stored at −20°C for further processing.

| Bioinformatics analyses
We used both VSEARCH v 2.8.1 (Rognes et al., 2016) and USEARCH v 11 (Edgar, 2010) to process raw sequences. Firstly, the read pairs were merged, primers were trimmed, and quality filtering excluded short and low-quality reads (below 25), using VSEARCH.
Then, we dereplicated the data (retaining abundance information) TA B L E 1 Tree species richness, composition, and caterpillar sample size of the study plots the OTUs were assigned taxonomic information using the Silva 138 SSU database by using the "sintax" command in USEARCH (threshold ≥ 90%), after which nonbacteria OTUs were removed (Pruesse et al., 2007). Bacteria were classified at the level of phylum, class, order, family, genus, and species.
Phylogenetic diversity (PD) of bacteria was incorporated as response variables. To construct the bacterial phylogeny, we used MAFFT v 7.0 (Misawa et al., 2002) to align the sequences, trimmed the alignment with MEGA v 7.0 (Kumar et al., 2016), and inferred the phylogeny using the ML software IQ-TREE v 1 (Nguyen et al., 2015).

| Leaf traits
We selected 11 morphological and chemical leaf traits which we considered prime candidates for determining leaf quality for herbivore insects, and to characterize plot conditions in accordance with nutritional quality and potential defense traits of the trees. We used leaf area, specific area, dry matter content and toughness as the main morphological traits, and leaf potassium, calcium, magnesium, sodium, phosphorus, carbon, and nitrogen content, and C: N ratio, as the chemical leaf traits (see Table S1 for abbreviations that will be used below). All of these traits were measured on sun-exposed, fully expanded, undamaged leaves from five to seven individuals per tree species according to standard protocols Pérez-Harguindeguy et al., 2003). More details on trait measurements can be found in Kröber et al. (2012).

| Community-weighted mean trait values, functional and phylogenetic diversity
We used the community-weighted mean (CWM) of each trait as well as the functional diversity of selected traits for each tree species, which is the mean value of each species' trait weighted by the species contribution to the plot wood volume. The CWM values of each trait in each plot were calculated by the following equation: is the relative tree wood volume of species i in plot p and t i is the mean trait value of species i (Garnier et al., 2004). Tree wood volume was estimated from basal area and tree height measured on trees in the center of each plot according to Fichtner et al. (2017). We used species-mean trait values as previous studies in BEF-China demonstrated that variability in trait-environment relationships was much more pronounced at the interspecific than the intraspecific level (Schuldt et al., 2012).
To characterize the plot conditions of the study sites, the following metrics were calculated at the plot level. The functional diversity of trees was calculated by the mean pairwise distance of trait values among tree species, weighted by relative wood volume, and expressed as Rao's Q (Ricotta & Moretti, 2011). We also depicted PD of tree communities by wood volume-weighted phylogenetic mean pairwise distance (MPD), which in the abundance-weighted case is equivalent to Rao's Q (Tucker et al., 2016). In addition, we calculated the mean nearest taxon distance (MNTD), a measure of the phylogenetic distance to the nearest taxon, which for each taxon quantifies the extent of terminal clustering (Webb, 2000;Webb et al., 2002).
Phylogenetic indices were calculated on a maximum likelihood phylogenetic tree of all woody species recorded in all plots (Purschke et al., 2017).

| Statistical analyses
Statistical analyses were conducted using the packages picante (Kembel et al., 2010), vegan (Oksanen, 2008), ape (Paradis & Schliep, 2019), edgeR (Robinson et al., 2010), phyloseq (McMurdie & Holmes, 2013), lavaan (Rosseel, 2012), and lulu (Frøslev et al., 2017) in R v 3.5.2 (http://www.R-proje ct.org). Firstly, to eliminate the impacts on differing read numbers across samples, the number of sequences of all samples was rarefied to the lowest read number by using the "rarefy" function of the vegan package (rarify depth: 11,000). The bacterial S obs (observed bacterial richness), Chao1 (nonparametric estimator for bacterial richness), Shannon diversity, and Pielou's evenness were calculated for each plot from bacterial abundance using the "diversity" function of the vegan package. To improve normality and variance in homogeneity of the model residuals, tree species richness, Lepidoptera richness, abundance, and bacterial S obs were log-transformed, and Chao1, Shannon diversity, and Pielou's evenness were square-root transformed. All continuous predictors were standardized before the analyses.
To avoid multicollinearity affecting our statistical analyses, we Path analyses were conducted to explore the potential causal relationships among tree species richness, Lepidoptera richness, leaf traits, plot covariables, and richness of caterpillar-associated bacteria. Based on prior and theoretical knowledge, we hypothesized that species richness of trees and Lepidoptera might directly influence bacterial communities. In addition, tree species richness could also indirectly influence bacterial communities via Lepidoptera richness.
Concurrently, leaf traits may have direct or indirect effects on bacterial communities. Nonsignificant pathways were gradually removed if their removal improved model fit (Scherber et al., 2010). The model fit was assessed by comparative p value, fit index value (CFI), Akaike Information Criteria (AIC), and root mean square errors of approximation (RMSEA). Adequate model fits are indicated by high CFI, low AIC, and low RMSEA.
As for bacterial Beta diversity, we defined it as Arrhenius exponent "z" (calculated with "betadiver" function) in our study and compared the variances of dissimilarities using the "adonis" function in vegan (Anderson, 2006), with a p value obtained from 9,999 per-

| RE SULTS
After merging reads and filtering for quality, we obtained 16,490,248 reads (out of 20,466,051 raw reads) in total from 634 caterpillar individuals, delineated to 7,909 bacterial OTUs. There were minor shifts in bacterial phyla among tree genera of plots ( Figure 2 and Table S2).
The three most abundant bacterial phyla were Proteobacteria, Firmicutes, and Actinobacteria. Among classes, Alpha-and Gammaproteobacteria classes predominated. We found six core OTUs (present in more than 90% of the samples), which belong to four phyla (Proteobacteria, Actinobacteria, Firmicutes, Deinococcus-Thermus).
Based on information of the taxonomic assignment, two of them belong to the genus Enterococcus, one was Enterococcus sp. and the other was Enterococcus faecalis, and the remains were common species exist in the environment. We chose the four most abundant Lepidoptera species (> 24 individuals and observed at all tree diversity levels) to examine distribution patterns of their bacterial phyla ( Figure S2). We found that the bacteria of a given Lepidoptera species varied greatly in composition among tree diversity levels, and the bacteria of different Lepidoptera species often varied substantially within a given tree diversity level. Further, based on the assigned taxonomic information, we inferred that only two of the core species are likely to have been derived from the gut of caterpillars, the remainder likely from the leaf surface or elsewhere in the environment.
Bacterial richness (both S obs and Chao1) and Shannon diversity differed across the study plots, and both were significantly correlated with tree species richness (Figure 3c-3f, Table 2). These results were confirmed when analyzing the relationship between tree species richness and bacterial diversity at the tree richness level (i.e., based on similar numbers of caterpillars in each richness level; Figure 3a and 3b). Furthermore, estimated bacterial richness (Chao1) at the plot level was also affected by Lepidoptera abundance, Lepidoptera richness, and the interaction of tree species richness and Lepidoptera richness (Table 2). Tree species richness also corresponded to the caterpillar microbiomes when using Pielou's evenness. Moreover, bacterial richness (both S obs and Chao1) were positively correlated with CWMs of several leaf traits of tree communities, especially LDMC and leaf K content. Shannon diversity was correlated with tree species richness and the CWMs of LT, leaf Ca content, and K content (Table 2). Pielou's evenness was also positively correlated with CWMs of Ca content and K content.
The path analyses (Figure 4; Table S3) showed that tree species richness directly influenced Lepidoptera richness, which in turn af- The analyses of differentially abundant bacterial OTUs between tree richness levels were conducted by fitting a generalized linear model with a negative binomial distribution to normalized values for each of the 7,909 bacterial OTUs, and testing for differential abundance using a likelihood ratio test. We first used bacterial counts from monoculture plots as a control and an adjusted P value cutoff of 0.01, and compared it with the bacterial counts from 2, 4, 8, 16, and 24 species mixed plots separately. As shown in Figure 6, the F I G U R E 2 Bacterial phyla present in each tree genus, with relative abundances averaged across tree individuals of the same genus. The bacterial phyla are listed in the legend. Analysis is limited to phyla with site relative abundance >= 0.1% enriched bacterial counts were always greater than depleted counts in comparison with the control. Then, we used the counts from 2, 4, 8, and 16 tree species mixtures as a control and compared successively with higher diversity mixtures. We found that the counts of the enriched species were higher than that of depleted species with increasing diversity when using monocultures and 2 species mixtures as controls. Although when using the 4-species mixtures as a control, the counts of depleted bacteria exceeded the enriched, and the counts of both reduced significantly. In addition, Figure S3 shows the numbers of differentially enriched and depleted bacterial OTUs between each tree richness level compared with different controls.

| D ISCUSS I ON
This study highlights the impacts of tree diversity on the diversity and community composition of herbivore-associated bacteria, and shows they are influenced by tree diversity and characteristics of the leaf, in what is both a direct and indirect interaction. The direct effect of tree diversity on bacterial diversity was found to predominate, whereas the composition of bacterial communities was to a large part determined by tree diversity and leaf functional traits, especially LDMC and LT, but also chemical leaf traits such as calcium and potassium concentrations. Considering that there is an increasing number of studies reporting that the larva of Lepidoptera recruit microbes from the environment, and that they lack a persistent gut microbiome (Hammer et al., 2017;Hammer et al., 2019), these are key findings that help to better understand which environmental factors determine these microbial communities and how such environmental effects may influence herbivore functioning.

TA B L E 2 (Continued)
Tree diversity affected the diversity of caterpillar-associated bacteria through influencing the abundance and diversity of lepidopteran larvae. Wang et al. (2019) reported that the impact of tree diversity on herbivore diversity is generally indirect, as tree diversity had strong effects on herbivore abundances, which in turn can affect herbivore diversity. The increase in bacterial diversity that follows increasing Lepidoptera diversity at the plot level is consistent with the expectation that more Lepidoptera individuals and species provide more niche opportunities for bacteria (Akiko et al., 2015).
Thus, diversity at one trophic level begets biodiversity at other trophic levels. Moreover, tree diversity was found to also directly influence the diversity of caterpillar-associated bacteria. The most common bacterial groups of the phyllosphere are Acidobacteria, Actinomycetes, Bacteroidetes, Firmicutes, and Proteobacteria (Bodenhausen et al., 2013;Bulgarelli et al., 2013), the latter both the most abundant taxonomic group observed in our study as well as generally associated with the phyllosphere reported by others (Humphrey et al., 2014). This suggests that the phyllosphere is one of the main sources of the herbivore microbiome (Hammer et al., 2017;Whitaker et al., 2016). Kembel et al. (2014) reported that phyllosphere microbial composition differs according to position and height of tree leaves, and physiological and biochemical features such as water content, leaf mass, nitrogen and phosphorus concentrations, leaf surface structure, and thickness. Moreover, both bacterial and fungal communities of the phyllosphere are seasonally dynamic (Jumpponen & Jones, 2010;Rastogi et al., 2012). We suspect that diversity in environmental drivers is one of the main reasons for the very high bacterial compositional dissimilarity observed between caterpillars in this study site, with merely 6 bacteria OTU of 7,909 that were commonly observed. Further, taxonomic analysis showed It is important to note the effect of leaf traits on the diversity and distribution of caterpillar-associated bacteria. As mentioned above, the phyllosphere appears to be a key source for herbivore microbiomes and is moderated by tree characteristics such as leaf structure (also, LDMC is directly affected by leaf thickness, structure, and specific leaf area and reflects the ability of plants to obtain resources). LDMC and LT are usually expected to negatively associate with herbivory because structurally robust leaves are relatively difficult to consume (Pérez et al., 2003).
However, both the results herein and some previous reports suggest a positive relationship between LDMC and leaf herbivory (Lepidoptera richness in this study; Schuldt et al., 2012), probably because there are herbivores specifically adapted to tough leaves (Pérez et al., 2003) and herbivores have to consume more of less nutritious foliage to gain the same nitrogen accumulation rates (Scriber & Slansky, 1981).
In addition to the two leaf traits (LDMC and LT) mentioned above, we also found that leaf potassium (K) content, calcium ( There is, however, considerable variation in mineral requirements of herbivore insects.
Tree species richness was found to be an important factor that affected caterpillar-associated bacteria community composition.
A remarkable result was that certain bacterial OTUs were more abundant in tree species mixtures compared to monoculture plots.
However, the accumulation rate of bacterial taxa in more speciesrich mixtures gradually decreased. That the increase of tree diversity might have a certain stabilizing effect on the herbivore-associated bacterial community was also supported by our finding that the bacterial species composition became more homogenous with tree species richness. From this, we would conclude that more tree species-rich forests might have richer but more stable and homogeneous bacterial communities.
We conclude that tree diversity and leaf traits of the tree community are strong drivers of the caterpillar-associated bacteria communities in our subtropical forest. Our study revealed the linkages between tree (leaves), herbivore insects, and herbivore-associated microbes, which contributes to develop a more comprehensive understanding of relationship between herbivores and their environment. Moreover, the driving and stabilizing effects of tree diversity on herbivore-associated bacteria suggest that future research should take effects of plants on herbivore-associated microbes into consideration, when studying the relationships between plants and herbivores.

F I G U R E 5
Distance-based redundancy analysis plot showing the relationships of CWM LDMC, CWM LT, and tree richness to the bacterial community structure. The plot represents db-RDA analysis based on Bray-Curtis distance with all of the plot covariables and CWM of leaf traits as explanatory variables. CWM LDMC, CWM LT, and tree richness were three significant explanatory variables (p < .05) F I G U R E 4 Path model of the effects of tree species richness (direct effect and indirect effects through Lepidoptera abundance and Lepidoptera richness), Lepidoptera richness, CWM LT (direct effect and indirect effect through Lepidoptera richness), CWM LDMC (indirect effect through Lepidoptera richness) on richness of bacterial community. The path coefficients next to the arrows represent the strength of the positive or negative effects of one variable on another (**p < 0.001; *p < 0.05). See Table S1 and S3 for abbreviations and statistical values

ACK N OWLED G M ENTS
We thank Xiaojuan Liu, Shan Li, Bo Yang, Yu-Xi Xue, and all members of the BEF-China consortium who coordinated and helped with the establishment and maintenance of the experiment. We thank Ren-Jie Zhang and several local assistants for their help in the sampling.
We are grateful to Wenzel Kröber for the assessment of tree traits.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Caterpillar-associated bacterial sequence data can be downloaded from NCBI, accession Number: PRJNA698461. The data of the plot information that support the findings of this study are available from the corresponding author upon reasonable request.