Warming Responses of Leaf Morphology Are Highly Variable among Tropical Tree Species

: Leaf morphological traits vary along climate gradients, but it is currently unclear to what extent this results from acclimation rather than adaptation. Knowing so is important for predicting the functioning of long-lived organisms, such as trees, in a rapidly changing climate. We investigated the leaf morphological warming responses of 18 tropical tree species with early (ES) abd late (LS) successional strategies, planted at three sites along an elevation gradient from 2400 m a.s.l. (15.2 ◦ C mean temperature) to 1300 m a.s.l. (20.6 ◦ C mean temperature) in Rwanda. Leaf size expressed as leaf area (LA) and leaf mass per area (LMA) decreased, while leaf width-to-length ratio (W/L) increased with warming, but only for one third to half of the species. While LA decreased in ES species, but mostly not in LS species, changes in LMA and leaf W/L were common in both successional groups. ES species had lower LMA and higher LA and leaf W/L compared to LS species. Values of LMA and LA of juvenile trees in this study were mostly similar to corresponding data on four mature tree species in another elevation-gradient study in Rwanda, indicating that our results are applicable also to mature forest trees. We conclude that leaf morphological responses to warming differ greatly between both successional groups and individual species, with potential consequences for species competitiveness and community composition in a warmer climate.


Introduction
Structural, chemical and physiological leaf traits have a decisive influence on the productivity of plants and their ecosystems [1,2].Leaf morphological traits, such as leaf mass per unit area (LMA), leaf size and shape, influence the interception of light and the exchange of gases and energy between plants and the atmosphere, thereby controlling plant photosynthesis, transpiration and leaf temperature [3][4][5][6].While it is well known that leaf traits are adapted to the environment where the species occur, the possibility to acclimate these traits to altered environmental conditions is crucial to the success with which plants and their ecosystems will respond to climate change [7][8][9].
Forests 2022, 13, 219 3 of 24 acquisition, carbon gain and fast growth, while LS species typically grow slower and use available resources more efficiently under unfavourable conditions (e.g., drought, shade) [52][53][54].Previous experiments with tropical trees have indicated that ES species respond more favourably (or less negatively) to warming than LS species.This has been shown for plant growth [55] as well as for photosynthesis [56].The plasticity of leaf traits, including LMA, in changes in light environment is also thought to be greater in ES compared to LS species [57].Whether or not leaf morphological traits also show larger acclimation capacity to warming in tropical ES compared to LS species remains unknown.Furthermore, leaf morphology may also change with ontogenetic development, i.e., between juvenile and mature trees.LMA is normally higher in mature compared to juvenile trees [21], while leaf size is more ontogenetically variable [58].Whether there are differences in leaf morphological responses to warming between young and mature trees remains to be investigated.
The East African region, with high mountains and great lakes, located between the monsoon domain of West Africa and the Indian ocean [59], faces some of the largest interannual rainfall variations in the world [60].It is projected that precipitation extremes will increase between 2030 and 2052 in Sub-Saharan Africa and that East Africa will experience increases in precipitation in some regions and decreases in others [61].Warming projections for the region exceed the global projection of 2 • C above pre-industrial levels by the end of the century [61,62].In the future, these shifts in seasonal variations in rainfall and temperature are likely to affect tree species composition, phenological patterns (i.e., deciduous vs. evergreen) and carbon balance of tropical forests [63].
Afromontane tropical forests are highly productive and store large amounts of carbon compared to tropical montane forests (TMF) worldwide and to central-east Amazonian tropical lowland forests [64][65][66].However, they are severely understudied and very little is known on how leaf morphological traits are influenced by climate and how they vary among species and successional strategies.To address this knowledge gap, we explored how leaf morphological traits respond to a warmer climate in a broad range of tropical tree species native to central and east Africa.The main study was conducted on juvenile trees in three multi-species common garden plantations along an elevation gradient.The results from these trees were compared to results from mature trees of the same species growing at different elevations in another elevation gradient.The elevation gradients were used to simulate different global warming scenarios, and the following hypotheses were tested for well-watered trees: (i) LMA decreases at warmer sites, with larger contributions from shifts in LD than from shifts in LVA (leaf thickness); (ii) Leaf size increases at warmer sites; (iii) LMA and leaf size of early successional species are more responsive to warming compared to late successional species; (iv) LMA and leaf size differ between juvenile and mature trees, but the responses to warming are similar.

Experimental Sites
This study was conducted during 2018 and 2019 in the TRopical Elevation Experiment in Rwanda (Rwanda TREE; see www.rwandatree.com(9 December 2021)) where 20 tropical tree species native to central and east Africa were grown at three experimental sites located at different elevations to expose the trees to different climates.A step down in the gradient was meant to represent a possible future climate-warming scenario.The sites were: (i) Sigira, high-elevation (HE) site, located in Nyamagabe district at 2 • 30 54 S, 29 • 23 44 E at the edges of Nyungwe national park in the tropical montane forest (TMF) vegetation zone at 2400 m a.s.l.; (ii) Rubona, mid-elevation (ME) site located in Huye district at 2 • 28 30.2 S, 29 • 46 49.0 E in the Lake Victoria transitional rainforest (LVTF) vegetation zone at 1600 m a.s.l.; (iii) Makera, low-elevation (LE) site located at Ibanda Makera in Kirehe district at 2 • 6 31 S, 30 • 51 16 E in the evergreen and semi-evergreen bushland and thicket vegetation zone at 1300 m a.s.l.(Table 1).
Table 1.Site, weather and soil characteristics at Rwanda TREE experimental sites.Weather data are annual mean ± SD for the period 1 February 2018-31 January 2020, except for wind, which is given for 1 February 2019-31 January 2020.Soil data are mean ± SD of three blocks (1 block = 6 subplots).HE, high elevation; ME, mid elevation; LE, low elevation; MAT, mean annual temperature; MAP, mean annual precipitation; T air day and night, mean air temperature in light and darkness (< and >2 µmol m −2 s −1 PPFD), respectively; VPD day, daytime mean vapour pressure deficit; PPFD day, daytime mean photosynthetic photon flux density; gust wind speed, average of maximum half-hourly wind speeds; T soil, mean soil temperature at 10-20 cm depth; SWC, soil water content at 10-20 cm depth; SBD, soil bulk density; P, phosphorus, N, nitrogen; Org C, organic carbon; All soil parameters from SBD and downwards are average between 0-30 cm depth.

