An integrated leaf trait analysis of two Paleogene leaf floras

Performing the sampling

• 1) Dicot angiosperm leaves that are taxonomically determined or have at least sufficiently well-preserved morphological characteristics have been selected for studying TCTs and insect damage types (Dataset 1). At least 50% of the leaf laminae should be preserved.

• 2) Leaves at least 70% preserved were documented (Dataset 2), of which 400 leaves per site were then randomly selected for the graphical reconstruction of the leaf outline and the determination of quantitative leaf traits (Dataset 3).

• 3) For each sampling step, the taxonomic composition was documented (Fig. 2).

• 4) Multivariate analyses were performed on subsets of Dataset 3 having no missing values (Datasets 4 and 5). A description of the datasets is available in Table 2.

The qualitative specimen data are available in Table S2 and the quantitative in Table S3.

Determining the trait combination types

The TCT approach by Roth-Nebelsick et al. (2017) categorizes dicot leaf architectural types in an assemblage. The classification of fossil leaves is made by using a hierarchic scheme based on four morphological characteristics: (1) lobation (lobed vs. unlobed), (2) leaf margin type (entire margins vs. toothed margins), (3) primary venation type (pinnate vs. palmate primary venation), and (4) secondary venation type (looped vs. non-looped secondary venation). In total, 16 TCTs are possible, each named with a capital letter from A to P (Fig. 3). The morphological characteristics were determined following the Manual of Leaf Architecture (Ellis et al., 2009). The assignment of a TCT to a specimen was made from direct observations of the characters as much as possible (specimen-based TCT approach; Roth-Nebelsick et al., 2017). In cases of poor preservation or visibility of secondary venation characteristics for specimens of well-determinable fossil-species (e.g., Engelhardia orsbergensis, Platanus neptuni, or Rosa lignitum), TCTs were determined based on the fossil-species diagnoses (combined specimen- and taxonomy-based TCT approach; Kunzmann et al., 2019). The specimen-based, taxonomy-based, and combined TCTs are stated in Table S2. Specimens in which TCTs could not be determined directly or by taxonomy (e.g., in cases where the fossil-species diagnosis implies multiple possible TCTs) were excluded from the analysis. The chipping of the coalified material at the Seifhennersdorf specimens did not affect the TCT determinations (e.g., characters are visible in better-preserved areas of the lamina or from the impression in oblique light). Poor visibility of secondary venation types was usually related to very dark coalified leaf laminae.

Reconstructing the leaf outlines

Because the lamina of fossil leaves is often incompletely preserved, digital reconstruction techniques to replenish missing parts of the leaf lamina outline are necessary before quantitative traits can be measured. A reliable reconstruction is based on the taxonomic description of a fossil-species and its morphological variability. According to Moraweck et al. (2019), the replenishment of leaf lamina should be at most 30% (meaning at least 70% of the original leaf lamina should be preserved). The ratio between the preserved and replenished leaf lamina is given by the area index (IA) or preservation index (Traiser et al., 2018; Moraweck et al., 2019). Image acquisition, processing, and graphical reconstruction were made according to Traiser et al. (2018). The 400 selected leaves per site (Dataset 3) were photographed with a Canon Power Shot G1 X Mark III. After image processing, the minimum (preserved) and maximum (replenished) leaf lamina outlines were generated as separate layers and saved in a shapefile format using the open-source geographic information system QGIS (version 3.10.14-A Coruña and 3.12.1-Bucuresti; QGIS Development Team, 2020). The method allows storing leaf lamina reconstructions, checking reliability, and ensuring data reproducibility.

Measuring the leaf quantitative traits

The maximum leaf lamina outlines were used for measurements of leaf length (in mm), leaf width (in mm), and leaf area (in mm²). The leaf length/width ratio was calculated. LM_A (in g/m²) was calculated by using the Eq. (1) for dicot leaves by Royer et al. (2007). Because this formula requires measurements of the petiole width (in mm) at the insertion point of the lamina and that petioles are often not preserved, LM_A could only be retrieved for 175 fossil leaves from Seifhennersdorf and 143 fossil leaves from Suletice-Berand (Table S3). Means for leaf area and LM_A were calculated for a fossil-species with at least five leaves.

