Latitudinal pattern in community‐wide herbivory does not match the pattern in herbivory averaged across common plant species

The latitudinal herbivory hypothesis (LHH) predicts that plant losses to herbivores decrease from low to high latitudes. Although the LHH is a community‐level hypothesis, it has been rarely tested with data on community‐wide herbivory, i.e. the percentage of annual production of foliar biomass consumed by insects from all plant species at a given site. Therefore, we asked whether community‐wide leaf herbivory follows the same latitudinal pattern as observed for an unweighted average of herbivory across common plant species. We selected 10 study sites in boreal forests from 60 to 69°N along a 1,000‐km long latitudinal gradient in NW Russia. We measured relative foliar losses to insect herbivores in seven woody plant species (jointly comprising over 95% of the community‐wide above‐ground biomass) and estimated their current‐year foliar biomass. We averaged leaf herbivory for all seven species and calculated community‐wide leaf herbivory by weighting the relative foliar losses of each plant species against the contribution of that species to the annual foliar biomass production. Leaf herbivory was five‐fold higher in deciduous species than in conifers. Latitudinal patterns in herbivory varied from a significant poleward decrease in all deciduous species to a significant poleward increase in Norway spruce. Herbivory values, averaged across seven plant species, decreased with latitude and followed the pattern observed in deciduous plants due to their higher foliar losses compared with conifers. By contrast, community‐wide herbivory did not change with latitude. This discrepancy emerged because the proportion of deciduous plant foliage in the community increased with increasing latitude, and this increase counterbalanced the simultaneous poleward decrease in losses of these species to insects. Synthesis. The herbivory measured by averaging relative losses of individual plant species and community‐wide herbivory is likely to show different latitudinal patterns in various plant communities. The contributions of plant species to the total foliar biomass production should be taken into account in studies of spatial patterns of herbivory which test community‐level hypotheses. This approach may provide new insight into macroecological research on biotic interactions and improve our understanding of the role of insect herbivores in ecosystem‐level processes.

Although the LHH was originally formulated as a community-level hypothesis (Anstett, Nunes, Baskett, & Kotanen, 2016;Lim, Fine, & Mittelbach, 2015), it has been rarely tested with data on community-wide herbivory, i.e. on the percentage of annual production of foliar biomass by all plant species consumed by insects at a given site. More generally, the studies providing accurate estimates of community-wide herbivory remain extremely rare (but see Rheubottom et al., 2019;Souza, Santos, Wirth, Leal, & Tabarelli, 2013;Suzuki, Kitayama, Aiba, Takyu, & Kikuzawa, 2013;Zhang, Adams, & Zhao, 2011). Consequently, they comprise only a very low fraction of the original data used for generalizations in the published reviews and meta-analyses (e.g. Anstett et al., 2016;Moles et al., 2011;Schemske et al., 2009). Instead, studies aimed at testing the LHH are most commonly based on data collected from one or a few plant species selected from plant communities for different reasons (e.g. species that is most common or dominating in the community and/or can be found in all sites along the gradient); and this non-random selection is likely to bias resulting conclusions regarding both absolute levels of plant losses to insects and spatial patterns in herbivory (Zvereva & Kozlov, 2019). Furthermore, species-level measurements of insect herbivory cannot be directly translated into a community-wide estimate of herbivory because foliar biomass is unevenly distributed among plant species in the community.
Latitudinal studies of insect herbivory apply two major methodological approaches. The first approach involves measurements of herbivory in a single or several (usually the most common) plant species along a latitudinal gradient. These studies (reviewed by Anstett et al., 2016;Moles et al., 2011) were usually restricted to a single biome, because only a few plant species grow naturally in several biomes. Importantly, within each study, the data on herbivory from several localities were collected simultaneously (or nearly simultaneously) using the same methods. Meta-analysis of these studies (Moles et al., 2011) neither supported the LHH nor identified any other macroecological pattern. The main reason for this lack of pattern lies in the great variability between the primary studies, as some supported the LHH (Garibaldi, Kitzberger, & Ruggiero, 2011;Kozlov, 2008), whereas others found no latitudinal trends (Adams, Rehill, Zhang, & Gower, 2009;Andrew & Hughes, 2005) or even demonstrated an opposite pattern to one predicted by the LHH del-Val & Armesto, 2010). The variation in latitudinal trends in herbivory among individual study systems was claimed as one of the most important sources of controversy in studies testing the LHH (Anstett et al., 2016).
The second approach used in latitudinal studies derives latitudinal trends in herbivory from global databases, which include data on insect herbivory that have been collected in multiple localities and from multiple plant species (Kozlov, Lanta, et al., 2015;Lim et al., 2015;Zhang et al., 2016). The larger part of these data was extracted from studies that were not designed to explore latitudinal patterns in herbivory; therefore, the databases contain information from all biomes.
The analyses of these data also yielded variable conclusions regarding latitudinal patterns in herbivory. For example, Zhang et al. (2016) found a poleward decrease in herbivory in the Northern hemisphere, but not in the Southern hemisphere, whereas Kozlov, Lanta, et al. (2015) revealed a hump-shaped relationship between herbivory and latitude. Lim et al. (2015) found a decrease in herbivory with latitude for evergreen, but not for deciduous, tree species.
These two approaches both implicitly assume that data on the relative losses of foliage to insects that are collected from one plant species or averaged across a few species can be used as proxies of community-wide foliar losses in studies addressing spatial patterns of herbivory. To test this assumption, three questions should be answered: (a) whether the common plant species demonstrate similar latitudinal trends in their relative foliar losses to insects (leaf herbivory hereafter); (b) whether an unweighted average of the leaf herbivory measured in common plant species is similar to the communitywide herbivory estimates that account for the contribution of foliar biomass production should be taken into account in studies of spatial patterns of herbivory which test community-level hypotheses. This approach may provide new insight into macroecological research on biotic interactions and improve our understanding of the role of insect herbivores in ecosystem-level processes.
individual plant species to the annual production of foliar biomass and (c) whether an unweighted average of herbivory measured in common plant species and community-wide estimate of herbivory demonstrate similar changes within the same latitudinal gradient.