Planting Design and Plant Material
At each site, 18 plots of 15 × 15 m, spaced by 2.5 m paths, were established on a 50 × 102.5 m area where all vegetation was cleared before planting.The 18 plots allowed for a full factorial experimental design, with three water levels and two fertility levels and a replication of three plots for each of the six treatment combinations.However, the treatments started after this study was conducted.Within each plot, 20 different tree species with a replication of 5 (i.e., 100 trees per plot) were planted using 1.5 × 1.5 m spacing.The positions of trees and species inside plots were randomized.
In this study, 18 tree species were included, representing early (ES) and late (LS) successional strategies, as well as TMF and LVTF vegetation zone origins as follows: 5 ES/TMF, 5 LS/TMF, 5 ES/LVTF and 3 LS/LVTF species (Table 2).Note that, from literature, the classification of two of the species into distinct successional groups is ambiguous (Table 2), but according to our assessment Ficus thonningii Blume and Markhamia lutea (Benth.)K.Schum Forests 2022, 13, 219 5 of 24 belong mostly to LS and ES groups, respectively, and therefore we will use that classification in this study.
Table 2. Taxonomy for species and their main forest type of origin, classification into successional group as well as their leaf type and shape.FT, forest type of origin; TMF, tropical montane forest (~>2000 m a.s.l.); LVTF, Lake Victoria Transitional Forest (~1500-2000 m a.s.l.); SG, successional group (ES, early LS, late); K, leaf area factor; *, species with mixed successional group features.

Code
Scientific Name and Author  4 Semideciduous species drop variable amounts of leaf depending the severity of drought, but are rarely completely defoliated; 5 Leaf type and shape classifications follow Ellis et al. [68], see also Figure S3.
The trees were propagated from seeds, cuttings or wildlings, depending on speciesspecific difficulties of propagation, in poly-pots in a nursery at Rubona research station during 2017.The germplasms were collected from Nyungwe TMF (high-elevation), or Rubona research station located in the LVTF vegetation zone (mid-elevation), depending on main species distribution range.After six to twelve months in the nursery and having reached a height of ~10-75 cm, depending on species, each plant was randomly assigned to site, plot and plot position (as constrained by the experimental design) and transplanted to the experimental plots within a period of approximately one month at the turn of the year from 2017 to 2018.All tree species received water when needed, irrespective of the planned water treatment, especially during the first dry periods in 2018.

Environmental Conditions
Weather stations were installed at all three Rwanda TREE sites to record ambient air and soil temperature, relative humidity, precipitation, solar radiation, soil water content, and wind speed and direction at a frequency of 30 min.The recording of most parameters started late January 2018, while wind and soil measurements started late August 2018 (Table 1).Soil temperature and water content sensors were installed in the center of 6 selected plots at each site, while other parameters were measured in an open area next to the plantations.The mean air temperatures (MAT), measured at 1.8 m above ground during the study period at the HE-site (Sigira), ME-site (Rubona) and LE-site (Makera) were 15.2, 20.0 and 20.6 • C, respectively.The daytime and extreme temperatures (expressed as 99%ile) were 17.1/23.1 • C at HE-site, 22.4/28.4• C at ME-site and 24.0/31.2• C at LE-site.The main reason for the larger difference in the daytime and extreme temperatures compared to MAT between the ME and LE sites was the larger difference between daytime and nighttime temperatures at the LE site.The sites also differed substantially in mean annual precipitation (MAP), decreasing progressively from HE (Sigira, c. 2100 mm) to ME (Rubona, c. 1700 mm) and LE (Makera c. 1100 mm; Table 1) sites.However, the relative seasonal distribution of precipitation was similar at all sites, with the highest rainfall in March-May and a dry period in June-August.Solar radiation was similar at ME and LE sites while the HE-site received less radiation, probably due to higher cloudiness.Soil temperatures were closely related to the site MAT, although it was probably also affected by radiation and canopy cover.Soil water content (SWC) was similar at LE and ME sites in spite of different MAP.The mean SWC at HE was substantially higher compared to the two other sites, probably because of both higher MAP and higher water-holding capacity due to higher soil clay content (Table 1).
In November 2017, previous to the planting of the trees, soil samples were collected at 0-10 cm and 20-30 cm depth in each plot.Soil bulk density (SBD), NH 4 + , NO 3 − and available P were analysed from composite samples of paired plots (→ 9 samples per site), while the remaining soil parameters were analysed from composite samples of six plots (→ 3 samples per site).All values presented here are averages for 0-30 cm soil depth (Table 1).The main soil texture differences between sites were a larger proportion of sand and less silt at ME site (Rubona, 53-62% and 5-9%, respectively) compared to the soil at HE (Sigira) and LE (Makera) sites (35-45% and 15-27%, respectively) while the clay content was relatively high at all sites (30-50%).The soil pH (water) was 4.2 at the HE-site and increased with approximately one unit for each step down the elevation gradient.Soil fertility expressed as total N and P content tended to decrease with decreasing elevation.

Leaf Sampling and Morphological Traits Measurements
Leaves or leaflets were collected in two campaigns (August to December 2018 and June to August 2019).For species with compound leaves (Table 2), only leaflets were sampled but for simplicity, we will denote them as leaves.In the first campaign, two randomly selected trees per species and plot were sampled, while one tree per plot and species was generally selected in the second campaign.However, for some species and plots, one additional tree was sampled during the second campaign, but sampling was always balanced between sites.Only plot averages were used for statistical analysis of site effects.In each campaign, one or more (2-3 from species with small leaves) mature and sun-exposed leaves from the upper half of the tree crown were sampled from each selected tree, resulting in: 1 sample * 2 trees * 18 plots * 18 species * 3 sites = 1944 leaves in 2018 and 1 sample * 1 tree * 18 plots * 18 species * 3 sites + 190 extra samples = 1162 leaves in 2019, giving a total number of 3106 leaf samples.
Leaf length (LL), leaf width (LW) and leaf thickness (LT, only in 2019 campaign) were measured directly after leaf collection, using a ruler and an electronic calliper (resolution 0.01 mm).Three to five leaf discs of known diameter (18 mm or 10 mm, depending on the size of the leaf) per leaf and tree were sampled using punchers.First-and second-order veins were avoided for LT measurements, as well as first-order veins for sampling of discs.For species with narrow leaves (Afrocarpus falcatus (Thunb.)C.N.Page and Faurea saligna Harv.) discs were not taken; instead, a photo was taken of the whole leaf next to a ruler for subsequent leaf-area determination.Both discs and photographed leaves and the remaining leaf material were brought in separate envelopes to the laboratory for later determination of LMA and nutrient contents, respectively.