(1) $$\mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{L}{\mathrm{M}}_{\mathrm{A}})=3.070+0.382\times \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{p}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{w}\mathrm{i}\mathrm{d}\mathrm{t}{\mathrm{h}}^2/\mathrm{f}\mathrm{o}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{a})$$

Inferring the fossil-species phenology

Since the taxonomic affiliation of most fossil-species in the studied assemblages is known (Kvaček & Walther, 1995; Walther & Kvaček, 2007), comparisons with modern relatives were used to infer the phenology of the fossil-species and thus to distinguish between evergreen and deciduous. Additionally, for fossil-species with at least five measurable leaves (Tables S4 and S5), the mean LM_A was used as an approximation of phenology as introduced by Royer et al. (2007): LM_A < 87 g/m² = deciduous, LM_A between 87 and 129 g/m² = intermediate, and LM_A > 129 g/m² = evergreen. Phenology information derived from modern relatives was compiled and used in multivariate analysis to assess the congruence of these estimates and the LM_A data.

Determining the insect damage types and herbivory metrics

First, the leaves in Dataset 1 were scored for the presence or absence of insect herbivore traces, further called insect damage types (DTs). Each DT was classified using Labandeira et al. (2007) DT catalog. A short description and an identification number define a DT. Further, DTs are assigned to terrestrial arthropod functional feeding groups (FFGs), such as external foliage-feeding, leaf-mining, or galling. Three criteria are used to distinguish DTs from detritivory or mechanical destruction: (1) the presence of distinctive, thickened tissues (callus), (2) structural features of herbivore-induced alterations such as resistant veinal stringers, necrotic tissue flaps or areas of secondary fungal infection and (3) distinctive stereotypy of a feeding pattern on a particular fossil-species (Labandeira, 2013). As mentioned above, the Seifhennersdorf specimens were affected by chipping coalified/compressed material. The extent of chipping varies from specimen to specimen and may also have affected the herbivory results. External foliage-feeding DTs were securely determinable since they were also recognizable as imprints. A loss of galls and leaf mines due to the chipping of material cannot be ruled out, which could reduce these FFGs’ occurrences and diversities. Therefore, herbivory metrics are probably partly affected by the state of leaf preservation. However, the scale is difficult to estimate.

Second, the sample-based metrics damage frequency, DT occurrence, and DT richness were used to describe, quantify, and compare insect herbivory patterns. The damage frequency (DAM%) represents the proportion between damaged and undamaged leaves in a sample (query DTs present or absent). The DT occurrence (DTO) states a sample’s number of DTs (absolute frequency). Since a leaf can have multiple occurrences, we did not standardize this metric by the number of leaves. The DT richness (DTR_{number of leaves}) reports the number of different DTs in a sample. Sample-based rarefaction procedures (Gotelli & Colwell, 2001) were applied using the free statistical software R (R Core Team, 2020) to standardize DT richness and ensure comparability between the fossil sites. In detail, we used the extension of analytical rarefaction by Gunkel & Wappler (2015). The total, external foliage-feeding (hole-feeding, leaf margin-feeding, leaf surface-feeding, leaf-skeletonization), leaf-mining, and galling DT richness were computed. The standard deviation (SD) for the resamples was calculated. The DT matrices for calculating the figured rarefaction curves are provided as supplemental information (Data S1–S6). The three metrics were applied to the two entire leaf assemblages and fossil leaves with assigned TCTs, fossil-species phenology (deciduous, evergreen, or semi-deciduous), and taxonomic affiliation (Dataset 1). The cut-off for rarefaction was 20 leaves.

Third, for the specimens showing DTs in Dataset 3, the damaged area (in mm²) was measured and related to the fossil leaf area (in mm²) to determine the herbivory index (HI). Damaged area and HI were used as specimen-based herbivory metrics for generalized linear model analyses in addition to DTs’ presence/absence data.

Performing the statistical analyses