| Study region and study sites
The study was performed in the boreal (taiga) forest zone, from the northern limit of broadleaved forests to the northern tree line. The climate of our study region is highly seasonal. Average daily temperatures drop below +5°C in the second half of September, when deciduous trees shed their leaves. No activities of herbivorous insects occur from September to May, when new leaves start to flush.
The data were collected from 10 unevenly aged, unmanaged old-growth forests located from 60°N near St. Petersburg to 69°N close to Murmansk in NW Russia (Table S1). The sites were selected as being those closest to the rounded degrees of latitude along the road that links the above-mentioned cities. The mean annual temperature within this gradient decreases from +4.28°C in the southernmost site to −0.94°C in the northernmost site.
All sites had similar forest type; the regional tree species pool (

| Measurements of forest community characteristics
At each site, all samples were collected within a plot 100 m × 100 m in size. In these plots, we measured the basal area of the tree stands (by tree species) using a relascope at five points selected at 20-m intervals along a predefined straight line. We used an increment borer to measure the age of the five largest coniferous trees in each of our plots, and a TruPulse 200 Laser Rangefinder (Laser Technology Inc.) to measure their heights. The longevity of the needles in the conifers was quantified as the number of needle age classes simultaneously occurring on a branch. These assessments were performed for 10 mature individuals of both Scots pine and Norway spruce on a first-found, first-assessed basis and averaged to obtain site-specific values ( Table S3). The cover of the field layer vegetation (by plant species) was evaluated visually in ten 1 m × 1 m plots selected at 5-m intervals along a predefined straight line. All these measurements were performed in August 2011 using standard protocols (Kozlov, Zvereva, & Zverev, 2009;West, 2003).

| Calculation of foliar biomass of forest community
The total foliar biomass of the trees in our study plots (d.w.), reflecting the amount of food available to defoliating insect herbivores, was calculated from the following allometric model: where F is the dry mass of foliage (×1,000 kg/ha), A is the stand age (years), H is the average height of a stand (m) and B is the basal stand area (m 2 /ha). For calculation of the regression coefficients, the model was transformed to a linear form involving log-corrections (Baskerville, 1972). The regression coefficients were obtained by analysing data from 615 stands in Northern and Central Europe (for summary information, consult Table S4), which were extracted from the database published by Usoltsev (2020). We used data from all stands located between 50°N and 70°N from Central Europe to Western Siberia. The vast majority of primary studies included in this database did not distinguish between Betula pubescens and B. pendula; therefore, regression analysis combined these two species.
Shortage of data on Salix spp. did not allow for regression analysis; therefore, we used average data from four available plots as a rough estimate of the foliar biomass of willows.
Importantly, addition of the latitude to our regression models did not improve their performance. For all study species, the regression coefficients for latitude did not reach the conventional significance level (p = 0.05), and Akaike criterion was smaller for models that did not include latitude (data not shown). Therefore we concluded that our models could be applied to our data without correction for latitudes of study sites.
The obtained species-specific models (Table S5) explained, on average, 67.5% of the total variation. This value is typical for models predicting foliar biomass from tree/stand measurements (e.g. Bi et al., 2010). The species-specific models were fitted by the plot-specific values of A, H and B, and the obtained estimates of foliar biomass were multiplied by the proportions of each tree species in the basal stand area in respective plots (Table S1). For two evergreen (coniferous) tree species, the current-year needle biomass (F c ) was calculated from the total needle biomass produced by the allometric models (F) using the plot-specific number of needle age classes (NAC, full years; Table S2), as follows: F c = F × NAC/ (1 + 2 + ⋯ + NAC). For deciduous tree species, F c = F.
The foliar biomass of bilberry was calculated from its cover (Table S2) using the regression equations published by Lehtonen et al. (2016). The obtained values of the above-ground biomass were multiplied by the coefficient 0.236, which was obtained by separating the leaves from the woody parts in five samples of bilberry from each of three study sites (located at 60, 64 and 69°N), drying the samples at +105°C for 48 hr and weighing them to 0.1 mg.
The annual production of foliar biomass at each site was obtained by summing the species-specific current-year foliar biomass across all seven plant species (listed in Section 2.1) that jointly comprised over 95% of the foliar biomass of the forest community across our study sites. We accounted for all trees growing in our study sites and for about 2/3 of the field layer vegetation. The latter contributes only 1%-2% to the community-wide biomass in forests (Bolte, Lambertz, Steinmeyer, Kallweit, & Meesenburg, 2004;Gilliam, 2007); therefore, the obtained values for foliar production are sufficiently accurate.

| Assessment of herbivory
We assessed leaf herbivory in the same woody plant species as used to evaluate the current-year foliar biomass (listed above). Both our observations and published data on plant-feeding organisms of Northern Europe indicate that all types of damage recorded in the course of our study were imposed by insects.
On 17-22 August 2014 (when the majority of insect herbivores had completed their feeding), at each site we haphazardly selected five (deciduous plants) to 10 (coniferous plants) mature individuals of each of our study species. Each plant was located at least 10 m apart from others of its species. From each of the deciduous individuals, we collected a haphazardly selected branch with approximately 80-120 leaves. In conifers, we collected a haphazardly selected branch with 120-200 current-year needles. To avoid the impact of unconscious biases on the values of herbivory, the tree branches and bilberry stems were selected while standing at a distance of 5-10 m away, which does not allow visual evaluation of leaf damage by insects.
Leaf herbivory was measured in the laboratory within two days from sampling. In broadleaved species, the leaves on each branch were counted, and each leaf was carefully examined for the presence of damage imposed by chewing insect herbivores (both miners and defoliators). In conifers, we examined 100 current-year needles, starting from the tip of the branch, for traces of insect feeding.
During this examination, we also counted needles that were missing from the shoot. We attributed the loss of entire needles on currentyear shoots to herbivory, because needle damage triggers its premature abscission (Zvereva & Kozlov, 2014). Mechanical damage of needles rarely occurs in our study region; it generally impacts multiple needles on a certain part of a shoot and therefore can be easily distinguished from insect herbivory.
Following a widely used methodology (Alliende, 1989;Kozlov, Lanta, et al., 2015), each leaf/needle (leaf hereafter) was assigned to one of the damage classes according to the percentage of the area of the leaf that was consumed or otherwise damaged by insects: 0 (intact leaves), 0.01%-1%, 1%-5%, 5%-25%, 25%-50%, 50%-75% and 75%-100%. Petioles of fully consumed leaves were also recorded and attributed to 100% damage. The plant-specific percentage of foliage lost to insects (i.e. leaf herbivory level) was calculated as follows: the numbers of leaves in each damage class were multiplied by the respective median values of the damaged leaf area (i.e. 0 for intact leaves, 0.5% for the damage class 0.01%-1%, 3% for the damage class 1%-5%, etc.); the obtained values were summed for all damage classes and divided by the total number of leaves (including undamaged ones) in a sample.

| Data analysis
For species-specific relative leaf herbivory we used ln(1 + √x) transformation, which allows the inclusion of the zero values in our analysis and yields a distribution of residuals that is close to a normal distribution. Sources of variation in species-specific herbivory were explored with ANCOVA (SAS GLIMMIX procedure, type III tests: SAS Institute, 2009). We considered the tree species as fixed effect, latitude as a covariate, and the study site as a random intercept effect. The simultaneous involvement of both site and latitude in our analysis is justified by the fact that our sites differ not only in latitude, but also in a number of other characteristics (plant community structure, in particular). To facilitate accurate F tests of the fixed effects, we adjusted the standard errors and denominator degrees of freedom by the latest version of the method described by Kenward and Roger (2009). The significance of a random factor was evaluated by calculating the likelihood ratio and testing it against the chi-squared distribution (as described in Littell, Milliken, Stroup, Wolfinger, & Schabenberger, 2006). The deviations of the regression slopes from zero were explored by t test.
The average leaf herbivory levels were calculated as unweighted means of untransformed species-specific relative foliar losses to insects across all seven study species. The community-wide leaf herbivory was calculated as follows: the foliar biomass of each of seven woody plant species was multiplied by the respective untransformed percentage of the species-specific foliar losses to insects; the obtained values were then summed across all plant species and divided by the total annual production of foliar biomass in a site (for more details, consult the description of community-weighted estimates of biomass loss to insects in Rheubottom et al., 2019). The site-specific average values of leaf herbivory were compared with the respective values of community-wide leaf herbivory across all study sites by a paired t test. The latitudinal patterns in foliar biomass production, community-wide leaf herbivory and leaf herbivory averaged across common plant species were explored by calculating Pearson product-moment correlation coefficients.

| Latitudinal changes in leaf herbivory in individual plant species
Leaf herbivory was significantly associated with the latitude of the study site (F 1,9.69 = 26.2, p = 0.0005), but the latitudinal trends varied among plant species (latitude × species interaction: F 6,402.6 = 13.1, p < 0.0001; Figure 1). All deciduous woody plants demonstrated a poleward decrease in herbivory (Figure 1a-d,g), and the rate of this decline showed marginally significant variation among species (latitude × species interaction: F 4,211.5 = 2.26, p = 0.06). The two coniferous species showed different latitudinal trends (F 1,183.5 = 7.25, p = 0.0078): herbivory on Scots pine did not correlate with latitude (Figure 1f), while herbivory on Norway spruce demonstrated a significant poleward increase (Figure 1e).
Overall, herbivory levels in deciduous and coniferous plants showed opposite latitudinal patterns (latitude × leaf longevity interaction: F 1,408.3 = 61.4, p < 0.0001; Figure 2). The leaf herbivory, averaged across all seven plant species, significantly decreased with an increase in latitude (Figure 3a).

F I G U R E 1
Latitudinal changes of foliar losses to insects (leaf herbivory) in seven common woody plant species. Values (transformed as ln(0.1 + √x)) are the estimated marginal means from the SAS GLIMMIX procedure; t tests refer to the deviations of the regression slopes from zero. Note the variation in the scales of the vertical axes

| Latitudinal changes in foliar biomass
Current-year foliar biomass (per area unit; Table 1) decreased with an increase in latitude (r = −0.84, n = 10 sites, p = 0.0023), mostly due to the steep poleward decrease in foliar biomass of conifers (r = −0.83, n = 10 sites, p = 0.0027), while the foliar biomass of deciduous woody plants did not correlate with latitude (r = 0.18, n = 10, p = 0.61). As a result, the proportion of deciduous foliar biomass in the total foliar biomass of the community increased towards the north (Figure 4).

| Community-wide herbivory versus herbivory averaged across common species
The community-wide foliar losses to insects, averaged across 10 sites, were 1.90 ± 0.39%, and their site-specific values differed significantly from the respective levels of herbivory averaged across study species (paired test: t 9 = 2.44, p = 0.04). Community-wide foliar losses did not correlate with either latitude (Figure 3b) or leaf herbivory averaged across study species (r = 0.35, n = 10 sites, p = 0.31).

| D ISCUSS I ON
Our study was conducted within one biome-a boreal coniferous forest. This choice allowed us to compare insect herbivory among sites with nearly the same plant species composition and, consequently, to avoid problems related with differences in plant community composition among biomes.
The leaf herbivory measured in our sites demonstrated nearly 10-fold differences between the most-damaged and least-damaged plant species. The levels of variation in these values are consistent with those reported earlier from roughly the same study area: from below 1% to 20% (Galasjeva & Pisareva, 1991;Kozlov, 2008;Kozlov, Filippov, Zubrij, & Zverev, 2015;Kozlov, Lanta, et al., 2015;Larsson & Tenow, 1980). Thus, both the measured herbivory levels in our sites and the magnitude of their variation among species are typical for background insect herbivory in boreal forests.
We found a significant variation in the latitudinal changes of leaf herbivory among the seven common species of woody plants.
This result is in line with several earlier studies, which reported different patterns in herbivory within the same latitudinal gradient in co-occurring plant species (Kozlov, 2008;Kozlov, Skoracka, Zverev, Lewandowski, & Zvereva, 2016;Kozlov, Stekolshchikov, et al., 2015). This variation may be at least partly explained by differential responses of plant-herbivore systems to ambient temperature, which is an important factor driving latitudinal trends in many ecosystem processes (De Frenne et al., 2013), herbivory in particular (Bale et al., 2002;Kozlov et al., 2013;Kozlov, Stekolshchikov, et al., 2015).
For example, the relationships between the loads of sap-feeding herbivores and the mid-summer temperatures along latitudinal gradients differ among Scots pine, Norway spruce, downy birch and F I G U R E 2 Latitudinal changes of foliar losses to insects (leaf herbivory) in woody plant species groups that differ in leaf longevity. (a) Deciduous species (Betula pubescens, B. pendula, Populus tremula, Salix caprea, Vaccinium myrtillus); (b) evergreen species (Pinus sylvestris, Picea abies). Values (transformed as ln(0.1 + √x)) are estimated marginal means from SAS GLIMMIX procedure; t tests refer to the deviations of the regression slopes from zero. Note the variation in scales of vertical axes F I G U R E 3 Latitudinal changes of foliar losses to insects (leaf herbivory) calculated by different methods. (a) Losses averaged across the seven common woody plant species; (b) communitywide losses calculated by dividing the site-specific losses of foliar biomass by the total site-specific foliar biomass production. r, Pearson correlation coefficient silver birch (Kozlov, Stekolshchikov, et al., 2015), i.e. tree species considered in our current study.
The larger part of the among-species variation was due to the great differences between deciduous and coniferous woody plants, whose foliar losses to insects showed opposite latitudinal patterns.
Meta-analysis of the experimental studies revealed that temperature elevation leads to an increase in defensive compounds in gymnosperms, but to a decrease in these compounds in angiosperm woody plants (Zvereva & Kozlov, 2006). The mid-summer (July) temperatures differ by 3.9°C between the ends of our gradient, and this difference is similar to the temperature elevation that is usually applied in the experimental studies summarized in the cited meta-analysis. Therefore, the opposite latitudinal changes in herbivory on angiosperm and gymnosperm plants observed in our gradient might be related to differential effects of temperature on the antiherbivore defences of deciduous (angiosperm) and coniferous (gymnosperm) species.
Our deciduous species experienced, on average, five-fold higher leaf herbivory than did the coniferous species. This result is in line with the conclusions by Turcotte, Davies, Thomsen, and Johnson (2014), who found that, at the global scale, herbivory in gymnosperms is less than one-seventh of that observed in angiosperms and explained this difference by the notable toughness of the needles and the presence of terpenoid resins in gymnosperms. In our sites, conifers contributed over a half of the annual foliar production; therefore, they have a stronger influence on community-wide herbivory than deciduous plants do. As a result, the community-wide levels of leaf herbivory amounted to only a half of the herbivory values averaged across all the study plants.
The average species-specific values of leaf herbivory will be equal to the community-wide losses only when all species within plant community either experience equal levels of herbivory, or contribute equally to the annual production of foliar biomass. These hypothetical situations evidently never occur in nature. Therefore, a discrepancy between community-wide estimates of herbivory and unweighted estimates obtained by averaging across plant species is likely to occur in other plant communities in various biomes. Thus, the earlier studies, which did not account for the proportions of individual plant species in community-wide biomass production, in all probability yielded biased estimates of community-wide foliar losses to insects and, consequently, of the overall effects of background insect herbivory on ecosystem-level processes.
The community-wide foliar losses will only follow the same latitudinal changes as observed for average species-specific values of herbivory when all plant species show similar latitudinal changes in herbivory, and when plant community composition, in terms of the proportions of foliar biomass of individual plant species, does not change with latitude. None of these conditions was met in our study, and they could hardly be met in any natural ecosystem. Therefore, the community-wide herbivory and unweighted average herbivory, which show different latitudinal changes in our study, could be expected to demonstrate different patterns also in many other latitudinal gradients both limited to one biome and spanning across several biomes.
In our gradient, the proportion of deciduous plant foliage increases towards the north (Figure 4), and this increase counterbalances the Note: Zero value means that tree species was absent from five point counts but usually present in our study plot, allowing to measure its foliar losses to insect herbivores (shown on Figure 2). TA B L E 1 Biomass (d.w., kg/ha) of current-year foliage of common plant species in study sites, calculated from characteristics of forest stands (Tables S1 and S3) and field layer vegetation (Table S2) F I G U R E 4 Latitudinal changes in the proportion of foliar biomass of deciduous plants in current-year foliar biomass of all woody plants. r, Pearson correlation coefficient simultaneous decrease in herbivory on deciduous species (Figure 2a).
The end result is a lack of significant latitudinal trends in herbivory at the community level (Figure 3b). At the same time, the average herbivory in our latitudinal gradient follows the pattern demonstrated by the deciduous species, which experience the highest levels of herbivore damage, but this average herbivory does not follow the pattern exhibited by the species comprising the largest part of the community biomass, i.e. coniferous trees.
The dramatic differences between coniferous and deciduous plants in the level of insect damage (Kozlov, Lanta, et al., 2015;Turcotte et al., 2014) is one of the reasons behind the inconsistency between both absolute values and latitudinal changes of community-wide and unweighted average values detected in our study. These differences may also be one of the main sources of biases in ecosystem-wide estimates of herbivory levels in some regions of the globe. The broadleaved tree species are better represented than needle-leaved species in studies of spatial patterns in insect herbivory .
Assessment of foliar damage in needle-bearing plants is more difficult, not only due to complications in evaluating the missing area of a needle, but also because of premature abscission of damaged needles (Zvereva & Kozlov, 2014). This explains why, despite the substantial ecological and economic importance of conifers and their dominance in some ecosystems, only 0.38% (Kozlov, Lanta, et al., 2015)  We conclude that the data on relative losses of foliage to insects, collected from one plant species or averaged across a few species, cannot serve as proxies for community-wide herbivory. This conclusion is likely to be valid for other plant communities in various biomes, because the situation necessary for correlation between these two estimates-equal levels of herbivory for all plant species in the community, and equal contribution of all plant species to the annual production of foliar biomass-evidently never occurs in nature.
This finding stresses the need for measuring community-wide herbivory at multiple sites around the globe by using adequate research methodology. This methodological advancement may provide new insight into macroecological studies of biotic interactions. A community-wide approach will also improve the accuracy of the data on absolute (per area unit) losses of plant biomass to insects. These data are crucial for assessment of the contribution of insect herbivores to ecosystem-level processes, such as carbon and nutrient cycling, at both global and regional scales.

ACK N OWLED G EM ENTS
We are grateful to A. Popova for assistance in data collection, to Y. Zhang for clarification of some details of his research methodology, to T. Klemola for statistical advice, and to A. Moles and an anonymous reviewer for inspiring comments to an earlier version of the manuscript. The study was supported by the Academy of Finland