Leaf Shape and Size Estimation
Leaf shape was determined as LW to LL ratio (leaf W/L ratio) and leaf size was expressed as leaf area (LA) and estimated from LW and LL using an allometric function: where K is the leaf-area factor, also known as the Montgomery parameter [69].Speciesspecific K-values were determined for the 18 species based on 20 leaves of varied sizes collected from each species at the ME site (Rubona) in June 2020.Collected leaves were scanned in a flatbed scanner (CanonScan LiDE300, Canon Inc., Tokyo, Japan) and thereafter the LA, LL and LW were analysed using ImageJ software 1.50i (Rasband, W.S., U. S. National Institutes of Health, Bethesda, MD, USA).Petioles and petiolules were excluded from the measurements.The K-values were obtained by regression analysis of LA versus LW*LL, setting the intercept to zero.Independent of fixed or variable intercept, the R 2 values obtained were 0.97-0.99,except for the two species with the largest and lowest LW to LL ratio, where the R 2 values were 0.93 (Dombeya torrida (J.F.Gmel.)Bamps) and 0.89 (A.falcatus), respectively.The K-values among species varied between 0.66 and 0.85 (Table 2) and were used to estimate the LA of all sampled leaves.

Determination of LMA, LD and LVA
Each collected leaf disc or photographed leaf sample was oven-dried at 70 • C for at least 48 h and then weighed using a laboratory balance with 0.1 mg resolution.The area, width and length of the photographed leaves were determined using ImageJ software 1.50i.Total projected areas of discs or photographed leaves were used to calculate LMA [70].Leaf area (LA) and leaf thickness (LT) were used to calculate leaf volume per area (LVA), while leaf mass (LM), LA and LT were used to calculate leaf density (LD) as follows: The following units are used throughout the paper: Leaf size (LA), cm −2 , LMA, g m −2 , LVA, cm −3 m −2 , LD g cm −3 ).

Leaf Nutrients Analysis
Composite samples of leaves (from six plots into one sample per species and site, constituting one block) were oven-dried at 70 • C for at least 48 h, and thereafter ground into a fine powder with a ball mill grinder (MM 301, MM 200, Retsch, Germany).Mass-based leaf nitrogen (N M ) content was analysed using an elemental analyser (EA 1108, Fison Instruments, Rodano, Italy).Ground samples were sent for analysis of 37 non-N elements using inductively coupled plasma mass spectrometry after digestion in HNO 3 and then aqua regia (Method VG101, Bureau Veritas Mineral laboratories, Vancouver, BC, Canada).Out of these non-N elements, only data for mass-based leaf phosphorus (P M ) are reported.The mass-based leaf nutrients were converted to area-based contents (N A and P A ) by multiplying with LMA.

Comparison of Juvenile Planted Trees with Mature and Natural Regenerated Trees
The leaf size and LMA of two ES species (Macaranga Kilimandscharica Pax and Polyscias fulva (Hiern) Harms) and two LS species (Carapa grandiflora Sprague and Syzygium guineense (Willd.)DC) in this study were compared to values obtained from two other studies on mature trees in south-west Rwanda: (i) an elevation gradient with five sites ranging 1700-2700 m a.s.l.[71] and (ii) 11 permanent monitoring plots in Nyungwe TMF [66].
In the mature-tree elevation-gradient study, leaves from trees with an average diameter at breast height (DBH) of 22-32 cm for all species and a total range from 8 to 82 cm were collected in three campaigns: (i) February-March 2017; (ii) September 2017 and (iii) January-February 2018.The leaves were collected at five sites along the elevation gradient with Nyungwe Bigugu Mountain (~2700 m) being the highest and coolest site, followed by Nyungwe East (~2500 m), Nyungwe West (~1950 m), Cyamudongo (1800-1900 m) and the Ruhande Arboretum (1700 m).The three highest sites are located in Nyungwe National Park (TMF) while the Cyamudongo site is located in the adjacent Cyamudongo forest (TMF/LVTF) and Ruhande Arboretum (LTVF), close to Huye town.In total, 1080 leaf samples were collected: 3 leaves * 6 trees * 4 species * 5 sites * 3 campaigns.For detailed information, see [71].
For the Nyungwe TMF permanent plots, leaves were collected from trees of all species having an average DBH of 31-51 cm for all species and a total range of 7 to 96 cm during September-December 2013 and February-April 2015.The leaves were collected from 11 out of 15 permanent plots along a 32 km east-west transect at an elevation ranging between 1950 and 2500 m a.s.l established during late 2011 and early 2012.In total, 445 leaf samples were collected: 5 leaves * 89 trees.For detailed information about the plots, see [66].
All leaf samples were collected from mature and sun-exposed leaves from the upper canopy, following the same protocol for determination of leaf size and LMA as in the Rwanda TREE project.

Statistics
The site and species effects on variables measured both in 2018 and 2019 (leaf W/L, LA and LMA) were analysed using a mixed linear model, using site and species as fixed between-subject factors and year as within-subject factor (repeated measure).The site and species effects on variables measured in only one year (LVA, LD, N M , N A , P M and P A ) were analysed using a general linear model (GLM), using site and species as fixed factors.In both types of analysis, tree height was used as a covariate when it was significant.Since tree height changed with time (Table S3) we used a repeated covariance and applied a first-order autoregressive structure with heterogeneous variance in the mixed model analysis.For significant interactions of year by species and site by species, a one-way ANOVA was used to analyse the site effect on different species individually, with a Bonferroni test for post hoc comparison of individual sites.For species where the response variable was significantly correlated to tree height, a regression analysis with tree height as a covariate was used.The site effect was analysed by using plot means (18 per site) as replicates for morphological parameters, while block means (3 per site, i.e., 6 plots per block) for each species were used as replicates for nutrient parameters.To test for the differences between successional groups at the HE (control) site, a one-way ANOVA (repeated measure if both years was included) was used.Site means of each species were used as replicates within the early (10 species) and late (8 species) successional groups.The effect of tree age (juvenile trees Forests 2022, 13, 219 9 of 24 in this study compared to mature trees from other studies) on LA and LMA was tested for a common elevation of 2000 m a.s.l., using a GLM with elevation as covariate and tree age as a fixed factor.The significance of the relationship between LMA or LA versus elevation was tested using regression analysis.A test of normality (Shapiro-Wilk) was performed to test the distribution of data, and Levene's tested for homogeneity of variance.Homogeneity of variance was obtained for analysis of individual species, but not always when all species were analysed together.Effects were considered statistically significant at p < 0.05 if homogeneity of variance was obtained and at p < 0.01 for main effects when homogeneity was not obtained.All statistical tests were made using the SPSS 27.0 software package (SPSS, Inc., Chicago, IL, USA).