Several statistical uni- and multivariate analyses were made using R 4.0.3 (R Core Team, 2020) to (1) test for the significance of differences in single leaf traits between Seifhennersdorf and Suletice-Berand, to (2) describe leaf characteristics from multiple quantitative traits and highlight possible trait co-variation (principal component analysis, PCA) and (3) test for the relationship between specific leaf characteristics (e.g., quantitative leaf traits and insect herbivore traces) with assemblage and vegetation properties (several generalized linear models, GLM). The R script is available in the supplemental information (Data S7).

For univariate analysis, Shapiro-Wilk tests were used to evaluate the data normality (R Core Team, 2020). This assumption being often not met in our data, differences for specific leaf traits were tested with the non-parametric Mann-Whitney test (R Core Team, 2020). Kruskal-Wallis post-hoc tests were used to compare multiple group values (e.g., leaf size or LM_A variation among plant families and TCTs; function “kruskal”, Package agricolae; Mendiburu, 2020). Quantitative leaf traits were analyzed with a PCA (i.e., five traits: LM_A, leaf size, leaf width, leaf length, leaf length/width ratio) using the packages FactoMineR (Lê, Josse & Husson, 2008) and factoextra (Kassambara & Mundt, 2020). This method is commonly used to characterize the functional diversity of present-day vegetation and is frequently applied in paleontological research (e.g., Díaz et al., 2016; Segrestin et al., 2021). Morphospaces were computed based on the two first principal components to illustrate leaf morphological disparity among the specimens of both assemblages, taxonomic units, and TCTs (Fig. 4). Five GLMs were made to highlight the possible relationship of leaf size (M1), LM_A (M2), and herbivory interaction (presence/absence of interaction (2), the area damaged (4), and HI (5)) with ecological parameters (i.e., locality, taxonomy, phenology, growth form, and TCTs), preservation (IA), and for herbivory interactions, leaf traits. Plant-insect interactions were modeled in two steps, with a first GLM exploring the probability of showing DTs (presence or absence of interactions, binomial model) and then two GLMs investigating the relationship of the leaf surface damaged with leaf and assemblage properties, using two different metrics: the area damaged or HI. All GLMs were parametrized with a Gamma distribution and identity link, except for M3 being a binomial model (link = logit). The best-fitting models, whose formulas are given in the results, were selected based on the Akaike Information Criterion (AIC, function “step,” package Stats; R Core Team, 2020). For these, the significance of explanatory variables was tested with a type-II ANOVA (function “Anova,” package car; Fox & Weisberg, 2019), whose validity was checked by investigating normality and homogeneity of the residuals with the function “check_model” (package performance; Lüdecke et al., 2021). R-squared were obtained to assess the different model’s explanatory power; they correspond to Nagelkerke’s R2 or Tjur’s R2 (for the binomial presence/absence model, M3; function “R2”, package performance; Lüdecke et al., 2021).

The presence of missing data in a dataset significantly reduces the number of specimens used in multivariate analyses (i.e., observations for all selected explanatory variables are mandatory). Within our data, the low availability of LM_A estimates was strongly limiting (i.e., it was only possible to calculate it for the 228 leaves with their petioles preserved over the 800 graphical reconstructed leaves). Therefore, multivariate analyses were computed on two datasets (Fig. 1 and Table 2): a first one, containing all specimens having calculated LM_A values (n = 228 leaves, Dataset 4), and a second dataset, an enlarged version of the first one, for which some of the missing LM_A values could be simulated (n = 506 leaves, Dataset 5). LM_A completion was made by using multivariate imputation by chained equations (function “mice” of the package MICE; Van Buuren & Groothuis-Oudshoorn, 2011) with the method “norm” (Bayesian linear regression). This method has the advantage of not biasing the variance (contrary to the common use of species mean values; Van Buuren & Groothuis-Oudshoorn, 2011). Note that it was not possible to simulate LM_A values from regression models because the number of specimens per fossil-species was too low, and no sufficiently strong correlation (R² > 0.4) was found between LM_A and other continuous variables per taxa (fossil-species as genus being considered, not shown). Missing LM_A values were simulated for leaves belonging to a fossil-species with at least five calculated LM_A values in the respective assemblage. The leaves’ locality and size were indicated as explanatory variables in the simulation process. Each predicted LM_A value averages 10 values (10 imputation cycles). Different analyses were made to assess the imputation quality (Fig. S1).

For both datasets, leaf size or LM_A outliers (> or <1.5 times the interquartile distance) and observations from fossil-species with less than five leaves were removed to avoid bias in multivariate analyses (false positives). Similarly, rare TCTs (less than five observations) were clustered or deleted. A significance level, α, of 0.05 was used for all analyses. Because predicted LM_A values are consistent with the observed values (Article S2), presented are the results obtained from the dataset with partly simulated LM_A values (Dataset 5). That allowed the inclusion of many more specimens in the multivariate analyses.

Leaf morphological characteristics and TCTs

Seifhennersdorf leaves (n = 4,935) are mainly unlobed (88.89%) and toothed (82.53%), have pinnate primary venation (87.50%) and non-looped secondary venation (74.14%). Suletice-Berand leaves (n = 1,349) are mainly unlobed (94.51%) and toothed (58.49%), have pinnate primary venation (92.06%) and more frequently looped secondary venation (64.05%). As shown in Fig. 3, the assemblages are dominated by five TCTs (>20 leaves): A, B, E, F, and P. Differences between the assemblages exist in the TCT proportion. The Seifhennersdorf assemblage is characterized by a high frequency of TCT F leaves (64.78%), coming from diverse deciduous trees and shrubs like Ailanthus prescheri, Alnus phocaeensis, Betula alboides, Carpinus grandis, Carpinus roscheri, Carya fragiliformis, Cyclocarya sp., Ostrya atlantidis, Populus zaddachii, Salix varians, Ulmus fischeri and Zelkova zelkovifolia (Table S4). Leaves of TCTs E (8.17%), P (6.36%), A (3.55%), and B (2.84%) are less frequent. TCT E is mainly related to Platanus neptuni and Rosa lignitum, TCT P to fossil Acer species, TCT A to evergreen trees and shrubs like Laurophyllum acutimontanum and Magnolia seifhennersdorfensis, and TCT B to diverse legume morphotypes and Cornus sp. Twenty fossil-species show leaves of TCT F, ten of TCT E, nine of TCT A, six of TCT P, and a minimum of two fossil-species shows TCT B (e.g., it is uncertain how many fossil-species are represented by the legume morphotypes).

TCTs E (40.62%) and A (21.50%) represent most leaves with looped secondary venation in Suletice-Berand. Conversely, leaves with non-looped secondary venation are less frequent in that assemblage (TCT F: 12.08%, TCT P: 3.56%), except leaves of TCT B (11.19%) that are frequently represented by legume morphotypes (Fig. 2). Leaves of TCT E often come from Engelhardia orsbergensis, an evergreen tree, and less frequently from P. neptuni, a deciduous tree. Leaves of TCT A are from smaller evergreen trees to shrubs like Laurophyllum acutimontanum, L. pseudoprinceps, and Oleinites maii (Table S5). Like TCT E, in Suletice-Berand, leaves of TCT F can either belong to evergreen (e.g., Sloanea olmediifolia) or deciduous fossil-species (e.g., Carpinus grandis, Carya serrifolia, or Zelkova zelkovifolia)—with evergreens predominating proportionally (Table S4). Lobed and palmately-veined leaves (TCT P) are from fossil Acer species, too. Thirteen fossil-species are included in TCT F, 11 in TCT A, eight in TCT E, four in TCT P, and a minimum of three in TCT B. Both assemblages contain TCTs that are represented by less frequently preserved fossil-species (less than 20 leaves) like Cercidiphyllum crenatum (TCT G), Tilia gigantea (TCT H), and cf. Crataegus sp. (TCT M) in Seifhennersdorf or Pungiphyllum cruciatum (TCT J) and Dombeyopsis sp. (TCT K) in Suletice-Berand (Fig. 3).

Leaf quantitative traits