Leaf Trait Interrelationships and Variation between Species
Leaf size, leaf W/L ratio, LMA, LVA, LD, N M , N A , P M and P A varied greatly across species, as shown for the HE (Sigira) site (Figure 1a-f, Table S1).The difference between the two years was small and only significant for a relatively small difference in LMA within the ES group (11% lower in 2018 compared to 2019).The successional strategy of the species contributed largely to the inter-specific variation, as a significant difference between ES and LS species was observed for most variables except LVA, N A and P A content.In general, ES species had significantly higher leaf W/L ratio, LA, N M and P M, and lower LMA and LD compared to the LS species.D. torrida was often in the extreme position within the ES group while A. falcatus and F. saligna (depending on variable) shared the extreme position within the LS group.The species within a successional group that overlapped with species within the other group were more variable, but F. thonningii, classified as LS species, often had ES features that sometimes even exceeded the range of the ES species (Figure 1a,c,e,f).
The relationships between LMA and LA and other leaf traits were characterized by plotting LMA versus leaf W/L ratio, LA, N M , N A , P M and P A (Figure 2a-f) and LA versus N M (Figure S1).The relationship followed the same pattern at all sites, but here we only present the details from the HE (Sigira) site.LMA decreased significantly with increasing leaf W/L ratio (p < 0.001), LA (p < 0.001), N M (p < 0.001) and P M (p < 0.05), while it increased with increasing N A (p < 0.05) and P A (p < 0.005).These patterns were consistent both within and across successional groups (Figure 2a-f).Leaf size increased significantly with increasing N M across ES species (p < 0.001), while no correlation was observed for LS species (Figure S1).
The degree to which LVA and LD contributed to the variation in LMA was investigated by plotting both the slopes and R 2 values of these relationships against each other (Figure 3, Figure S2 and Table S2).Based on these plots, three groups of species were identified where LMA was mostly affected by (i) LVA (4 species), (ii) LD (9 species) or (iii) equally by both LVA and LD (5 species).Species from both successional strategies and both forest types of origin were represented in all three identified groups, suggesting that neither successional strategy nor forest type are main factors for explaining the different influences of LVA and LD on LMA.     1).Blue squares, LMA mainly influenced by LD; green triangles, LMA mainly influenc by LVA; red circles, LMA influenced by both LD and LVA.p-value refers to the regression line.

Leaf Shape, LMA and Leaf Size Responses to Sites and Warming
To separate from possible confounding influences of non-temperature site diffe ences, we defined warming effects as a significant change in the same direction (↗ increa ing or ↘ declining effect) at both ME and LE sites compared to HE, or a significant chang at LE compared to HE sites but no difference between ME and HE (see Table 3).For sp cies with a decline or increase at ME but not at LE (→ no systematic temperature effec compared to HE, the change was not attributed to warming. A significant site effect on leaf shape was observed (Table 4), with leaf W/L ratio b ing on average 5.4% and 7.5% higher at the warmer ME and LE sites, respectively, com pared to the HE site.However, the site effect varied between species as well as betwee years (Table 4), and therefore the site effect of each species and year was analysed (Tab 3, Figure 4a,b).Positive warming effects were observed for four species each year, b only Chrysophyllum gorungosanum Engl and S. guineense were similarly affected in bo years.The only negative warming response of leaf W/L ratio was observed for Prun africana (Hook.f.) Kalkman in 2018.No main site effects were observed for leaf size, LM or LVA, but LD overall declined at warmer sites (Table 4).However, significant site * sp cies, year * species and/or year * site * species * tree height interactions indicated variab  2).Blue squares, LMA mainly influenced by LD; green triangles, LMA mainly influenced by LVA; red circles, LMA influenced by both LD and LVA.p-value refers to the regression line.