Based on quantitative trait analyses (Dataset 3), leaves tend to be larger and with a lower LM_A in Seifhennersdorf than in Suletice-Berand. These differences in leaf size and LM_A are also significant in Dataset 5, used for multivariate analyses (Table 2). The calculated mean LM_A of seven fossil-species from Seifhennersdorf is below 87 g/m² (Fig. 5 and Table S4)—the threshold for deciduousness stated by Royer et al. (2007). The calculated values are consistent with inferred phenology from modern relatives. In Suletice-Berand, the calculated mean LM_A of six fossil-species falls into the transition between deciduous and evergreen (i.e., between 87 and 129 g/m², Royer et al., 2007; Fig. 5 and Table S5). These LM_A values are lower than expected from the phenology information inferred from modern relatives, which suggest Daphnogene cinnamomifolia, Engelhardia orsbergensis, Laurophyllum cf. acutimontanum, and Symplocos deichmuelleri most likely to be evergreen trees or shrubs. The opposite is observed for Platanus neptuni and Zelkova zelkovifolia, which were most probably deciduous by their modern relatives but showed mean LM_A greater than 87 g/m² (Fig. 5 and Table S5).

Multivariate analyses of leaf quantitative traits

Results of multivariate analyses (PCA and GLM) are derived from Dataset 5, which contains calculated and simulated LM_A values (Fig. 1). Figures showing the quality of the LM_A imputation are available in the supplemental information (Fig. S1). The same analyses were also conducted on the smaller Dataset 4 (without imputing LM_A missing values) and are available in Article S2.

The morphospace obtained from the PCA on quantitative leaf traits is illustrated in Fig. 4. The two first axes, namely PC1 and PC2, explain 83.9% of the variance (54.2% and 29.7%, respectively). PC1 indicates the leaves’ size, width, length, and LM_A (each parameter explaining 33.1%, 28.7%, 18.1%, and 16.4% of specimen distribution along the axis, respectively). Thereby, small-sized leaves with higher LM_A are associated with negative values. In contrast, large-sized leaves with lower LM_A but higher widths and lengths are associated with positive values. PC2 mainly reflects circularity (l/w ratio explains 57.6% of point distribution along the axis), with leaves of high (resp. low) l/w ratio associated with positive (resp. negative) values.

To better understand leaf area and LM_A variabilities among leaves, GLM analyses were run. Among the variables tested, those that best explain leaf size and LM_A variations are contained in the following models:

$$\mathrm{M}1:\mathrm{L}{\mathrm{M}}_{\mathrm{A}}\stackrel{~}{}\mathrm{S}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{s}+\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}+\mathrm{p}\mathrm{h}\mathrm{e}\mathrm{n}\mathrm{o}\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{y}+\mathrm{I}\mathrm{A}(+\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}),$$

$\mathrm{M}2:\mathrm{L}\mathrm{e}\mathrm{a}\mathrm{f}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\stackrel{~}{}\mathrm{S}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{s}+\mathrm{I}\mathrm{A}+\mathrm{p}\mathrm{h}\mathrm{e}\mathrm{n}\mathrm{o}\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{y}(+\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}).$

The parameters indicated after “~” are the significant explanatory variables as given by the function “step” (see Material & Methods). The variables which best explain LM_A variations (model M1) are the taxonomic affiliation (fossil-species), leaf size, and phenology (Table S6). LM_A is negatively correlated to leaf size but positively associated with evergreen and, to a lesser extent, semi-deciduous species (Figs. 4A and 4C). A positive relationship with leaf preservation (IA) is also observed, leaves of higher LM_A tend to be more complete. This could be due to a higher preservation potential of these leaves than leaves with low LM_A (i.e., lower fragmentation before fossilization and during the fossil collection). Similarly, leaf size variations are also partly explained by taxonomic affiliation (fossil-species), leaf preservation (IA), and phenology (model M2). Negative relationships are observed between leaf size, leaf preservation, and phenology (i.e., the smallest leaves are the best preserved and are related to evergreen fossil-species).