Leaf Shape, LMA and Leaf Size Responses to Sites and Warming
To separate from possible confounding influences of non-temperature site differences, we defined warming effects as a significant change in the same direction ( increasing or declining effect) at both ME and LE sites compared to HE, or a significant change at LE compared to HE sites but no difference between ME and HE (see Table 3).For species with a decline or increase at ME but not at LE (→ no systematic temperature effect) compared to HE, the change was not attributed to warming.
A significant site effect on leaf shape was observed (Table 4), with leaf W/L ratio being on average 5.4% and 7.5% higher at the warmer ME and LE sites, respectively, compared to the HE site.However, the site effect varied between species as well as between years (Table 4), and therefore the site effect of each species and year was analysed (Table 3, Figure 4a,b).Positive warming effects were observed for four species each year, but only Chrysophyllum gorungosanum Engl and S. guineense were similarly affected in both years.The only negative warming response of leaf W/L ratio was observed for Prunus africana (Hook.f.) Kalkman in 2018.No main site effects were observed for leaf size, LMA or LVA, but LD overall declined at warmer sites (Table 4).However, significant site * species, year * species and/or year * site * species * tree height interactions indicated variable site responses among species and years for LA, LMA, LD and leaf W/L ratio.Therefore, species-specific tests for individual years were conducted (Table 3), with tree height (TH) included as covariate for species where LA and LMA significantly correlated to tree height (Table 3 and Table S3).In 2018 and 2019, significant site effects on LA were observed for 13 and 10 species, respectively, of which seven species declined with warming in each year (Table 3, Figure 4c,d).Declines were more common in ES compared to LS species; five and seven out of the seven species with declines were ES species in 2018 and 2019, respectively.For LMA, significant site effects were observed for 11 and 12 species in 2018 and 2019, respectively, of which six and nine, respectively, declined with increasing warming (Table 3, Figure 4e,f).Contrary to the warming effects on LA, the effects on LMA were relatively evenly mixed between ES and LS species.A significant increase by warming was only observed for LMA of Harungana montana Spirlet in 2018 (Table 3).Although no significant main site effects were observed for LA and LMA, the across species average decreased with 23 and 35% for LA and 4.4 and 17.5% for LMA at the ME and LE sites, respectively, when compared to the HE site.
Table 3. Significance of warming effect on leaf shape (width-to-length ratio, W/L) leaf size (leaf area) and leaf mass per unit area (LMA) in 2018 and 2019, and on leaf volume per unit area (LVA) and leaf density (LD) in 2019.Different letters (a, b, c) denote significant differences between sites.H, M and L represent high-(Sigira), mid-(Rubona) and low-(Makera) elevation sites, respectively.Explanation of species codes is given in Table 2; SG, successional group (ES, early; LS, late); FT, forest type (TMF, tropical montane forest; LVTF, Lake Victoria transitional forest); C, covariate used in statistical analysis when significant; h, tree height.The arrows indicate the direction of warming effect of species with significant site effect (p < 0.05).Decrease with warming (down at both M and L or no/small effect at M and down at L compared to H elevation sites).Increase with warming (up at both M and L or no/small effect at M and up at L compared to H elevation sites).→ No warming effect (up at one site and down at one site or down/up at M but no effect at L compared to H elevation sites).* a statistically significant site effect was obtained in the ANOVA analysis, but not in the Bonferoni post-hoc test.

Leaf L/W Ratio
a positive relationship between the LMA effect size and PA effect size when the LE site (Makera) was compared to the HE site (Sigira).Since the PA effect size was more negative at the ME site compared to LE site, while the opposite was true for LMA, it is not likely that variable leaf PA could explain the warming effect.3 for detailed statistical results.Note the broken axis for leaf size.
The soil content of N and P (Table 1), as well as both mass-and area-based leaf N and P content, differed significantly between sites (Table 4).Linear regressions were used to test if these site differences of leaf N and P content could explain the site and 'warming' effect on LMA and leaf size (Figures S4 and S5).The only significant effect observed was a positive relationship between the LMA effect size and P A effect size when the LE site (Makera) was compared to the HE site (Sigira).Since the P A effect size was more negative at the ME site compared to LE site, while the opposite was true for LMA, it is not likely that variable leaf P A could explain the warming effect.Table 4. p-values for effects of year, site and species, and for leaf width-to-length ratio (W/L), leaf size, leaf mass per area unit (LMA) and effects of site, species for leaf volume per unit area (LVA), leaf density (LD) and leaf nitrogen per unit mass (N M ) and area (N A ), and leaf phosphorus per unit mass (P M ) and area (P A ). Tree height was used as covariate when significant (Table S3).df represents degree of freedom.We observed significant main effects of both species and sites as well as an interaction between species and sites for LD (Table 4).The species-specific analysis showed a decline in LD with warming in six species (Table 3, Figure 5b).Most (four species) of these six species were included in the group where LMA was primarily affected by LD (Figure 3, Table S2).For LVA, a significant main effect was only found for species (Table 4, Figure 5a).   1 for full names) measured at three sites of different elevation (HE, ME and LE, high, mid and low elevations, respectively) in 2019.The species are divided into four groups, depending on successional stage and forest type of origin.ES and LS, early and late successional, TMF, tropical montane forest and LVTF, Lake Victoria transitional forest.* indicates statistical differences between sites within species (p < 0.05); see Table 4 for detailed statistical results.Note that there was no significant species x site interaction for LVA, and therefore the site effect for individual species was not analysed.4 for detailed statistical results.Note that there was no significant species x site interaction for LVA, and therefore the site effect for individual species was not analysed.
To test if LD and/or LVA contributed significantly to the site effects on LMA, regression analyses between sites effects on LMA versus site effects on LVA and LD, respectively, were conducted (Figure 6).These analyses showed a significant positive relationship between site effects on LMA and LVA, when both the ME and the LE sites were compared to the HE site, while no significant effects were observed for LD.These results partly contradict the results from analysis of the main site effect using ANOVA (Table 4).A likely explanation is that the site effect on LD is strong but only for a few species, while the site driven effect on LVA is more subtle and has a more general influence on the site effects on LMA.The significant effect of tree height on LVA may also have contributed to conceal the site effect.

Leaf Size and LMA in Juvenile and Mature Trees at Different Elevations
We compared the average leaf size and LMA from 2018 and 2019 of the juvenile trees in present study with those from mature trees from other studies of four species (two ES and two LS) growing at different elevations (1700 to 2700 m a.s.l) (Figure 7a-h, Table 5).When standardised to an elevation of 2000 m a.s.l., a significant main age (juvenile vs. mature) effect on leaf size or LMA was observed in M. kilimandscharica, but not in the other three species (Table 5).The elevation effect was also analysed using regression analysis, separating age classes for M. kilimandscharica, but not for the other species (Figure 7).Similarly to the analysis of the juvenile trees alone (Table 3), there was a significant decline in LMA, with declining elevation for C. grandiflora, M. kilimandscharica (mature data only) and P. fulva, but not for S. guineense (Figure 7e-h).However, the LMA-to-elevation relationship for P. fulva was only significant when the observation at 2700 m a.s.l.(Figure 7g, Table 5) was excluded.The reason to exclude the highest elevation is that the site was on the edge of the elevation range for that species.For leaf size, there was no significant elevation effect except for mature leaves of M. kilimandscharica, which declined slightly in size with increasing elevation (Figure 6b, Table 5).Table 5. p-values for effect of tree age (juvenile and mature trees) on leaf size and leaf mass per area unit (LMA) of four species along elevation gradients using elevation as a covariate and standardised to 2000 m a.s.l (common for all gradients).The elevation effect was analysed by regression analysis, and for significant age * elevation interaction, a separate regression analysis was conducted for each tree age class (Figure 7a-h).For full name of species, see Table 2

Discussion
We report results from an experimental elevation-gradient study where leaf morphological traits were collected from 18 Afrotropical tree species with different successional strategies and originating from different vegetation zones.All responses to a warmer climate were different among species and between early-and late-successional species groups.For species and traits that responded to warming, we generally observed decreased leaf size and LMA and increased leaf width-to-length ratio at warmer sites.Below, we discuss sources of variation in leaf morphology in tropical trees, including inter-trait relationships, acclimation to warming, and effects of species identity, successional strategy and ontogeny.

Characterisation of Leaf Traits in Different Species and Successional Groups
The LMA range (~50-150 g m −2 ) observed across species in our study is comparable to those reported in studies of trees in tropical rain and deciduous forests reviewed by Poorter et al. [21] and in a study along an elevation gradient in the Neotropics [37].The variation in species leaf width-to-length ratio and leaf size varied 20-fold (~0.05 to 1) and 15-fold (30-440 cm 2 ), respectively.The interspecific variation in LMA as a function of leaf width-to-length ratio, leaf size and leaf N and P per unit mass also broadly agreed with observation from other studies [1,72].Furthermore, the ranges of LD and LVA in our study (see Figure S2) are comparable to those found in previous studies of evergreen and deciduous species [21].Similarly to other studies [21][22][23][24], the contributions of LD and LVA to within-species LMA variation differed greatly among species; in some species, LD contributed more than LVA, in other species LVA contributed more than LD, and in a third group of species, LD and LVA equally contributed to the LMA variation (Figure 3).Our observed differences between successional groups are in agreement with previous studies, where ES species showed more acquisitive traits (higher leaf W/L ratio, LA, NM and PM, and lower LMA and LD) compared to LS species [51][52][53][54].Overall, the range of leaf-morphology characteristics across the selected species in this study is typical for humid tropical evergreen and deciduous trees, although extreme values are lacking.

LMA Responses to Warming
Regardless of successional group, LMA declined with warming (i.e., increased with elevation; Table 3, Figure 4) for some species in our study, which confirms the first part of hypothesis #1.This is consistent with observations from both field studies and controlled experiments in the tropics as well as in several other biomes [1,13,15,21,24,36,73].Within a given species, LMA usually correlates with photosynthesis, since thicker leaves with more palisade parenchyma tissue tend to have both higher LMA and photosynthetic capacity [21,74,75].Indeed, the lower LMA at warm sites agrees with the decline in photosynthetic capacity observed in LS species in the present elevation gradient [76].
The results on the contribution by LD and LVA to warming-induced shifts in LMA are less consistent with previous research.We hypothesised that LD contributed more than LVA to shifts in LMA at higher temperatures, based on the global meta-analysis [21].Indeed, there was a species-dependent effect of site on LD (declines in six species) but not on LVA (Tables 3 and 4).However, when the site effects on LVA and LD were compared with site effects on LMA using regression analysis, we found a significant positive correlation for LVA only (Figure 6).This suggests that shifts in both LVA and LD contributed to warming-induced changes in LMA, but that the contribution by shifts in LVA is consistent across species, while the contribution by LD is strong, but only for a few species.Therefore, the second part of hypothesis #1, that mostly LD affected the warming response of LMA, was not confirmed.
LVA shifts were found to be important for LMA changes in a study on 42 tropical tree species from different habitats in south-east Asia, regardless of environmental factors causing these changes [26], as well as in an elevation-gradient study on six tropical tree species in south China [73].The conflicting results on LVA and LD contributions to warminginduced changes in LMA in these studies on tropical trees (including our study) compared to global datasets [21] may reflect mixed effects of warming and other factors, especially effects of VPD along elevation gradients.The effect of high VPD on leaf size (i.e., leaf size declines) by reducing the total number of epidermal cells per leaf and per leaf area [6] likely altered the relationships between LMA, LD and LVA at our warmer sites.Furthermore, de la Riva et al. [77] concluded that both LD and LVA varied in a study of 34 Mediterranean woody species along a water-availability gradient.Several studies have also attributed both LD and LVA to LMA variation as a result of several environmental factors, including warming, cooling, drought, aridity, light intensity and nutrients availability, which complicate the causal interpretation of LD and LVA contributions to LMA variation in response to warming for our three groups of species identified in Figure 3 [17,21,25,37,78,79].

Leaf Size and Leaf Width to Length Response to Warming
Leaf size significantly decreased with warming in several species, but not in all (Tables 3 and 5).This result was opposite to hypothesis #2.It is in contrast to a study on trees in a temperate rainforest along an elevation gradient in New Zealand [45] as well as an investigation of herbarium specimens from eight tree species in China sampled at different latitudes (~20 to 40 • N) and prevailing temperatures [46].The conflicting results on leaf size may be a result of interacting effects of temperature and water limitations on leaf size [44].Although the plants in this study were watered at all sites, the VPD was higher at lower elevations (Table 1).High VPD causes declines in xylem and leaf water potential, production of abscisic acid and stomatal closure responses, even if soil water is abundant [80].Our result is probably explained by hydraulic constraints on water transport and leaf growth [23,46,81].This is also supported by Li et al. [46], who concluded that precipitation was more important than temperature for the effect on leaf size.Smaller leaf size in a warmer climate also implies thinner leaf boundary layer and more efficient heat dissipation, mitigating heat stress under sunny and hot conditions [40][41][42][43].However, with the rather small changes in leaf size observed here, this effect was likely of less importance than that of hydraulic constraints.
The study by Li et al. [46] also observed that temperature was more important for leaf W/L ratio, and contrary to our observation, reported a decrease with increasing temperature.This suggests, together with the results from our 18 tree species, that the warming response of leaf W/L ratio is species-specific.

Responses to Warming in Species of Different Successional Groups
We found that leaf size typically declined with warming in ES species, but not in LS species, while the effects on both leaf W/L and LMA were fairly evenly distributed between ES and LS species (Table 3).Hypothesis #3 was therefore only confirmed for leaf size and not for LMA and leaf W/L ratio.However, there is much information supporting the theory that ES trees have higher plasticity in response to warming compared to LS trees [55,56].The opposite effect (i.e., lower plasticity in ES compared to LS trees) has been found for changes in light conditions [82].Whether the plasticity of leaf morphology in response to warming is higher in ES compared to LS species remains to be more extensively investigated.