The significant effect of fossil-species on leaf area and LM_A shows that the taxonomic composition of assemblages is important to consider during sampling (Fig. 2) and to study these traits at a paleo-community level. Thus, part of the differences in leaf size and LM_A between Seifhennersdorf and Suletice-Berand likely result from the different fossil-species compositions and proportions in the two assemblages (Fig. 4B). Some families with the largest specimens are more abundant in Seifhennersdorf, especially Sapindaceae and Betulaceae. Conversely, families with smaller specimens are more represented in Suletice-Berand (e.g., Fabaceae, Juglandaceae, Lauraceae, and Platanaceae; Fig. 5A). Similarly, families with specimens tending to higher LM_A (e.g., Fabaceae, Lauraceae, Platanaceae, Symblocaceae, and Ulmaceae) are more abundant in Suletice-Berand than in Seifhennersdorf (Fig. 5B). While significant differences were found when testing the correlation between size or LM_A and locality alone (Table S6), the locality was not found to be a significant parameter in models M1 and M2 (p-values of 0.14 and 0.09, respectively). Interestingly, this shows that the locality factor does not explain most of the variation in size (or LM_A) between specimens and highlights the interest in multivariate approaches to understand trait variations. The results are fairly consistent with those obtained on the smaller Dataset 4 without LM_A imputation (Fig. 1), although no relationship was found between leaf preservation (IA) and LM_A (Article S2).

General comparison

Herbivory metrics compared between the assemblages (Dataset 1) indicate that insect herbivory was more pronounced on leaves from Suletice-Berand (Fig. 6). Representative examples of DTs are shown in Figs. 7 and 8. In detail, 4.78% of 4,935 leaves from Seifhennersdorf and 12.08% of 1,349 leaves from Suletice-Berand were affected by insect damage. The DT occurrences between the assemblages are compared in Fig. 6B1. Regarding the number of sampled leaves, DT occurrence was denser in the Suletice-Berand assemblage. Moreover, the richness of DTs, standardized by sample-based rarefaction, is more significant in Suletice-Berand than in Seifhennersdorf (Fig. 6C1). On 400 sampled Seifhennersdorf leaves 2.59 FFGs (SD = 0.67) and 7.92 different DTs (SD = 1.68) are present, whereas 4.60 FFGs (SD = 0.81) and 15.27 different DTs (SD = 2.01) are present on the same number of Suletice-Berand leaves. In total, 19 different DTs of five FFGs are present in the Seifhennersdorf assemblage, and 23 different DTs of six FFGs in Suletice-Berand.

Most herbivory traces were caused by external foliage-feeding (Fig. 6). Common are hole feeding (DTs 2, 3, and 5), leaf margin feeding (DTs 12 and 14), and leaf skeletonization (DT 20). The assemblages differ regarding the richness of hole-feeding DTs (Seifhennersdorf: DTR_{400 leaves} = 4.13 different DTs (SD = 1.14), Suletice-Berand: DTR_{400 leaves} = 6.34 different DTs (SD = 1.21)) and leaf-skeletonization DTs (Seifhennersdorf: DTR_{400 leaves} = 0.31 different DTs (SD = 0.50), Suletice-Berand: DTR_{400 leaves} = 3.51 different DTs (SD = 0.96)). Comparable between the assemblages is the low frequency (DAM% < 0.5%), occurrence (DTO < 5 DTs), and richness (DTR_{400 leaves} < 0.5 different DTs) of galling and leaf-mining DTs (Figs. 7F, 7K, 8C, 8H, and 8I), and the absence of DTs caused by piercing and sucking insects.

Insect herbivory regarding the TCT of leaves