LMA and Leaf Size in Juvenile and Mature Trees at Different Elevations
Mature and juvenile trees exhibited similar relationships between elevation (and thus temperature) and leaf size or LMA in three out of four species (Figure 7), thus mostly supporting hypothesis #4.When combining juvenile and mature tree data, LMA decreased with declining elevation, on average over the four species by 24 g m −2 per 1000 m.This is of similar magnitude to the decrease of 17 g m −2 per 1000 m observed for the communityweighted LMA of sun leaves down an Andes-Amazon elevation gradient [13].This is a surprisingly low difference in elevation response when accounting for the variability in the LMA responses among the species in our study, and the differences in species compositions between the study sites.
When standardised at the same elevation (2000 m a.s.l.), no age effect was found on LMA and leaf size, except for M. kilimandscharica.Our findings do not support the common observation of ontogenic changes in leaf traits in other studies, with mature trees from both temperate and tropical biomes having higher LMA and thicker, larger leaves than juvenile trees [83][84][85].Only our observation of lower LMA in juvenile compared to mature M. kilimandscharica trees is consistent with most previous research, while the larger leaf size in juvenile compared to mature trees of this species were instead opposite to the common observation.However, in a study of 51 tree species in south-east Asian rain forests, it was shown that juvenile trees could have larger, smaller and similar sized leaves compared to mature trees [58], which is more in line with our results.The larger variability of leaf size compared to LMA in our study, as evidenced by on average higher within-site coefficient of variation (33% compared to 15% for leaf size and LMA, respectively), could contribute to potential inconsistences.However, the overall similarity between the responses from young and mature trees provides support for using more frequently available data from young trees to predict climate change effects on leaf morphology of mature forest trees, but with caution, as the responses are very species-specific.

Conclusions
This elevation-gradient study demonstrated that effects of a warmer climate on leaf morphology are highly variable among species, and sometimes also between successional groups.Leaf size decreased at warm sites in ES species, but mostly not in LS species, while changes in LMA (decreases) and leaf W/L ratio (increases) were common in both groups.The similarity in the warming response of leaf morphology in juvenile and mature trees suggests that our results are also applicable to mature forest trees.The reduction in leaf size at warmer sites found in most ES species is contrary to results in most controlled warming experiments.This may be linked to water economy, with hydraulic constraints on plant water-transport capacity acting to decrease leaf size in air with higher temperatures and VPD in the field, but not in chamber experiments with small and constantly well-watered plants.The warming-induced reduction in leaf size in ES species only may imply that these are better than LS species in preventing overheating in a warmer climate.Lower LMA in a warmer climate in some species is likely to make their leaves less resistant to external forces, such as wind and herbivory.The large variability in leaf morphological responses to a warmer climate may therefore shift the competitive balance among species and between successional groups, potentially affecting the tree community composition of tropical montane forests in a warming world.

Supplementary Materials:
The following supporting information can be downloaded online at: https://www.mdpi.com/article/10.3390/f13020219/s1,Table S1: Data of leaf morphology at the high-elevation site during 2018 and 2019; Table S2: Regression information for leaf volume per unit area (LVA) and leaf density (LD) versus leaf mass per unit area (LMA); Table S3: Tree height of sampled trees and the relationship of leaf morphology and nitrogen contents versus tree heights;

Forests 2022 , 24 Figure 1 .
Figure 1.(a-f).Leaf characteristics of 18 species at high-elevation site (Sigira).(a) Leaf shape (i.e., leaf width to leaf length ratio: W/L), (b) leaf size, (c) leaf mass per unit area (LMA), (d) leaf volume per area (LVA), leaf density (LD), (e) leaf nitrogen (N) per unit mass (NM) and area (NA), and (f) leaf phosphorous (P) per unit mass (PM) and area (PA).Leaf shape, size and LMA (a-c) were measured in both 2018 and 2019.Leaf N and P were measured in 2018 and LVA and LD in 2019 only (d-f).The boxplot boundary indicates the 25th and 75th percentile, the error bars indicate the 95th and 5th percentiles, the black solid line and the red dotted lines inside the boxes indicate median and mean, respectively.The black dots show the highest and lowest species means, for which species codes (see Table 1) are indicated in the lower and upper part of the figure.Different letters (a-c) inside boxes indicate significant differences (p < 0.05) of the mean between early (ES) and late (LS) successional species and years when relevant.

Figure 1 .
Figure 1.(a-f) Leaf characteristics of 18 species at high-elevation site (Sigira).(a) Leaf shape (i.e., leaf width to leaf length ratio: W/L), (b) leaf size, (c) leaf mass per unit area (LMA), (d) leaf volume per area (LVA), leaf density (LD), (e) leaf nitrogen (N) per unit mass (N M ) and area (N A ), and (f) leaf phosphorous (P) per unit mass (P M ) and area (P A ). Leaf shape, size and LMA (a-c) were measured in both 2018 and 2019.Leaf N and P were measured in 2018 and LVA and LD in 2019 only (d-f).The boxplot boundary indicates the 25th and 75th percentile, the error bars indicate the 95th and 5th percentiles, the black solid line and the red dotted lines inside the boxes indicate median and mean, respectively.The black dots show the highest and lowest species means, for which species codes (see Table 2) are indicated in the lower and upper part of the figure.Different letters (a-c) inside boxes indicate significant differences (p < 0.05) of the mean between early (ES) and late (LS) successional species and years when relevant.

Figure 2 .
Figure 2. (a-f) Leaf mass per unit area (LMA, y axis) in relation to (a) leaf width-to-length ratio (leaf W/L ratio); (b) leaf size; (c) leaf nitrogen per unit mass (Leaf N M ); (d) leaf nitrogen per unit area (Leaf N A ); (e) leaf phosphorus per unit mass (Leaf P M ); (f) leaf phosphorus per unit area (Leaf P A ).Each data point represents species mean at each site.HE, High-elevation site (Sigira); ME, Medium elevation site (Rubona); LE, low elevation site (Makera); ES, early successional trees; LS, late successional trees.Lines, equation and R 2 represent the regression line for all sites and successional groups.p-value refers to the regression line.

Figure 3 .
Figure 3. Slope (a) and R 2 (b) for leaf density (LD) versus leaf mass per unit area (LMA) in relati to slope and R 2 for leaf volume per unit area (LVA) versus LMA for 18 species (for species code see Table1).Blue squares, LMA mainly influenced by LD; green triangles, LMA mainly influenc by LVA; red circles, LMA influenced by both LD and LVA.p-value refers to the regression line.

Figure 3 .
Figure 3. Slope (a) and R 2 (b) for leaf density (LD) versus leaf mass per unit area (LMA) in relation to slope and R 2 for leaf volume per unit area (LVA) versus LMA for 18 species (for species codes, see Table2).Blue squares, LMA mainly influenced by LD; green triangles, LMA mainly influenced by LVA; red circles, LMA influenced by both LD and LVA.p-value refers to the regression line.

Figure 4 .
Figure 4. (a-f).Site (i.e., warming) effect on (a,b) leaf width-to-length ratio (W/L), (c,d) leaf size and (e,f) leaf mass per area unit (LMA) in 18 species (see Table 1 for full names) measured at three sites of different elevation (HE, ME and LE, high, mid and low elevations, respectively) in (a,c,e) 2018 and (b,d,f) 2019.The species are divided into four groups, depending on successional stage and forest type of origin.ES and LS, early and late successional, TMF, tropical montane forest and LVTF, Lake Victoria transitional forest.* indicate statistical differences between sites within species (p < 0.05); see Table3for detailed statistical results.Note the broken axis for leaf size.

Figure 4 .
Figure 4. (a-f) Site (i.e., warming) effect on (a,b) leaf width-to-length ratio (W/L), (c,d) leaf size and (e,f) leaf mass per area unit (LMA) in 18 species (see Table 2 for full names) measured at three sites of different elevation (HE, ME and LE, high, mid and low elevations, respectively) in (a,c,e) 2018 and (b,d,f) 2019.The species are divided into four groups, depending on successional stage and forest type of origin.ES and LS, early and late successional, TMF, tropical montane forest and LVTF, Lake Victoria transitional forest.* indicate statistical differences between sites within species (p < 0.05); see Table3for detailed statistical results.Note the broken axis for leaf size.

Figure 5 .
Figure 5. (a,b).Site (i.e., warming) effects on (a) leaf volume per area unit (LVA) and (b) leaf density (LD) in 18 species (see Table1for full names) measured at three sites of different elevation (HE, ME and LE, high, mid and low elevations, respectively) in 2019.The species are divided into four groups, depending on successional stage and forest type of origin.ES and LS, early and late successional, TMF, tropical montane forest and LVTF, Lake Victoria transitional forest.* indicates statistical differences between sites within species (p < 0.05); see Table4for detailed statistical results.Note that there was no significant species x site interaction for LVA, and therefore the site effect for individual species was not analysed.

Figure 5 .
Figure 5. (a,b) Site (i.e., warming) effects on (a) leaf volume per area unit (LVA) and (b) leaf density (LD) in 18 species (see Table2for full names) measured at three sites of different elevation (HE, ME and LE, high, mid and low elevations, respectively) in 2019.The species are divided into four groups, depending on successional stage and forest type of origin.ES and LS, early and late successional, TMF, tropical montane forest and LVTF, Lake Victoria transitional forest.* indicates statistical differences between sites within species (p < 0.05); see Table4for detailed statistical results.Note that there was no significant species x site interaction for LVA, and therefore the site effect for individual species was not analysed.

Forests 2022 , 29 Figure 6 .
Figure 6.(a,b) The effect of (a) leaf volume per unit area (LVA) and (b) leaf density on leaf mass per unit area (LMA) between the mid-(Rubona) and low-(Makera) elevation versus high-elevation sites (Sigira).Each symbol represents one species, yellow, Rubona; red, Makera.Open and closed symbols, early (ES) and late (LS) successional trees, respectively.R 2 and p-values are given for the regression line of all species at a site.

Figure 7 .
Figure 7. (a-h) Leaf size and leaf mass per area unit (LMA) at different elevations for Cgr, Carapa grandiflora (a,e); Mki, Macaranga kilimandscharica (b,f); Pfu, Polyscias fulva (c,g); Sgu, Syzygium guineense (d,h) in three independent studies in Rwanda.Each data point shows mean ± SD at each site included in the studies.Blue triangle, mature-tree elevation-gradient study; green square, permanent monitoring plots in Nyungwe forest with mature trees; red circle, elevation-gradient study with juvenile trees (Rwanda TREE, this study); lines, equation and R 2 represents the regression line for all tree studies.p-value refers to the slope coefficient.

Figure S1 :
Leaf size in relation to leaf nitrogen per unit mass; Figure S2: Leaf volume per area (LVA) and leaf density (LD) in relation to leaf mass per area unit (LMA); Figure S3: Scanned fresh leaves of the 18 studied species; Figure S4: The effect of leaf nitrogen and phosphorus on leaf size effects between sites; Figure S5: The effect of leaf nitrogen and phosphorus on LMA effect between sites.Author Contributions: Conceptualisation, A.M., E.B., D.N., J.U. and G.W.; Data curation, A.M., E.Z. and G.W.; Formal analysis, A.M., M.E.D. and G.W.; Funding acquisition, G.W., AM and J.U.; Investigation, A.M., B.N. (Bonaventure Ntirugulirwa), E.Z., B.N. (Brigitte Nyirambangutse), M.M. and G.W.; Methodology, A.M., J.U. and G.W.; Project administration, B.N. (Bonaventure Ntirugulirwa), D.N., J.U. and G.W.; Supervision, E.B., D.N. and G.W.; Visualisation, A.M. and G.W.; Writing-original draft, A.M. and G.W.; Writing-review & editing, A.M., J.U., D.N., M.E.D. and G.W. All authors have read and agreed to the published version of the manuscript.Funding: We acknowledge funding from the Swedish Research Council (VR; grant 2015-03338) and the Swedish Research Council for Environmental, Agricultural Science and Spatial Planning (Formas; grant 2015-1458).AM acknowledges University of Rwanda-Sweden programme grant, Central Research Grants sub-programme managed under the UR Directorate of Research and Innovation.GW and JU acknowledge the Swedish strategic research area 'Biodiversity and Ecosystem services in a Changing Climate' (BECC; http:www.becc.lu.se (9 December 2021)).

p-Values for Two Years Data p-Values for One Year Data
3.3.Effect of Warming on LVA and LD, and Their Contribution to the Effect on LMA