Overall, DTs are most common on the abundant TCTs (n > 20 leaves, Fig. 3) but with intra- and inter-site differences (Fig. 9). First, TCTs differ regarding the damage frequency, especially in comparison between the assemblages. In Suletice-Berand, a TCT is more frequently affected than in Seifhennersdorf due to the overall higher incidence of insect herbivory. Within the assemblages, the differences in the damage frequency of TCTs are more pronounced in Suletice-Berand than in Seifhennersdorf. Damage frequencies vary between 2.85% (TCT B) and 5.89% (TCT F) in Seifhennersdorf and between 4.64% (TCT B) and 18.75% (TCT P) in Suletice-Berand (Fig. 9). Leaf types with toothed margins (TCTs E, F, and P) are more frequently damaged than those with entire margins (TCTs A and B). The same applies to leaf types with non-looped secondary venation (TCTs F and P, except TCT B) compared to those with looped secondary venation (TCTs A and E). Second, TCTs differ in the occurrence of DTs (Fig. 9). In both assemblages, DT occurrences tend to be higher on leaves with toothed margins (TCT E, F, and P) than on entire-margined leaves (TCTs A and B). In Seifhennersdorf, most DTs occur on unlobed and toothed leaves with non-looped secondary venation—TCT F. In Suletice-Berand, the number of DT occurrences is higher on leaves with looped secondary venation, especially on TCT E, but also on TCT A. Third, TCTs differ regarding DT richness. According to sample-based rarefaction, DT richness for a given TCT tends to be higher in Suletice-Berand than in Seifhennersdorf, except for TCT B (Fig. 9). In Seifhennersdorf, the DT richness is the highest on leaves of TCT F with in total 15 different DTs of three FFGs. For comparison, there are only five different DTs of two to three FFGs on TCTs A, B, E, and P leaves. Based on sample-based rarefaction, the DT richness in Suletice-Berand is also highest on leaves of TCT F, followed by TCTs E and A. In total, 15 different DTs are present on TCTs E and F, 11 on TCT A, and five on TCTs B and P.

In both assemblages, hole-feeding DTs such as DTs 2–5 and leaf margin-feeding DTs such as DTs 12–14 frequently occur on leaves of different TCTs. Thus, they are not TCT-specific. There are only a few exceptions that indicate a relationship between the occurrence of certain DTs and TCTs. For example, TCT F leaves frequently showed traces of herbivory that were concentrated between secondary leaf veins (DT 78, Figs. 7C, 7D, and 8B) or at the junction of secondary and primary leaf veins (DT 50 and 57, Figs. 7H and 7J). Characteristic leaf skeletonization was observed on TCTs E and F in Suletice-Berand (DT 20, Figs. 8A and 8E). Rare evidence of leaf-mining was found on leaves of TCT E (Suletice-Berand, Fig. 8C) and of galling DTs on leaves of TCT P (Figs. 7F, 7K, 8H, and 8I).

Insect herbivory regarding fossil-species phenology and taxonomic affiliation

In the following, the relevance of the phenological and taxonomic background to the observed relationship between TCTs and insect herbivory is reviewed. The analyzed assemblages show distinct differences in the proportion of leaves originating from deciduous or evergreen fossil-species (Figs. 6A3 and 6A4). Also, herbivory metrics indicate meaningful intra- and inter-site differences regarding the relationship between insect herbivory and deciduous or evergreen leaves (Fig. 6).

The Seifhennersdorf assemblage contains a high proportion of leaves from various deciduous fossil-species, and with 76.73%, most of these leaves are assigned to TCT F (Table S4). Damage frequency, DT occurrence, and DT richness indicate an affinity of insect herbivores to fossil-species with these traits (Figs. 6 and 9). However, herbivory metrics also show a taxonomy-related influence on the observed relationship. Damage frequency varies between 2.61% (Ulmus fischeri) and 6.81% (Carya fragiliformis), DT occurrences between one (Salix varians) and 99 DTs (C. fragiliformis), and DT richness between one (S. varians) and 12 different DTs (C. fragiliformis). Fossil-species of Juglandaceae and Betulaceae differ remarkably from others by higher damage frequencies, DT occurrence, and DT richness (Table S4). Therefore, the DT record on leaves of Carya fragiliformis, Carpinus grandis, and Carpinus roscheri (Figs. 7A–7E and 7G–7J) is mainly responsible for the observed relationship between deciduousness, TCT F, and insect herbivory in Seifhennersdorf.

Unlike Seifhennersdorf, the Suletice-Berand assemblage is characterized by a higher proportion of leaves originating from evergreen fossil-species (Table S5). According to their frequencies, most of these leaves are of TCTs E (47.55%), A (27.04%), and F (7.17%). Leaves of deciduous fossil-species are mainly of TCTs E (33.23%), F (29%), and P (14.50%). Leaves of evergreen and deciduous fossil-species are more frequently damaged in Suletice-Berand than in Seifhennersdorf and show distinctly higher DT richness (Fig. 6). Within the Suletice-Berand assemblage, leaves of evergreen and deciduous fossil-species differ regarding DT occurrence and DT richness. Comparable to Seifhennersdorf, the evidence of insect herbivory in Suletice-Berand also varies between fossil-species with similar TCT or phenology—indicating the importance of taxon-specific traits for insect herbivory. The damage frequencies between the fossil-species varies between 5% (Conus studeri) and 28.57% (Acer tricuspidatum), the DT occurrences between one (C. studeri) and 49 DTs (Engelhardia orsbergensis), and the DT richness between one (C. studeri) and 11 different DTs (E. orsbergensis). Like in Seifhennersdorf, Juglandaceae leaves (E. orsbergensis) are most abundant in Suletice-Berand. Insect herbivores also frequently affected these leaves, accounting for a significant portion of the damage reported for TCT E (Fig. 8D and Table S5). TCT A leaves in Suletice-Berand also show stronger evidence of herbivory than in the Seifhennersdorf. That is related to frequently damaged evergreen Lauraceae leaves (Table S5). Leaves with non-looped secondaries (TCTs F and P) differ from leaves with looped secondaries (TCTs A and E) by showing higher damage frequencies and DT richness (Fig. 9). Again, significant portions of the DT record are related to individual fossil-species. For instance, the DTs on leaves of Sloanea artocarpites, an evergreen tree, and Zelkova zelkovifolia, a deciduous tree, significantly determined the herbivory metrics regarding TCT F (Fig. 9 and Table S5).

Multivariate analyses of insect herbivory

The different herbivory metrics could be partly explained with multivariate analysis (Dataset 5, Fig. 1). First, the presence/absence of DTs (interaction) varies significantly between the localities, among classes of TCTs (when secondary venation is not considered) and depends on the LM_A (2). It confirms the abovementioned observations: traces are often found on leaves with relatively low LM_A, from Suletice-Berand, and belonging to TCT class E/F. However, from this data (Dataset 5), no significant effect of TCTs is found when secondary venation is considered. Although included in the list of explanatory variables, leaf size does not explain the presence/absence of herbivory traces (2). Dataset 1, which does not provide quantitative trait measurements but enables to test of the effect of TCTs, phenology, fossil-species, and locality, reflects a relationship of herbivory (presence/absence) with the TCTs, especially TCT F (3). The models focusing on the presence or absence of DTs have extremely low R2, probably due to the low proportion of leaves showing DTs (i.e., ~10%, Table S6). That indicates that it is not possible to predict insect herbivory based on the quantitative traits and environmental/ecological variables used in the models.

Second, models analyzing the leaf surface damaged (area damaged and herbivory index in Dataset 5) reveal an important effect of phenology (as assessed with the modern relatives, models M4 and M5, Table S6). A positive correlation between the herbivorized leaf surface and the leaf size is also visible when the area damaged is the herbivory metric used (M4). Conversely, analyses made on the smaller Dataset 4 (without LM_A imputation, Fig. 1) did not provide any information regarding DTs and their relationship to leaf traits and/or assemblage properties. The detailed results are provided in Article S2.

(2) $$\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\stackrel{~}{}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}+\mathrm{L}{\mathrm{M}}_{\mathrm{A}}+\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{T}\mathrm{C}\mathrm{T}\mathrm{s}(+\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e})$$

(3) $$\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\stackrel{~}{}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}+\mathrm{T}\mathrm{C}\mathrm{T}\mathrm{s}$$

(4) $$\mathrm{A}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{d}\stackrel{~}{}\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}+\mathrm{p}\mathrm{h}\mathrm{e}\mathrm{n}\mathrm{o}\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{y}$$

(5) $$\mathrm{H}\mathrm{e}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{v}\mathrm{o}\mathrm{r}\mathrm{y}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}\stackrel{~}{}\mathrm{p}\mathrm{h}\mathrm{e}\mathrm{n}\mathrm{o}\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{y}(+\mathrm{L}{\mathrm{M}}_{\mathrm{A}})$$