**Drought Impact on Leaf Phenology and Spring Frost Susceptibility in a** *Quercus robur* **L. Provenance Trial**

### **Ivica Cehuli´ ˇ c 1, Krunoslav Sever 2, Ida Katiˇci´c Bogdan 2, Anamarija Jazbec 2, Željko Škvorc <sup>2</sup> and Saša Bogdan 2,\***


Received: 20 November 2018; Accepted: 8 January 2019; Published: 11 January 2019

**Abstract:** Research highlights: The susceptibility of oaks to late spring and early autumn frosts is directly related to their leaf phenology. Drought may alter the leaf phenology and therefore frost tolerance of oaks. However, the effects of drought on oak leaf phenology and frost resistance have not been thoroughly studied. Background and objectives: One of the consequences of climate change is an increase in the frequency of dry episodes during the vegetation period. Pedunculate oak (*Quercus robur* L.) is an economically and ecologically important forest tree species that prefers humid habitats. Therefore, knowledge of the impact of drought on this species is of great importance for the adaptation of forestry strategies and practices to altered environmental conditions. The aim of this study was to determine the impact of drought on leaf phenology and spring frost susceptibility in nine provenances. Materials and methods: One-year-old saplings originating from nine European provenances were used in the trial. The saplings were exposed to experimental drought and then re-watered in two subsequent years. Spring and autumn leaf phenology were scored. The trial was impacted by a late spring frost in the third year, and the resulting leaf frost injury was scored. The effects of drought treatment on the phenology and frost susceptibility of plants from the provenances were analysed. Results: Leaf phenology of plants from most of the studied provenances was significantly influenced by the drought treatment (*p* < 0.001). Drought induced a carry-over effect on flushing phenology, which was observed as delayed bud burst (from 0.6 to 2.4 days) in the second year and as advanced bud burst (from 0.1 to 6.3 days) in the third year. Therefore, opposite shifts in flushing phenology may be induced as a result of differences in the time span when plants sense water deficits. In contrast to flushing, autumn leaf phenology was unambiguously delayed following the drought treatments for all studied provenances (from 2.1 to 25.8 days). Differences in late frost susceptibility were predominantly caused by among-provenance differences in flushing phenology. However, the drought treatment significantly increased frost susceptibility in the plants (the rate of frost-injured plants per provenance increased from 3% to 78%). This higher susceptibility to spring frost was most likely caused by the advanced flushing phenology that resulted from the drought treatment in the previous year.

**Keywords:** pedunculate oak; drought; stress; memory; flushing; autumn leaf senescence; phenological shift; carry-over effect

#### **1. Introduction**

The hazards associated with climate change, including rising temperatures, decreasing precipitation, and an increasing frequency of extreme climatic events, are expected to intensify [1–3]. In general, the productivity of forest ecosystems is severely impaired by water availability, and drought

may induce episodes of large-scale tree decline in temperate forests [4]. Most likely, such declines will substantially increase the necessity for artificial regeneration of temperate forests. However, the increasing severity and frequency of droughts may also impact forest regeneration, causing unacceptably high seedling mortality rates. Late spring frost is an additional abiotic factor that strongly limits forest regeneration [5]. Therefore, it is necessary to enhance our understanding of the mutual impacts of drought and frost on forest reproductive material.

Pedunculate oak (*Quercus robur* L.) is a widespread European temperate forest tree species that generally prefers fertile and moist habitats [6], although it is also tolerant to arid habitats such as the forest-steppe ecosystems in eastern Europe (e.g., see [7]). Moreover, this species is one of the most economically valuable European hardwood tree species and is a climax species in forests that harbor high biodiversity. Because the species is mainly adapted to moist habitats, its drought tolerance will become increasingly important for its survival, especially in the presence of additional stresses such as frost, competition, pests, and diseases. Theoretically, pedunculate oak may adapt to stressful environmental changes at the individual, population, and plant community levels [2] through phenotypic modifications (i.e., plasticity), natural selection, and hybridization with related xerophilous species, such as sessile oak and pubescent oak [8,9]. In general, oaks are relatively tolerant to drought [10], as they have a xeromorphic leaf structure and root structure that can cope with variability in soil water availability; additionally, oaks display the ability to rapidly resume assimilation after periods of water deficiency [11–15]. However, intraspecific variation in responses to drought occurs among oaks from different provenances. Some of this variation has been attributed to adaptability driven by natural selection in the original provenance habitats [16,17], but in other cases, the original site-level climate of provenances has not been correlated with the responses of oaks to drought (e.g., [18]). Reduced above-ground biomass production, together with a shift towards root growth, is a drought response often observed in oak species [18–20]. Other drought responses by oak include changes in the morphological (e.g., reduced leaf area and chlorophyll content) and physiological (e.g., reduced stomatal conductance and water use efficiency) properties of leaves [21]. Although pedunculate oak does not appear to possess efficient responses to leaf frost injury [22], oak populations display adaptation to frost events in the form of differences in leaf phenology [23]. Additionally, occasional frost events (together with other factors) may have led to high within-population genetic variation in pedunculate oak populations, which are often characterized by various phenological types [24]. For deciduous tree species such as *Q. robur*, the timing of bud burst in spring and leaf senescence in autumn is very important in defining the period of carbon fixation and growth during the vegetation season as well as determining survival due to the increased frequency of climatic extremes such as frosts. Among the most important adaptive traits, phenology is thought to be one of the most affected by climate change [25]. However, it is necessary to further enhance our understanding of the impacts of interactive stressful climatic events on the phenology of deciduous forest tree species. The effects of drought on oak phenology (especially that of the pedunculate oak) have not been thoroughly studied and have been mainly described as earlier growth cessation under dry conditions [17,26], which is also visible as earlier radial growth cessation [27]. Interestingly, drought may induce an after-effect (i.e., a carry-over effect) in the subsequent spring that is indicated by advanced bud burst [26,28]. The same effect of water stress on leaf phenology was also recorded in beech [29]. In one report, sessile oak from different provenances showed a drought-induced carry-over effect in the subsequent spring, which was observable as a delay in the bud burst date [30]. To the best of our knowledge, this discrepancy in observed leaf phenological shifts induced by drought in oaks has not been explained, nor reported for other forest tree species. Here, we propose that bud burst shifts that occur as a carry-over effect induced by drought have opposite directions depending on the time span over which the plants sensed the water deficit. That is, if plants sensed the water deficit early in the vegetation period (i.e., in spring), then flushing would likely be delayed in the subsequent year; in contrast, if drought were sensed late (i.e., in summer), then advanced flushing would likely occur in the subsequent year. The effects of drought on oak frost resistance have been

even less thoroughly studied. In general, few papers address between-provenance genetic variation in oak frost susceptibility (e.g., pedunculate oak [17] and holm oak [31]). One study examined the effect of drought on oak frost hardiness (measured as resistance to winter freezing periods) [32], but it did not evaluate frosts that occurred in late spring or early autumn. As drought may alter phenology and therefore the susceptibility of forest tree species to late spring and early autumn frosts, it is necessary to improve our understanding of these interactive effects.

We established a short-term pedunculate oak provenance trial with seedlings originating from nine provenances along the eastern European north-south gradient, i.e., from Estonia to Italy. Plants from these provenances were experimentally treated to induce drought conditions during two growing seasons. As the trial was affected by a late spring frost in the third year, we had the opportunity to examine variation in oak frost susceptibility among provenances as well as the impact of drought on the frost susceptibility of these plants. The aims of the study were to determine the following: (1) the impact of drought on the flushing phenology of plants from various provenances; (2) the impact of drought on the autumn leaf senescence of plants from these provenances; and (3) among-provenance differences in the impact of drought on spring frost susceptibility.

#### **2. Materials and Methods**

#### *2.1. Plant Origin and Trial Design*

In autumn 2013, acorns were collected below crowns of at least ten randomly chosen mature trees in each of nine allegedly natural pedunculate oak stands. The provenances are located along a latitudinal gradient from Estonia to Italy (Table 1). Acorns were sent to the Croatian Forest Research Institute and sown in 0.5-litre polyvinyl chloride (PVC) pots filled with nursery substrate. In spring 2014, germinated seedlings were individually transplanted into 50-litre PVC pots filled with natural soil extracted from a local pedunculate oak forest (soil type: gleysol; pH = 7.6; soil texture: silty loam). The soil was homogenized and cleared of stones, leaves, sprouts, and twigs prior to filling the pots. The pots were kept outdoors during the first vegetation period; thus, they were exposed to natural local weather conditions at the nursery of the Croatian Forest Research Institute (45◦40 08.27" N; 15◦38 31.96" E; 141 m a.s.l.). In spring 2015, the pots were transferred to a greenhouse and arranged following a strip-plot experimental design. The first factor (water availability) had two levels: a control plot and a drought treatment plot (DT plot). The second factor (replication) had three levels. Each provenance was represented within each combination of the factors by various numbers of randomly placed plants (3–9 per provenance due to differences in acorn germination and seedling survival).


**Table 1.** Location and climatic characterization of the provenance sites.

Mean annual temperature (T), mean annual precipitation (PR) and average beginning of the frost-free period in spring (bFFP) were generated with ClimateEU software (version 4.63) (http://www.ualberta.ca/~{}ahamann/data/ climateeu.html) for the period of 1901 to 2009. The provenance abbreviations indicate their country of origin.

The trial included a drop-in irrigation system that allowed water availability in the main plots (i.e., control and DT plots) to be controlled. Plants in the control plot were constantly maintained in soil with a high volumetric soil moisture (VSM) content (i.e., 40%–50% VSM) during the growing period.

The VSM content was monitored by soil moisture sensors, which were set up in six representative pots (one in each combination of the factors) and connected to automatic meteorological stations.

#### *2.2. Treatments*

In the first growing period, the trial was established in a greenhouse that was automatically ventilated, thus preventing excessive heat and air humidity during the summer. Plants in the DT plot were deprived of water on March 30th, 2015 (Figure 1). The drought treatment was stopped, and the plants were re-watered on July 21st (Figure 1), a date when approximately 20% of the plants in the plot had visible leaf injury symptoms. All plants in the trial continued to regularly receive water until the end of the growing season. During the dormant period of the year (i.e., beginning on November 1st, 2015), the pots were transferred outside the greenhouse and exposed to normal local winter conditions. Because the plants had unexpectedly shown vigorous growth in the greenhouse and threatened to reach the ceiling, the trial was performed outside for the remainder of the study. To prevent the soil from receiving natural precipitation, each pot was covered with polystyrene panels that were shaped to fit around the plant stem and additionally sealed with polyurethane foam and duct tape.

**Figure 1.** The beginning and end of growing seasons in the control and drought treatment (DT) plots in three consecutive years (2015, 2016 and 2017). Open symbols indicate the mean bud burst date in the control (square) and DT (circle) plots, whilst closed symbols indicate the mean date of autumn leaf senescence in the control (square) and drought treatment (circle). The symbol × indicates the late frost event (on 21st April 2017). Grey polygons indicate the duration of drought treatments in 2015 (from 30th March to 21st July) and 2016 (from 20th June to 16th July).

In the subsequent growing period (in 2016), the plants in the DT plot were deprived of water beginning on June 20th. The drought treatment was stopped on July 16th (Figure 1), when almost 25% of the plants showed visible leaf injury symptoms. Afterwards, the plants in the DT plot continued to receive regular watering until the end of the growing season, similar to the way plants in the control plot were treated.

#### *2.3. Leaf Phenology Scoring*

Leaf flushing and autumn leaf senescence were studied from 2015 to 2017 (flushing only in the last year).

Leaf flushing phenological phases were scored on all plants in the trial on a 1–8 ordinal scale:


4—folded leaf visible;

5—leaf unfolding but not yet flattened, small;

6—leaves still relatively small but with flattened blades, blade edges bent downwards, withered, lighter green or reddish;

7—leaves appear developed, larger but still more tenderly structured (thinner) than fully developed leaves and lighter green or reddish;

8—leaves fully developed, darker green, thicker.

Autumn leaf phenological phases were scored on a 0–5 ordinal scale:

0—leaves completely green with no visible yellowing;

1—up to 25% of plant leaves brown;

2—up to 50% of plant leaves brown;

3—more than 50% of plant leaves brown;

4—more than 75% of plant leaves brown;

5—leaves shed.

The day of the year when each plant reached flushing phase 3 was recorded as the bud burst date, and the day of the year when each plant reached autumn leaf phenological phase 3 was recorded as the autumn senescence date. These individual plant data were used for subsequent statistical analyses.

#### *2.4. Spring Frost Injury Scoring*

A late spring frost affected the trial (as well as a much wider geographic area) on April 21st, 2017. After a warm period (the mean air temperature from April 1st to April 20th was 12 ◦C, with a maximum temperature of 29.5 ◦C), the air temperature suddenly dropped from the maximum of 13.6 ◦C at 3 p.m. on April 20th to −0.1 ◦C at 8 p.m. on the same day, and the temperature continued to decrease until it reached a minimum of −5.8 ◦C at 5 a.m. on April 21st, after which the temperature rose to −0.3◦C by 7 a.m. and continued to rise.

Several days later, there were visible leaf injuries on the plants in the trial. Leaf injuries appeared as necrotic leaf tissue visible as various shades of brown stretching from the margins of the leaf blade towards the midrib. On some leaves, these injuries were located closer to the apex, but on others, the injuries were located closer to the base. In general, there were no cases in which the whole leaf blade surface was necrotic. Additionally, in all cases, a rather low percentage of plant leaves were injured (on average, 5% of all plant leaves were partly injured). As a result, it was difficult to grade the plants based on frost injury; thus, we decided to record frost damage using a binary scale (0—plant has no visible injuries, and 1—plant has visible injuries).

#### *2.5. Statistical Analyses*

Individual plant data on bud burst date and autumn leaf senescence date were used for statistical analyses. All statistical analyses were performed in SAS/STAT 15.1 software, a free version of SAS University Edition by SAS Institute Inc., Cary, NC, USA. Descriptive statistics were performed using the MEANS procedure. Analyses of variance (ANOVAs) were performed using the MIXED procedure to determine the statistical significance of the factors (irrigation treatments, blocks, provenances, and the provenance by treatment interaction) according to the following linear model:

$$\mathbf{y}\_{\text{iklį}} = \boldsymbol{\mu} + \mathbf{T}\_{\text{i}} + \mathbf{B}\_{\text{j}} + \mathbf{P}\_{\text{k}} + \mathbf{T}\mathbf{P}\_{\text{ik}} + \varepsilon\_{\text{ijk}} \text{ (ANOVA model)}\tag{1}$$

where yijkl—individual value of a trait; μ—overall mean; Ti—fixed effect of treatment i, where i = 1, 2; Bj—random effect of block j, where j = 1, 2, 3; Pk—fixed effect of population k, where k = 1, 2, ... , 9; TPik—population by treatment interaction; and εijkl—random error.

Assumptions of residual normality and variance homogeneity were tested by using the Shapiro–Wilk test and Levine's test [33] with the GLM and UNIVARIATE procedures in SAS. Residuals were plotted as a function of fitted values to test for variance homogeneity, and the distribution of residuals was also tested. The significance of the irrigation treatments was tested by the Wald test and with the Satterthwaite approximation [34]. Student's *t*-test was performed to establish the significance of differences in the bud burst date between the different groups of frost-injured and uninjured plants.

#### **3. Results**

#### *3.1. Flushing Phenology*

As expected, in spring 2015 (prior to drought treatment), there were no statistically significant differences in the timing of bud burst between the control plot and the DT plot (F = 2.52; *p* = 0.1131). However, in spring 2016 (the year after the first drought treatment), the mean bud burst date of plants from almost all (except CR (OT)) provenances in the DT plot was delayed (1.5 days on average), ranging from 0.6 to 2.4 days per provenance (Figure 2). The drought treatment effect was statistically significant (F = 12.90; *p* = 0.0004). In spring 2017, the mean bud burst date of plants in the DT plot was advanced by 3 days on average (compared with the control), ranging from 0.1 to 6.3 days per provenance (Figure 2). The only exception was provenance ES, for which plants exhibited a delay in mean bud burst of 1.3 days (Figure 2). The effect of drought on the bud burst date was again highly statistically significant (F = 31.44; *p* < 0.0001).

**Figure 2.** Provenance mean bud burst dates in the treatments (control and drought), expressed as day of year (DOY) in 2016 and 2017. Vertical bars indicate the standard deviation of the mean values. Asterisks above bars indicate significant differences between treatments for the same provenance. Provanance abbreviations: ES—Estonia; LI—Lithuania; PL—Poland; HU—Hungary; CR (RE)—Croatia (Repaš); CR (KO)—Croatia (Koška); CR (KA)—Croatia (Karlovac); CR (OT)—Croatia (Otok); IT—Italy.

#### *3.2. Autumn Leaf Phenology*

Autumn leaf phenology was scored in 2015 and 2016 after one and two years of successive drought treatments were imposed on plants in the DT plot, respectively. In 2015, the plants exhibited an average delay in senescence of 8 days (ranging from 2 (provenance ES) to 13 (provenance PL) days) (Figure 3). In 2016, leaf senescence was again delayed, in this case by an average of 15 days (ranging from 3 (LI) to 26 (CR (RE) days) (Figure 3). Based on ANOVA, the leaf senescence date of plants in the DT plot was significantly delayed compared to that in the control in both years (F = 179.00 and 318.85; *p* < 0.0001, respectively). Based on Tukey's test, the mean autumn senescence of plants from only the most northern provenances (ES and LI) did not significantly differ between the DT plot and the control plot.

**Figure 3.** Provenance means of autumn leaf senescence in the treatments (control and drought) expressed as day of year (DOY) in 2015 and 2016. Vertical bars indicate the standard deviation of the mean values. Asterisks above bars indicate significant differences between treatments for the same provenance. Provanance abbreviations: ES—Estonia; LI—Lithuania; PL—Poland; HU—Hungary; CR (RE)—Croatia (Repaš); CR (KO)—Croatia (Koška); CR (KA)—Croatia (Karlovac); CR (OT)—Croatia (Otok); IT—Italy.

#### *3.3. Spring Frost Susceptibility*

The trial was impacted by a late spring frost occurring on April 21st, 2017. The proportion of frost-injured plants varied among provenances. More than 90% of plants from provenances CR (KO), HU, and CR (RE) were injured in the DT plot, though no plants from provenance LI were injured (Figure 4). Additionally, a larger proportion of damaged plants were observed in the DT plots for most provenances, indicating an adverse impact of drought on frost tolerance.

**Figure 4.** Proportions of frost-injured plants by provenance and treatment. Provanance abbreviations: ES—Estonia; LI—Lithuania; PL—Poland; HU—Hungary; CR (RE)—Croatia (Repaš); CR (KO)—Croatia (Koška); CR (KA)—Croatia (Karlovac); CR (OT)—Croatia (Otok); IT—Italy.

Based on the relationship between the state of flushing on April 21st and the incidence of frost injuries, plants that had entered flushing phases 5, 6, and 7 were injured (Figure 5). All plants that had entered flushing phase 7 were injured in both the control and DT plots (only a single plant in the DT plot had no visible symptoms of frost injury). Additionally, all plants in the DT plot that had entered flushing phase 6 were injured, whereas only 44% of plants in the same flushing phase had visible leaf injuries in the control plot. There was an obvious difference in the proportion of injured plants that had reached flushing phase 5 between the control and DT plots (Figure 5).

**Figure 5.** Total proportion of frost-injured plants from all provenances (in the control and drought treatments) that were in flushing phase 5 (leaf unfolding but not yet flattened, small), 6 (leaves still relatively small but with flattened blades, blade edges bent downwards, withered, lighter green or reddish), and 7 (leaves appear completely developed, larger but still more tenderly structured (thinner) than fully developed leaves and lighter green or reddish) at the time of the frost event (on 21st April 2017).

Because the frost-injured plants were those that had entered flushing phases 5, 6, and 7 on April 21st, we analysed the proportions of these plants from the studied provenances separately for the control and DT plots. These results are presented in Figure 6. The provenances with the largest proportion of frost-injured plants in the control plot (see Figure 4) were also those with the largest proportion of plants that had entered more advanced flushing phases, and vice versa (compare Figures 4 and 6).

**Figure 6.** Proportions of plants in advanced flushing phases at the time of the frost event on 21st April 2017 by provenance and treatment: 5 (leaf unfolding but not yet flattened, small), 6 (leaves still relatively small but with flattened blades, blade edges bent downwards, withered, lighter green or reddish), and 7 (leaves completely developed, larger but still more tenderly structured (thinner) than fully developed leaves and lighter green or reddish). CO—control; DR—drought treatment. Provanance abbreviations: ES—Estonia; LI—Lithuania; PL—Poland; HU—Hungary; CR (RE)—Croatia (Repaš); CR (KO)—Croatia (Koška); CR (KA)—Croatia (Karlovac); CR (OT)—Croatia (Otok); IT—Italy.

Additionally, the increase in frost susceptibility in the DT plot may be explained by more advanced (i.e., earlier) flushing of plants in that plot. For example, the proportion of frost-injured plants increased by >50% in the DT plot for provenances PL (from 0% to 78%), HU (from 39% to 96%), and CR (KO) (from 43% to 96%) (Figure 4). The same provenances exhibited the greatest shifts towards more advanced flushing on April 21st (Figure 6).

These results led to the question of why all plants that had reached flushing phase 6 (FP6) suffered frost injuries in the DT plot while only 44% of the FP6 plants were damaged in the control plot. Thus, we tested for significant differences in the mean bud burst date between control plants not showing frost injury symptoms and drought-treated plants as well as between frost-injured control plants and drought-treated plants. These results are shown in Table 2.

**Table 2.** Student's *t*-test comparing the mean bud burst date between the two groups of control plants (injured or uninjured) and DT plants (all frost injured).


\* flushing phase 6. DT = drought treatment.

Therefore, we compared two groups of control plants with drought-treated plants. Student's *t*-test showed that uninjured control plants had started flushing significantly later, but injured control plants had not. Although all these plants were recorded as being in flushing phase 6 on April 21st, these results (Table 2) indicate that plants in the DT plot were at a more advanced flushing state similar to that of the damaged control plants.

#### **4. Discussion**

#### *4.1. Leaf Phenology*

In our study, pedunculate oak plants that were experimentally exposed to drought exhibited delayed bud burst in the subsequent spring (Figure 2). However, after the second drought treatment, the mean bud burst date of the plants was advanced (Figure 2). This discrepancy (i.e., opposite shifts in bud burst induced by water stress) was also revealed in other studies. For example, Mijnsbrugge et al. [30] reported delayed bud burst of oak exposed to drought and Kuster et al. [28] and Spieß et al. [26] reported advanced bud burst of oak exposed to drought stress in previous years. Thus, based on our results as well as those from the above-mentioned studies, drought stress has a carry-over effect on oak flushing phenology, which causes it to shift. This effect of water stress on flushing phenology was also recorded in *Fagus crenata* Blume [29]. Plant environmental responses are epigenetically regulated. Various environmental signals and stresses can induce epigenetic modifications, thereby creating a flexible "memory" system for short or prolonged periods of time [35]. Changes in epigenetic marks in trees are accompanied by morphological and physiological changes in various processes such as ageing, organ maturation, and bud set or burst [36–38]. Therefore, drought stress may have triggered an epigenetic response ("memory") resulting in the observed carry-over effect on bud burst date in the oaks from the studied provenances. Why bud burst shifts may occur in opposite directions (i.e., towards an earlier or a later date) is an open question. There may be more than one molecular mechanism behind this epigenetic "memory" phenomenon. The complexity of signalling events associated with sensing drought stress in plants is well known (e.g., see [39]). Various chemical signals (e.g., reactive oxygen species, calcium, and plant hormones) participate in the induction of stress tolerance via transduction cascades and the activation of genomic re-programming [40]. For example, abscisic acid (ABA) plays an important role in epigenetic processes

such as in ABA-mediated responses to abiotic stress (see [41]). Therefore, the complexity of drought sensing and signalling pathways likely results in different epigenetic modifications and therefore phenotypic changes such as opposite shifts in bud burst date. The different time spans over which plants sense water deficits in the vegetation period presumably cause opposite carry-over effects regarding flushing phenology. Accordingly, in the first year of the experiment, we stopped watering the plants in the DT plot on March 30th, allowing them to sense the water deficit early in the vegetation period, i.e., during or shortly after flushing. In contrast, in the second year of the experiment, the plants were deprived of water in the second half of June after they had completely developed leaves. Plants from the studied pedunculate oak provenances responded with delayed bud burst following the first drought treatment (Figure 2). Combined with advanced bud burst following the second drought treatment (Figure 2), this indicated that the hypothesis in this study was supported. Of course, the molecular mechanisms underlying the epigenetic "memory" phenomenon are still poorly understood [35]; thus, it is not clear why bud burst shifts occur in opposite directions. We can only speculate that drought in different time spans induced different sensing and signalling pathways, which resulted in the observed effects on the flushing phenology of the studied oak plants. Addressing the mechanisms underlying the observed effects of drought on flushing phenology was beyond the scope of this study.

Leaf senescence was significantly delayed by drought in plants from all provenances in both the first and second year of the experiment (Figure 3). Delayed autumn leaf senescence is a known physiological response to drought and re-watering in oaks [9,30] as well as other tree species [42]. This phenomenon has been explained as the photosynthetic and growth compensation of trees during the post-drought period of recovery. Notably, leaf senescence in the most northern provenances (i.e., Estonia and Lithuania) did not significantly differ between the DT and control plots. Plants from the northern provenances may be more sensitive to environmental signalling at the end of the growing season, which may prevent them from overextending their growing season. Vitasse et al. [43] reported a similar result in a study on clinal variation in the leaf senescence of *Quercus petraea* (Matt.) Liebl. from provenances along an altitudinal gradient. Additionally, Deans and Harvey [44] reported a negative correlation between leaf yellowing and latitude in sessile oak, although this correlation was not as strong as the positive correlation between budburst date and latitude. Nevertheless, natural selection in harsher (colder) habitats appears to favour genotypes that respond with more rapid growth cessation in autumn. This relationship may be why the delayed leaf senescence of plants from the northern provenances was not as significant as that of plants from the more southern provenances.

#### *4.2. Spring Frost Susceptibility*

Spring frost affected the trial (as well as a much wider geographic area) on the morning of April 21st, 2017. Within a few days, symptoms of necrosis were observable on plant leaves. This incident enabled us to study variation in frost susceptibility among provenances and between the treatments in the trial. There were significant differences among the studied oak provenances in terms of spring frost susceptibility. The mean proportion of injured plants was strongly related to their state of flushing development on the particular day when frost impacted the trial. Thus, differences among provenances in the proportion of frost-injured plants were caused by differences in flushing phenology. Additionally, the difference in the proportion of frost-injured plants between the control and DT plots may also be explained by differences in flushing phenology. Plants that were in more advanced flushing phases suffered leaf injury, whereas plants that had not started flushing or were in a state of less advanced flushing did not suffer damage (Figure 5). This result confirms the report by Utkina and Rubtsov [24], who found greater resilience to spring frosts of later flushing pedunculate oaks.

The most frost-susceptible plants were clearly those in flushing phase 7 (almost fully developed leaves with completely flattened blades but still tenderly structured, light green or reddish) because all but one of these plants suffered frost injuries in both the control and the DT plots. Additionally, all plants in phenological phase 6 in the DT plot were frost injured (Figure 5). Pedunculate oak do

not appear to possess effective active molecular responses to frost; instead, avoidance is their major anti-frost "strategy". This lack of active frost responses was true at least for plants from the early flushing provenances that were exposed to the frost in 2017 and suffered leaf injury. The use of escape as an adaptive response to spring frost risks has been reported in sessile oak populations along an elevation gradient [45] as well as in various pedunculate oak populations [24,46,47]. However, the presence of a few exceptions (such as a single plant in flushing phase 7 that did not show visible symptoms of leaf injury) suggests that at least some pedunculate oak genotypes possess an efficient response to freezing at the molecular level. Such a molecular response might be due to alterations in enzymatic activity [48] and the accumulation of soluble sugars, as has already been reported for various *Quercus* species [49,50]; additionally, such a response might be due to the more efficient mobilization of non-structural carbohydrates (NSCs), as reviewed by Villar-Salvador et al. [51].

There was a large increase in the proportion of frost-injured plants in the DT plot (for most of the provenances), indicating an adverse impact of drought on oak frost susceptibility (Figure 4). This result raises the question of why plants in the DT plot (drought stressed) showed increased spring frost susceptibility. Based on our results, the main reason for the increased frost susceptibility of these plants was the earlier flushing caused by the drought treatment in the previous year. Drought stress in the previous year induced a clear carry-over effect, which resulted in the earlier flushing of plants from most provenances (Figure 2; year 2017). Therefore, most of the provenances in the DT plot exhibited an increased proportion of plants at more developed flushing phases on the day when the spring frost occurred. Nonetheless, it was not as obvious why all plants at flushing phase 6 (100%) in the DT plot were harmed by the frost, while only 44% of plants in the same flushing phase in the control plot showed visible injuries. Drought stress may have impacted plants in the DT plot, resulting in an even more reduced efficiency of molecular responses to freezing temperatures. For example, the accumulation and re-mobilization of NSCs and nitrogen (N) reserves play a very important role in frost resistance (see [51] and references therein). Drought stress in the current growing season may hinder re-mobilization capacity in the next growing season by decreasing seedling N and NSC content. [51]. Oaks use NSC and N reserves from storage tissues for their early spring growth before bud burst [52], and significant reductions in these reserves in spring have been reported for oaks [48]. Therefore, the higher proportion of frost-injured saplings in the DT plot may be explained in the following way: summer drought impaired the capacity to accumulate, store, and re-mobilize NSC and N reserves (responsible for frost resistance and vessel formation in early spring), which elevated the frost susceptibility of the saplings in the DT plot.

Another possibility was that not all plants in flushing phase 6 were in the exact same state of flushing on a particular day. Specifically, injured plants may have started flushing slightly earlier (i.e., were at a slightly more advanced flushing state) than the uninjured plants, despite all of the plants being scored as in the same flushing phase on that particular day. This explanation for the difference in the proportions of damaged phase-6 plants between the DT and control plots is supported by the test results (Table 2). Uninjured phase-6 plants in the control plot started flushing significantly later (had an earlier bud burst date, on average), while injured phase-6 plants did not (Table 2). Therefore, injured phase-6 plants in the control plot and phase-6 plants in the DT plot were at a slightly more advanced state of flushing on the day of the frost, which may have rendered them more susceptible to frost injuries. Furthermore, the flushing phenology scoring method, in which a plant was graded based on the most frequent phase, may have resulted in two plants in somewhat different states of flushing being assigned the same score. For example, if most buds indicated flushing phase 6, a plant was graded as being in phase 6; however, some newly emerging leaves may have already been in phase 7. Such plants may have sustained frost injuries on those few phase-7 leaves. Although plausible, this hypothesis should be further tested.

#### **5. Conclusions**

The results of this study confirmed the previously reported effects of drought stress on pedunculate oak leaf phenology, i.e., delayed autumn leaf senescence and a carry-over effect on bud burst date. Although the effect of drought on autumn leaf senescence was unambiguous (i.e., delayed senescence), the mean bud burst dates were shifted in opposite directions in the two analyzed years. Our study supports the hypothesis that opposite shifts in bud burst dates occur due to the different times in the growing season at which plants sense water deficit signals. Accordingly, if oaks sense water deficits early in the growing season (i.e., during or shortly after flushing), then their bud burst date will be delayed in the subsequent spring. In contrast, if the oaks sense water deficit signals late in the growing season (i.e., summer drought), then their bud burst will shift towards an earlier date in the subsequent spring.

The results strongly indicated an adverse effect of drought stress on pedunculate oak spring frost susceptibility. The adverse effect was most likely caused by a shift in spring leaf phenology. In other words, summer drought stress induced a carry-over effect, resulting in earlier mean bud burst in the subsequent spring, which made plants from the majority of the studied oak provenances more susceptible to late-frost leaf injuries. Importantly, the impacts of drought on autumn leaf senescence, i.e., delayed leaf senescence, most likely increase the susceptibility of oaks to early autumn frosts.

However, plants from the latest-flushing provenances were the least sensitive to the above-mentioned adverse effects of drought on frost resistance. Therefore, these provenances may be considered valuable sources of forest reproductive material for the regeneration of oak forests made necessary by an increase in the frequency of drought and frost events.

**Author Contributions:** Conceptualization, K.S. and S.B.; Data curation, I.C., A.J. and S.B.; Formal analysis, I. ˇ C., ˇ A.J. and S.B.; Funding acquisition, S.B.; Investigation, I.C., I.K.B. and S.B.; Methodology, K.S. and S.B.; Project ˇ administration, S.B.; Resources, I.C. and A.J.; Supervision, S.B.; Validation, I. ˇ C., I.K.B. and A.J.; Visualization, I.K.B. ˇ and Ž.Š.; Writing—original draft, I.C. and S.B.; Writing—review & editing, K.S., A.J. and Ž.Š. ˇ

**Funding:** This study was fully supported by the Croatian Science Foundation under grant number IP-2014-09-4686 (Phenotypic and Epigenetic Response to Drought Stress and Adaptability of *Quercus robur* L. Populations along a Latitudinal Gradient).

**Acknowledgments:** Thank you to all those who kindly supplied acorns from their local provenances for the trial, including Dr. Moica Piazzai (Monte Rufeno Nature Reserve, Italy); Dr. Krista Takkis (Estonia); Dr. Virgilijus Baliuckas (Forest Tree Genetics and Breeding Department of the Lithuanian Research Centre for Agriculture and Forestry); Czesław Kozioł (Le´sny Bank Genów Kostrzyca); Tomasz Miedzyrzecki and Bartlomiej Laski (Mi ˛ekinia Forest District, Poland); Goran Peri´c (Forest Office Koška); staff of the Forest Office Repaš; staff of the Forest Office Karlovac and Mladen Šimuni´c (Forest Office Otok). The trial was maintained by employees of Nursery Production and the Division for Genetics, Forest Tree Breeding and Seed Husbandry of the Croatian Forest Research Institute. Many thanks to all of them. Thank you to the three anonymous reviewers for helpful comments and suggestions.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article*

### **Differentiation and Non-Linear Responses in Temporal Phenotypic Plasticity of Seasonal Phenophases in a Common Garden of** *Crataegus monogyna* **Jacq.**

#### **Kristine Vander Mijnsbrugge 1,\* and Astrid Janssens <sup>2</sup>**


Received: 16 February 2019; Accepted: 19 March 2019; Published: 27 March 2019

**Abstract:** Phenology in perennial plants implies the temporal occurrence of biological events throughout the year. Heritable phenotypic plasticity in the timing of the phenophases can be of importance in the adaptation of woody species to a quickly changing environment. We observed the timing of bud burst, flower opening, leaf senescence and leaf fall in two successive years in a common garden of *Crataegus monogyna* Jacq. in Belgium, consisting of six local and five non-local provenances. Data were processed with cumulative logistic mixed models. Strong auto-correlation was present among the spring phenophases as well as among the autumnal phenophases, with spring phenophases being negatively correlated with fall phenophases. The strongest between-provenance differentiation was found for the timing of bud burst in spring, followed by flower opening and finally by leaf senescence and leaf fall. Warmer spring temperatures in March 2017 advanced the timing of bud burst, and to a lesser extent of flower opening, in all provenances compared to 2016. However, the advancement was non-linear among the provenances, with the lower latitude provenances being relatively less early and the higher elevation provenances being more late than the local provenances in this year. It can be hypothesized that non-local provenances display larger temporal phenotypic plastic responses in the timing of their spring phenophases compared to local provenances when temperatures in the common garden deviate more from their home-sites.

**Keywords:** phenology; leafing out; flowering; senescence; cumulative logistic regression; hawthorn; provenance trial; non-local populations; variance analysis

#### **1. Introduction**

Woody plants are sessile and perennial organisms that are characterized by long generation times and slow migration rates [1]. Therefore, it is expected that woody species can adapt relatively quickly to changing local environmental conditions [2]. The capacity for evolutionary change depends on the standing genetic variation in tree populations [3] and climate is a major driver of evolutionary change over longer time scales [4]. Understanding the genetic basis of complex polygenic traits in woody plants that are clearly influenced by climate is therefore a current challenge in forest genetic research. Because all individuals in a common garden share the same environment, any average difference in a trait between provenances of the same species has a genetic origin. The genetic variation in fitness-related traits is typically estimated in open pollinated progeny tests in common garden experiments, including the estimation of differentiation between provenances [5,6], and the study of this variation has its place in the field of quantitative genetics [7].

The seasonal cycle in deciduous woody plant species of temperate regions is characterized by the timing of bud burst and flowering in spring and by the timing of bud formation, leaf senescence and leaf shedding in autumn. These phenophases can deviate between genetically differentiated populations and are believed to be adaptive, responding to selection induced by environmental change [5,8,9]. The timing of the phenophases marking the beginning and the end of the yearly growing season maximizes the annual growth while minimizing the risk of frost damage in spring and autumn. Frost can damage woody plants at high fitness costs [10]. Late frosts in spring can damage the soft tissues of young leaves, whereas early frosts in autumn may cause early leaf abscission, hindering the resorption of nutrients. In addition, early frosts are also known to damage the cambial zone in trees [11].

The timing of bud burst and leaf unfolding vary considerably among tree species, responding to divergent climatic conditions [3,5]. Bud burst can be the main trait that is affected by climate mediated selection [3]. Selection on timing of bud burst can be sufficiently strong to counteract the homogenizing effect of gene flow [12,13]. Rising spring temperatures advance the onset of the growing season in many woody species of the temperate zone, prolonging the growing season and thus affecting plant productivity and the global carbon balance [14–16]. On the other hand, late frosts in spring may hamper early flushing species in a global warming scenario [17]. Because of its adaptive nature and the easy assessment from an early age onward, spring flushing is frequently evaluated in common gardens [5,18]. Within the same woody species, population differentiation for bud burst, as observed in common garden experiments, typically follows clines along gradients of elevation and/or latitude of the home-sites of the populations [5,19,20] possibly due to different temperature requirements that are genetically determined [21]. The genetic variation in timing of bud burst, together with divergent selection, is believed to have allowed tree species in temperate regions to occupy large distribution ranges [22]. Within populations, individual trees display variable timing of bud burst, which is suggested to have a genetic cause [23]. Individual trees in a population can therefore be categorized as early, intermediate and late phenological forms [24]. Neutral genetic marker analysis showed that late bud burst forms of *Fagus sylvatica* L. in natural populations in Poland displayed higher within-population genetic variation in comparison to the early forms, suggesting that late spring frosts shape the neutral genetic structure of the populations [24]. Finally, no relationship was found between the timing of cambial activity and the timing of bud burst in *Quercus robur* L. [23].

Flowering and the subsequent fruit formation are part of the sexual reproductive cycle in plants. Reproductive phenology is sensitive to environmental cues such as temperature, moisture and herbivory [25]. Divergent timing of flowering can stimulate assortative mating in populations of woody plants, reducing gene flow and promoting population differentiation [25]. The timing of flowering varies strongly among woody angiosperm species and can occur before, during or after bud burst. For instance, pollen emission concurs with leaf unfolding in oaks. Assortative mating through long-distance pollen flow is therefore suggested to interact with local adaptation of bud burst [3,26]. In general, timing of bud burst and flowering are most likely auto-correlated in temperate tree species [26,27].

The emergence and growth of new spring foliage in temperate deciduous trees relies strongly on the nutrients that were resorbed during the preceding leaf senescence [28]. The timing of leaf senescence is affected by both photoperiod and temperature [29–32] while the timing of bud burst is primarily influenced by temperature [33,34]. It shows less year-to-year variability in comparison with timing of bud burst and is concomitant with less favorable conditions for photosynthesis [29]. In 59 tree species, the timing of autumnal leaf senescence displayed a pattern according to the climatic clines of the home-sites of the studied populations, which was clearer than the pattern observed in bud burst timing [5]. But, in a common garden experiment of *Quercus petraea* (Matt.) Liebl. composed of populations derived from the same geographic region but from deviating elevations, the population differentiation for autumnal leaf senescence did not correlate with the elevation of the origin [3]. There is no consensus concerning the factors controlling the leaf senescence process [29,35]. Leaves that

emerged after a late spring frost in beech and oak displayed higher photosynthesis rates and a delayed leaf senescence in autumn, compensating for spring frost damage and demonstrating that long-lived trees can adapt their autumnal phenology depending on preceding productivity [36]. The spring phenophases of bud burst and flowering, and the timing of leaf senescence and leaf abscission are most probably strongly auto-correlated. Therefore, it is suggested that leaf fall can be used as a proxy for leaf senescence, assuming a certain time delay [30].

The range of phenotypes that a plant can express as a function of the environment is called phenotypic plasticity. Genetically controlled, heritable phenotypic plasticity has the potential to influence plant evolution [37]. Species that are able to adjust their phenological responses to warming spring temperatures by earlier bud burst or earlier flowering show better performance when compared with less responsive species [38]. Because of their longevity, phenotypic plasticity may play an important role in the adaptation of woody perennials to the predicted climate change [39]. Studying responses of provenances of woody species in a common garden not only allows assessment of the differentiation among the provenances that is shaped by divergent selection, e.g. [40]. Also, repeated observations in successive years allow the estimation of phenotypic plastic reactions to variable meteorological conditions in the garden over time. We planted a common garden in Belgium consisting of *Crataegus monogyna* Jacq., a common shrub species in western Europe, including local and non-local provenances, and studied the seasonal phenophases that mark the growth cycle of woody plants. Using neutral (non-adaptive) molecular markers, high levels of genetic diversity within populations, but low levels of population differentiation, were found in this species [41]. We hypothesized that: (i) the timing of bud burst, flower opening, leaf senescence and leaf fall differ between local and non-local provenances; (ii) the timing of the four phenophases display varying degrees of within-provenance and between-provenance variation; (iii) spring phenophases are strongly auto-correlated, as well as autumnal phenophases; and (iv) the non-local provenances respond in a non-linear way, in comparison with the local provenances, to the variable meteorological conditions in the common garden environment in two successive years (variable temporal phenotypic plasticity among the provenances) which can be related to the home-site conditions of the non-local provenances.

#### **2. Materials and Methods**

#### *2.1. Common Garden*

A description of the provenances, the seed collection, the growth of the plants and the planting of the common garden have already been reported in [42]. In short, the common garden consisted of six local Flemish provenances (northern part of Belgium), two Walloon provenances (southern part of Belgium), one provenance from the UK, one from Italy and one from Hungary (Table 1). The Belgian (Flemish and Walloon) provenances were collected by the authors and grown in the nursery of the research institute in Geraardsbergen, Belgium, whereas the UK, Italian and Hungarian provenances were grown in adjacent nursery beds in a Flemish forest nursery, under a sales contract. For the collection of the seeds in the Belgian populations, care was taken to collect from *C. monogyna*, excluding *Crataegus laevigata* (Poir.) DC. individuals and individuals with a putative hybridogenic background, by visually assessing the morphology of the leaves and the berries. Seeds were collected from at least 30 individuals from each population. Information on the seed collection in the commercial provenances was not available as it concerned commercial plant material without a certificate of provenance. Certificates of provenance are not compulsory for *C. monogyna* according to the Council Directive 1999/105/EC [43]. Therefore, the exact location of these commercial seed sources was also not available. In the beginning of 2008, planting stock of the 11 provenances were planted in a common garden in Londerzeel, Flanders, Belgium. The provenances were randomly mixed in a single tree plot design and planted with a spacing of 1.5 × 1.5 m. Mean monthly temperatures for 2016 and 2017 were acquired from the weather station in Groenendaal which is located at a distance of 30 km from the common garden site (Figure 1).


**Table 1.** Descriptive data of the provenances of *C. monogyna* in the common garden. No precise data on seed-stock populations were available for the commercial provenances. -: no data available.

**Figure 1.** Mean monthly temperatures in 2016 and 2017 in Groenendaal, Belgium.

#### *2.2. Scoring of the Four Phenophases*

Observations of the phenophases were performed on the shrubs in the common garden in 2016 and 2017. Bud burst, flower opening, leaf senescence and leaf fall were scored following two 6-level and two 5-level protocols, respectively (Table 2). For all phenophases, the whole shrub was evaluated visually and a mean score level was given. Bud burst was scored in 2016 on 8, 23, 31 March and 16 April, and in 2017 on 12, 21, 26 March and 2 April. Flower opening was scored in 2016 on 28 April, 4 and 11 May, and in 2017 on 26 April, 1 and 11 May. Leaf senescence and leaf fall were scored in 2016 on 21 October and 4 November, and in 2017 on 8 and 29 October.



**Table 2.** *Cont.*


#### *2.3. Statistical Analysis of Phenological Data*

All statistical analyses were performed in the open source software R 3.5.1. [44]. Models were fitted to examine the timing of the four phenophases in the two successive years. Each phenological trait (*Tph*) was a response variable and was modeled using cumulative logistic regression in the "ordinal" package [45], as the recorded observations were in an ordinal scale. The command "clmm" in the "ordinal" package models the probability (p) of having reached maximally a given level of the ordinal response variable. The score levels for leaf senescence and fall were defined in increasing order, and the probability was modeled for having reached maximally a given score level; e.g., having reached maximally a leaf senescence score of 3 was to be interpreted as the probability of having reached scores of 1, 2 or 3. The score levels of bud burst and flower opening were defined in decreasing order, so that the probability of having reached maximally a bud burst score of 3, for example, included the probability of having reached scores of 6, 5, 4or 3. In this way, this could be interpreted as having reached a score of at least 3. Mixed models were fitted as the phenological data concerned repeated observations on the same plants. The year of observation (*Y*, categorical variable) and the provenance (*P*, categorical variable) were present in the fixed part of each model, including an interaction term between these two variables. In the four phenological models, the local Flemish provenance FL1 (Table 1) was taken as the standard provenance to which the timing of the other provenances in the common garden was compared. Day (*D*, numerical variable) was added in the fixed part to account for the different observation days. The random part (random intercept) consisted of a unique shrub identity code (*ID*). The latter accounted for the repeated observations on the same plants.

$$\log\left(\mathbf{p}\_{\mathrm{Tph}}/(1-\mathbf{p}\_{\mathrm{Tph}})\right) = \mathbf{a}\_i - \beta\_{\mathrm{Y}}.\\\text{Y (fixed)} - \beta\_{\mathrm{P}}.\\\text{P (fixed)} - \beta\_{\mathrm{Y}\mathrm{P}}.\\\text{Y (fixed)} - \beta\_{\mathrm{D}}.\\\text{D (fixed)} - \mathbf{r}\_{\mathrm{ID}} \text{ (random)}$$

where α<sup>i</sup> was an intercept value indicating the passing from one level of the ordinal phenological response variable to the next. β*Y*, β*<sup>P</sup>* and β*<sup>D</sup>* were the estimated coefficients for the fixed covariates *Y*, *P* and *D*, and r*ID* was the random effect coefficient for all levels of the variable *ID*.

To be able to compare the timing of the four phenophases for all provenances, the days were calculated for which the different provenances in each phenophase had reached the same stage of phenological development. The DOY (day of the year) was calculated for which the probability for having reached a bud burst or flower opening score of at least 4 and having reached a maximal leaf senescence and leaf fall score of 3, attained 50% (*D*50%*PY*) in the observation years 2016 and 2017 and for every provenance. A *D*50%*PY* for a given provenance and a given year therefore indicated the day that half of the plants of this provenance had reached minimally (bud burst and flower opening) or maximally (leaf senescence and fall) a given score level of the respective phenophase in the respective year. This calculation was based on log(p*Tph*/(1 − p*Tph*)) being 0 for p = 50%. With 2016 as the standard level for the variable *Y*, to which 2017 is compared, the following formulas were used:

$$D\_{\mathbb{50}\% P2016} = (\infty\_{\mathbf{i}} - \beta\_P) / (\beta\_D)\_{\mathbf{i}}$$

$$D\_{50\% P2017} = (\alpha\_i - \beta\_{2017} - \beta\_P - \beta\_{2017P})/(\beta\_D)\_i$$

When the timing of a given phenophase for a given provenance differed significantly from the standard provenance FL1, the time lag was calculated between this provenance and the standard provenance. For 2016, this was inferred from the model with 2016 as the standard level of the categorical variable *Y*, whereas for 2017 this was inferred from the model with 2017 as the standard level for the variable *Y*. The time span between a provenance and the standard provenance from which it differed significantly, was calculated by subtracting the calculated *D*50%FL1*<sup>Y</sup>* from *D*50%*PY*. Pearson correlation coefficients were calculated between the timing of the phenophases in the two observation years using the *D*50%*P*<sup>2016</sup> and *D*50%*P*<sup>2017</sup> values.

A significant interaction term between provenance *P* and year *Y* in the model statistics indicated that the time span between the timing of a phenophase for a given provenance and the standard provenance in 2016 differed significantly from the respective time span between these two provenances in 2017. The significant interaction terms were visualized in reaction norm figures, in which the timing of a phenophase for a provenance was compared between 2016 and 2017. The slope of the line connecting the timing of a phenophase between 2016 and 2017 for a provenance with a significant interaction term in the model statistics, differed significantly from the slope connecting the timings for the standard provenance FL1.

To examine the relative variance in the timing of bud burst, flower opening, leaf senescence and leaf fall among the different provenances in the common garden, in comparison with the relative variance among the shrubs within a provenance, the four phenological models were adapted by moving the provenance variable (*P*) from the fixed to the random part.

$$\log(\mathbf{p}\_{T\text{ph}}/(1-\mathbf{p}\_{T\text{ph}})) = \alpha\_{\text{i}} - \beta\_{Y}.Y - \beta\_{D}.D \text{ (fixed)} - \mathbf{r}\_{P} \text{ (random)} - \mathbf{r}\_{\text{ID}} \text{ (random)}$$

The relative variance between the provenances (σ<sup>2</sup> *<sup>P</sup>*) and the relative variance between the shrubs within a provenance (σ<sup>2</sup> *ID*) were obtained from these models.

#### **3. Results**

#### *3.1. Timing of Bud Burst*

All provenances in the common garden burst buds earlier in 2017 compared with 2016 (the covariate year was significant in the model statistics, Table 3, Figure S1). Significant differences in the timing of bud burst were observed between several provenances and the standard Flemish provenance FL1 in 2016 (Table 3) and in 2017 (Table S1). The time spans between the timing of these provenances and FL1 were calculated for both observation years (Table 4). The Italian and Hungarian provenances flushed earlier in comparison to the Flemish provenances, whereas the Walloon provenances with a higher elevation in the home-sites, flushed later (Figure 2a, Table 3).

**Table 3.** Model statistics for the timing of bud burst and flower opening. The provenance FL1 and the year 2016 are the standard levels for the categorical variables provenance and year, to which the other year, 2017, and the other provenances are compared. DOY: day of the year. Provenance abbreviations are given in Table 1.



**Table 3.** *Cont.*

Significant results: \*\*\* *p* < 0.001; \*\* *p* < 0.01; \* *p* < 0.05.

**Table 4.** Differences in timing between the standard local provenance FL1 and the other provenances in the common garden in 2016 and 2017 for the four phenophases. Time spans are shown only for the provenances that differed significantly from the standard provenance FL1 in the models (Table 3, Table S1). Negative values indicate earlier timing of the phenophase, positive values indicate later timing. Provenance abbreviations are given in Table 1.


**Figure 2.** Modeled timing of bud burst for the different provenances in 2016 and 2017: (**a**) Modeled probability of having reached at least bud burst score 4; (**b**) *D*50%*PY* values for the provenances (*P*) with a significant interaction term in the model statistics, for the years (*Y*) 2016 and 2017. *D*50%*PY* values indicate the modeled day of the year (DOY) when half of the plants of a provenance attain a bud burst score of at least four.

A significant interaction term between provenance and year in the model indicated a significant relative change in time span (between the timing of bud burst for the respective provenance and the timing for the standard provenance FL1) between the two observation years (Table 3). A significant interaction term for a provenance was visualized as the slope of a line connecting the timings (expressed as *D*50%*PY* values) for this provenance in the two observation years differing from the slope of the standard provenance FL1 (Figure 2b). In 2016, the southern European provenances burst their buds about 10 days earlier than the Flemish provenance FL1, whereas this time lag between the southern European provenances and FL1 was reduced to about 6 days difference in 2017, resulting in less steep slopes for the provenances HO and IT compared to FL1 in Figure 2b (provenance abbreviations in Table 1). Whereas in 2017 the difference between the timing of bud burst in the southern European provenances and FL1 was reduced in comparison with 2016, the time lag between the Walloon provenance WA1 and FL1 increased in 2017 in comparison with 2016, resulting in a less steep slope for WA1 compared with FL1 in Figure 2b (provenance abbreviations in Table 1). The UK provenance and the Flemish provenance FL3 also displayed significant interaction terms (Tables 1 and 3, Figure 2b).

#### *3.2. Timing of Flower Opening*

In general, all provenances in the common garden opened their flowers earlier in 2017 compared with 2016 (the covariate year was significant in the model statistics, Table 3, Figure S2). The southern European provenances and the Walloon provenances differed significantly in the timing of this phenophase when compared with the standard local provenance FL1 in 2016 and in 2017 (Table 3, Table S1 and Figure 3a). In addition, the Flemish FL4 and the UK provenance differed significantly from FL1 in 2017 (provenance abbreviations in Table 1). The time spans between these provenances and FL1 were calculated for both observation years (Table 4). Similar to bud burst, the Italian and Hungarian provenances opened their flowers earlier in comparison with the Flemish provenances, whereas the Walloon provenances tended to flower later (Figure 3a, Table 3). All time spans for flower opening between FL1 and the provenances that differed significantly from FL1, were smaller than for bud burst in both observation years.

**Figure 3.** Modeled timing of flower opening for the different provenances in 2016 and 2017: (**a**) Modeled probability of having reached at least a flower opening score of four; (**b**) *D*50%*PY* values for the provenances (*P*) with a significant interaction term in the model statistics, for the years (*Y*) 2016 and 2017. *D*50%*PY* values indicate the modeled DOY when half of the plants of a provenance attain a flower opening score of at least four.

Significant interaction terms between provenance and year in the model (Table 3) are visualized in Figure 3b. The pattern was comparable to bud burst. In 2016 the southern European provenances opened their flowers about 6 to 7 days earlier than the Flemish provenance FL1. In 2017, this time lag between the southern European provenances and FL1 was reduced to about 2 to 5 days, resulting in less steep slopes than FL1 in Figure 3b. Similar to bud burst, the time lag between the Walloon provenances and FL1 increased in 2017 compared with 2016 (from 2 to 3 days in 2016 to about 5 days in 2017), also resulting in less steep slopes than FL1 in Figure 3b. One extra local provenance, FL4, and the UK provenance differed significantly in timing of flower opening from FL1 in 2017 (Table S1). The UK provenance displayed a significant interaction term (Table 3, Figure 3b), whereas among the Flemish provenances only FL3 had a significant interaction term, visualized as a less steep slope than FL1 in Figure 3b.

#### *3.3. Timing of Leaf Senescence and Leaf Fall*

For the timing of leaf senescence and leaf fall, fewer provenances differed significantly from the standard provenance FL1 in comparison with the timing of bud burst and flower opening (Table 5, Table S2, Figure 4, Figure 5a, Figure S3 and Figure S4). For the southern European provenances, leaf senescence in 2016 occurred 9 to 16 days later than the local provenance FL1, and around 12 days later in 2017 (Table 4). Timing of leaf fall was modeled for these provenances 6 to 9 days later than the local provenance FL1 in both years (Table 4). In 2017, the timing of leaf senescence in the UK provenance was 7.5 days earlier than FL1, whereas in 2016, the timing of leaf fall in the Flemish provenance FL5 was 4 days later (Table 4).

**Table 5.** Model statistics for the timing of leaf senescence and leaf fall. The provenance FL1 and the year 2016 are the standard levels for the categorical variables provenance and year, to which the other year, 2017, and the other provenances are compared. Provenance abbreviations are given in Table 1.



**Table 5.** *Cont.*

Significant results: \*\*\* *p* < 0.001; \*\* *p* < 0.01; \* *p* < 0.05.

**Figure 4.** Modeled timing of leaf senescence for the different provenances in 2016 and 2017. Modeled probability of having reached a maximal leaf senescence score of three.

**Figure 5.** Modeled timing of leaf fall for the different provenances in 2016 and 2017: (**a**) Modeled probability of having reached a maximal leaf fall score of three; (**b**) *D*50%*PY* values for the provenances (*P*) with a significant interaction term in de model statistics, for the years (*Y*) 2016 and 2017. *D*50%*PY* values indicate the modelled DOY when half of the plants of a provenance attain a maximal leaf fall score of three.

#### *3.4. Correlation and Variance Analysis*

Pearson correlation coefficients between the timing of the different phenophases in the two observation years revealed, in general, high and significant correlations. The spring phenophases bud burst and flower opening displayed the highest correlation coefficients (0.98 \*\*\* and 0.97 \*\*\* in 2016 and 2017, respectively, Table 6). In comparison, the autumnal phenophases, leaf senescence and leaf fall, displayed lower correlation coefficients (0.79 \*\* in 2016 and 0.9 \*\*\* in 2017, Table 6). Both spring phenophases were correlated with leaf senescence (−0.87 \*\*\* and −0.86 \*\*\* for bud burst and flower opening in 2016, and −0.82 \*\* and −0.75 \*\* in 2017, respectively, Table 6), with later spring flushing and flowering occurring with earlier leaf senescence. In comparison, the spring phenophases were less correlated with leaf fall (−0.7 \* and −0.65 \* for bud burst and flower opening in 2016, and −0.78 \*\* and −0.65 \* in 2017, respectively, Table 6). Intra-phenophase correlation coefficients, between 2016 and 2017, were obviously high, with bud burst displaying the highest correlation coefficient (0.97 \*\*\*), followed by flower opening (0.92 \*\*\*), leaf fall (0.89 \*\*\*) and finally leaf senescence (0.86 \*\*\*).

**Table 6.** Correlations between the timing of the four phenophases bud burst (Bb), flower opening (Fo), leaf senescence (Se) and leaf fall (Fa) in 2016 and 2017. Pearson correlation coefficients and corresponding *p*-values are indicated above and below the diagonal respectively. A correlation coefficient with a corresponding p-value below 0.001 is indicated in bold and is underlined, between 0.001 and 0.01 is in bold and between 0.01 and 0.05 is underlined.


The relative variability in timing of the phenophases between the provenances was compared with the relative variability between the different shrubs within a provenance (Figure 6). In 2016, the between-provenance variance was relatively highest in the phenophase bud burst, and lower in decreasing order in flower opening, leaf senescence and leaf fall. The relative within-provenance variance increased accordingly. In 2017, the between-provenance variance for bud burst was lower compared with 2016 but was still higher than the relative within-provenance variance in this year. For flower opening, the between-provenance variance in 2017 was lower compared with 2016 and was as high as the relative within-provenance variance in this year. For leaf senescence and leaf fall in 2017, the relative variances attributable to the differentiation between the provenances were lower than in 2016, and thus displayed correspondingly higher relative within-provenance variances. In general, the spring phenophases displayed the largest relative between-provenance variances, in comparison with the autumnal phenophases.

**Figure 6.** Relative between-provenance (between prov) and within-provenance (between geno) variance for the timing of bud burst, flower opening, leaf senescence and leaf fall in 2016 and 2017.

#### **4. Discussion**

#### *4.1. Timing of the Phenophases*

Our results showed that differentiation between local and non-local provenances of *C. monogyna* is present in the phenological traits marking the seasonality in woody plants, as observed in a common garden. Differentiation in phenological traits in common gardens has been found in many tree species [46,47]. As shown before [42], the southern European provenances burst their buds earlier (lower latitude), the Walloon provenances are later (higher elevation), and the UK provenance is similar to the local Flemish provenances (lower longitude). We also detected, although in a lower order of magnitude, differentiation for bud burst between certain local provenances in both observation years. This finding may be related to results from Danish common gardens consisting of local populations of insect-pollinated shrub species (*Cornus sanguinea* L. *Malus sylvestris* (L.) Mill. and *Rosa dumalis* Bechst.), where differentiation in bud burst on a local scale, with very little spring temperature deviation between the home-sites of the populations, was suggested to be driven not only by natural selection but also by neutral processes [48]. Despite the fact that timing of bud burst and flower opening are strongly correlated, the time spans between the timing of bud burst from, on the one hand the standard local Flemish provenance and on the other hand the non-local provenances from the Walloon region (higher elevation), southern Europe (lower latitude) and the UK (lower longitude), were larger in comparison to the respective time spans between the timing of flower opening. This phenomenon has already been detected in the studied common garden [42] and proved to be consistent for two additional observation years. As hypothesized before [42], the timing of flower opening may be less sensitive to natural selection and local adaptation due to two reasons. Shrubs start flowering only after several crucial years of establishment and seedling development, and a year of reduced reproduction due to an improper timing of flowering may be less detrimental to a woody plant than reduced growth (accompanied by unfavorable competition with neighboring plants) because of an improper timing of bud burst. Therefore, the timing of flowering may be more responsive to the local micro-climate.

Leaf senescence and leaf fall in the southern European provenances were delayed compared with the local provenances. Together with an earlier bud burst, this implies a longer growing season for these provenances (spring phenophases were negatively correlated with autumnal phenophases). It is questionable whether the longer growing season is advantageous or disadvantageous. In a reciprocal common garden experiment of *Populus fremontii* S. Watson, southern populations planted in colder climates set buds relatively later in comparison with the same genotypes planted in a common garden at their home-site [39]. This later bud set is described as an inability to avoid early autumnal frosts and is interpreted as non-adaptive phenotypic plasticity, possibly caused by a lesser sensitivity to photoperiod as a cue to initiate bud set in more southern populations [39]. Bud flush in the southern population was found to be later in the colder common garden compared to the home-site common garden (but still earlier than the local populations from the cold environment) and was therefore interpreted as adaptive plasticity [39]. In our experiment, the higher elevation provenances of the Walloon region displayed a later bud burst in comparison to the local provenances, but no earlier nor later leaf senescence and leaf fall. This finding is likely in line with [34] who found no correlation between temperatures of the source sites along an altitudinal gradient and timing of leaf senescence in a common garden for *Acer pseudoplatanus* L. and *Fraxinus excelsior* L.. This may be due to a higher sensitivity of leaf senescence to photoperiod, compared to that of bud burst [29,30]. Higher elevations in the Walloon part of Belgium imply a generally colder climate compared to the local climate of the common garden site, but a negligible difference in photoperiod. When compared with the timing of bud burst and flower opening, a relatively smaller contribution of the between-provenance variation in the variance analysis indicates a weaker population differentiation for the autumnal phenophases. These results are in line with findings of [39] and [40] who both found a stronger influence of genotypic effects for bud set in poplar compared with bud burst, whereas bud burst showed stronger population-level effects relative to bud set. Although occurring in the same vegetative organ of the plant, the correlation between the timing of leaf senescence and leaf fall was smaller compared with the correlation between the timing of bud burst and flower opening, the latter implying a correlation between a vegetative (leaf bud) and a generative (flower) organ, which may suggest less tight genetic control for the autumnal phenophases.

#### *4.2. Non-Linear Temporal Responses in Timing of the Phenophases*

Inter-annual variation over long time periods in the timing of bud burst in temperate tree species has been widely modeled and discussed, e.g., [49]. Still, it remains difficult to accurately predict bud burst on smaller time scales. We studied the temporal phenotypic plasticity in the phenological responses on a small time scale, i.e., the responses to local meteorological conditions in the common garden, by observing the variability of the phenological responses on the same shrubs in two successive years. The non-linear inter-annual response of timing of bud burst was expressed in five significant interaction terms between the variables provenance and year in the modeling analysis, including four non-local provenances. For flower opening we found a comparable number of significant interaction terms (one extra Walloon provenance) with the same provenances involved as for bud burst, most probably giving expression to the auto-correlation of both phenophases. For leaf fall, only two provenances displayed significant interaction terms in the models, both being non-local, and for leaf senescence there were none. The lesser sensitivity of the autumnal phenophases to plasticity compared with the spring phenophases can be due to a higher sensitivity to photoperiod as a stable cue to initiate these processes [40]. The relatively high presence of non-local provenances among the significant interaction terms in the phenological models may be indicative of their non-local origin. The non-linear temporal responses in non-local provenances can be interpreted as a reaction to the prevailing growth conditions that deviate from the home-site conditions that they are adapted to. Considering the timing of bud burst, the time span between the southern European provenances and the local standard provenance in 2016 was larger than in 2017. For the Walloon provenances we observed an opposite response. As temperature is a well-known determinant for the timing of bud burst, the warmer temperatures in March 2017 have likely advanced the timing of bud burst in all provenances compared to 2016. It is by now well-known that woody plants shift to earlier bud burst dates upon increased warming [50]. However, the advancement that we observed in the non-local provenances was not linear in comparison with the local Flemish provenances, with the southern European provenances being relatively less early and the Walloon provenances being relatively later than the local provenances in 2017. In transplant experiments of *Populus fremontii*, the magnitude of phenotypic plasticity in bud flush and bud set was found to be correlated with the home-site climates [39]. Planting the southern European provenances of *C. monogyna* in the common garden

in Belgium implied a transfer to a cooler environment, with the cooler spring of 2016 inducing a longer time span in the timing of bud burst with the local provenances than the warmer spring of 2017. The warmer spring in 2017 drove the higher elevation provenances to a longer time span in the timing of bud burst with the local provenances, compared with the cooler spring in 2016. Together, these results suggest a larger plastic reaction in the timing of bud burst in the common garden environment may imply that prevailing meteorological conditions deviate more from the home-site conditions.

#### **5. Conclusions**

In Europe, populations of woody species at the southern limits of their natural range have maintained relatively large levels of genetic diversity and are therefore considered as appropriate gene pools for assisted migration towards the north as a climate adaptation strategy [51]. On the other hand, results from Danish (high latitude in Europe) common garden trials with local populations for several shrub species indicate substantial genetic variation and evolutionary potential, questioning the need for assisted migration specifically for widespread and generally occurring woody species [48]. In addition, *C. monogyna*, as a key component of old hedgerows, displayed large levels of genetic diversity in nuclear and chloroplast markers in the UK (medium latitude in Europe) [41]. Although limited in time, our two-year study showed that all provenances of the common shrub species *C. monogyna* adjust their phenological responses to the prevailing temperatures and that non-local provenances tend to react non-linearly relative to the local provenances, with larger temporal spring plasticity coinciding with a larger difference in climatic conditions between home-sites and the common garden environment. Our results can be interpreted as an extra argument in the debate over assisted migration, suggesting that for widespread species planted stock of non-local origin may tend to "over-react" to variable environmental conditions at the site to which they are transplanted. This over-reaction in plastic response may be caused by the environmental conditions that deviate more from the optima these provenances are adapted to at their home sites. Our results therefore stress the importance of carefully reflecting on assisted migration projects, and considering at least the addition of local planting stock when deciding to transport populations of trees and shrubs over longer distances in anticipation of the predicted climate change.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1999-4907/10/4/293/s1, Figure S1: Scores of bud burst in the common garden on four observation days in the years 2016 (a) and 2017 (b). Provenance abbreviations are in Table 1. Descriptions of score levels are in Table 2; Figure S2: Scores of flower opening in the common garden on three observation days in the years 2016 (a) and 2017 (b). Provenance abbreviations are in Table 1. Descriptions of score levels are in Table 2; Figure S3: Scores of leaf senescence in the common garden on two observation days in the years 2016 (a) and 2017 (b). Provenance abbreviations are in Table 1. Descriptions of score levels are in Table 2; Figure S4: Scores of leaf fall in the common garden on two observation days in the years 2016 (a) and 2017 (b). Provenance abbreviations are in Table 1. Descriptions of score levels are in Table 2; Table S1: Model statistics for bud burst and flower opening. The provenance FL1 and the year 2017 are the standard levels for the categorical variables provenance and year, to which the other year, 2016, and the other provenances are compared. Provenance abbreviations are in Table 1; Table S2: Model statistics for leaf senescence and leaf fall. The provenance FL1 and the year 2017 are the standard levels for the categorical variables provenance and year, to which the other year, 2016, and the other provenances are compared. Provenance abbreviations are in Table 1.

**Author Contributions:** K.V.M. conceptualized the study; K.V.M. and A.J. performed the field work and the statistical analysis and wrote the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** We specifically express our gratitude to Stefaan Moreels for general technical assistance throughout the whole project, Wim Stevens for maintenance of the common garden and Joris Vernaillen for help with the scoring.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### **E**ff**ects of Lime Application and Understory Removal on Soil Microbial Communities in Subtropical** *Eucalyptus* **L'Hér. Plantations**

### **Songze Wan 1, Zhanfeng Liu 2, Yuanqi Chen 3, Jie Zhao 4, Qin Ying <sup>1</sup> and Juan Liu 1,\***


Received: 7 March 2019; Accepted: 13 April 2019; Published: 16 April 2019

**Abstract:** Soil microorganisms play key roles in ecosystems and respond quickly to environmental changes. Liming and/or understory removal are important forest management practices and have been widely applied to planted forests in humid subtropical and tropical regions of the world. However, few studies have explored the impacts of lime application, understory removal, and their interactive effects on soil microbial communities. We conducted a lime application experiment combined with understory removal in a subtropical *Eucalyptus* L'Hér. plantation. Responses of soil microbial communities (indicated by phospholipid fatty acids, PLFAs), soil physico-chemical properties, and litter decomposition rate to lime and/or understory removal were measured. Lime application significantly decreased both fungal and bacterial PLFAs, causing declines in total PLFAs. Understory removal reduced the fungal PLFAs but had no effect on the bacterial PLFAs, leading to decreases in the total PLFAs and in the ratio of fungal to bacterial PLFAs. No interaction between lime application and understory removal on soil microbial community compositions was observed. Changes in soil microbial communities caused by lime application were mainly attributed to increases in soil pH and NO3 –-N contents, while changes caused by understory removal were mainly due to the indirect effects on soil microclimate and the decreased soil dissolved carbon contents. Furthermore, both lime application and understory removal significantly reduced the litter decomposition rates, which indicates the lime application and understory removal may impact the microbe-mediated soil ecological process. Our results suggest that lime application may not be suitable for the management of subtropical *Eucalyptus* plantations. Likewise, understory vegetation helps to maintain soil microbial communities and litter decomposition rate; it should not be removed from *Eucalyptus* plantations.

**Keywords:** lime application; understory removal; microbial community; forest management; *Eucalyptus*

#### **1. Introduction**

Soil microorganisms play essential roles in regulating ecosystem processes and maintaining ecosystem functions and services, such as litter decomposition [1], nutrient cycling [2], primary production, and climate regulation [3]. Both biotic and environmental factors drive the activity, structure, and diversity of soil microbial communities, which are controlled by many factors including

edaphic conditions [4]. Soil microbes generally respond quickly to environmental changes. Although numerous previous studies have explored the soil microbial responses to environmental changes in various forest types, more quantitative research estimating the effects of forest management practices on soil microbial community compositions is needed [5,6]. Information from such studies is valuable in both theoretical and practical aspects, which will improve our understanding of the mechanisms underlying soil carbon and nutrient cycling in forest ecosystems and may support sustainable forest management.

Acrisol is a common type of soil, distributed mainly in humid tropical/subtropical and even warm temperate regions, such as Southern Asia, South China, and Southeast US [7]. Lime application is one of the most widely applied management practices for amending Acrisol [8]. In addition to lowering soil acidity, lime application could improve soil structure, and thus typically promote agriculture and forest productivities [4,9]. Lime application also affects soil biological and biochemical properties, including microbial activities and carbon and nitrogen mineralization [7]. Previous studies investigating the effects of liming on soil microbial community composition reported mixed results. For example, both Pennanen et al. (1998) and Pawlett et al. (2009) found that lime application significantly reduced bacterial phospholipid fatty acids (PLFAs) and did not affect fungal PLFAs [10,11]. While Bruneau et al. (2010) reported that the density of bacteria was increased by liming [12]. Kamal et al. (2010) found that lime application led to increases in the relative abundance of prokaryotes, bacilli, and actinomycetes but decreases in micromycetes [13]. On the contrary, Wang et al. (2007) reported an inhibitory effect of lime on actinomycetes that could facilitate soil carbon mineralization [14]. Although the acid tolerance of individual microbial species varies widely, it is well recognized that low pH inhibits the growth of soil bacteria but favors fungi that are more resistant to acidity. Therefore, the soil microbial community generally shifts from fungi to bacteria as the pH increases after lime application [15]. However, most of these studies focused on the soil microbial communities in agricultural soil [8,16]. We still have a poor understanding of how fungal and bacterial communities respond to lime application in forests, especially for *Eucalyptus* plantations in South China where soils are strongly acidic.

*Eucalyptus urophylla* S.T.Blake is a fast-growing species and is widely planted in South China. Due to the high demand for wood and fiber, the planting area of *E. urophylla* is expanding [17]. Within an *E. urophylla* plantation, the understory vegetation is typically dominated by a native fern *Dicranopteris dichotoma* (Thunb.) Bernh., which usually forms a mat-like dense layer beneath the open-canopy areas [6,17]. Several studies have revealed that understory vegetation plays important roles in maintaining both the aboveground and belowground biodiversity, the microclimate, and the nutrient cycling [1,6,18]. Furthermore, understory fern biomass can occupy up to 20% of the total biomass in *Eucalyptus* plantations [19], indicating that understory plants could make a substantial contribution to the regulation of soil processes. Although the number of such studies is growing recently, the importance of understory vegetation in driving ecosystem structure and functions remains unfocused in most ecological studies [20,21]. In practice, understory vegetation is usually removed to reduce competition for nutrients against plantation trees, which could facilitate the growth of tree seedlings [22] and cause significant changes in soil microclimates [20] or nutrient availabilities [23]. To date, however, only a few studies quantitatively determined the effects of understory removal on soil microbial community composition, especially in combination with other forest management practices, and the effects tend to be conflicting. Some studies reported significant alterations in the composition of soil microbial communities [6,24], while others did not [25]. Therefore, the effects of understory removal on soil microbial communities are highly context dependent and need to be further studied.

In the present study, a lime application experiment combined with understory removal was conducted in a subtropical *E. urophylla* plantation to examine the main and interactive effects of lime application and understory removal on soil microbial community composition. We hypothesized that: (1) Lime application would increase soil pH, leading to negative effects on fungal biomass but positive effects on bacterial biomass, hence, changing the microbial composition, because most fungi are moderately acidophilic but bacterial species could grow within wider pH ranges; and (2) understory removal would result in negative effects on soil microbial biomass and composition, due to reductions of belowground resource inputs and changes in soil microclimates.

#### **2. Materials and Methods**

#### *2.1. Site Description*

This study was conducted in three six-year old *E. urophylla* monoculture plantations at Heshan National Field Research Station of Forest Ecosystem (112◦50 E, 22◦34 N), Chinese Academy of Science (CAS), Guangdong Province, China. Each of the three plantations occupied about 1 ha. The distance between any two plantations exceeded 1 km. The climate in this region is typically subtropical monsoon, with a humid hot season from April to September and a dry cold season from October to March. The mean annual temperature and rainfall are 21.7 °C and 1700 mm, respectively. The soil is classified as an Ultisol developed from sandstone [26]. *Eucalyptus* plantations were reestablished in 2005 after all previous trees (*Pinus elliottii* Engelm.) had been cut, with saplings planted at a spacing of 3 m × 2 m (1650 trees per ha). The understory vegetation within the plantation was dominated (almost 100%) by *D. dichotoma*, accompanied by some rarely occurring species such as *Rhodomyrtus tomentosa* (Aiton) Hassk. and *Miscanthus sinensis* Andersson.

#### *2.2. Experiment Design*

A split-plot design was used and four treatments were replicated three times, with the lime application as the main-plot factor and understory removal as the sub-plot factor. In December 2011, two paired main-plots were established in each plantation, and within each main plot, a sub-plot was distinguished. The main-plot size was 10 m × 10 m, while the sub-plot size was 5 m × 5 m. Lime application and understory removal were arranged for the main plot and the sub-plot, respectively. Therefore, the four treatments were: Control without lime application and with intact understory (CK), understory removal (UR), lime application (LA), and lime application combined with understory removal (LUR). Around each plot, polyvinyl chloride (PVC) boards were embedded at a depth of 80 cm to exclude disturbance from roots outside the plots and to eliminate the movement of water and nutrients. For plots receiving lime application, lime was added to the soil surface on 15 December 2011 and 15 August 2012 in the form of pulverized lime (60 kg per 100 m2 as CaO in total per year). Understory vegetation (including the above- and belowground components) was physically removed from plots for the understory removal treatments (i.e., UR and LUR) prior to the start of the experiment. Any understory plant germinated during the experimental period was removed on a monthly basis.

#### *2.3. Soil Sampling and Analysis*

Soils were sampled on 0, 10, 60, 147, 275, and 377 days after the treatments were applied, with the day 0 sample as the baseline. The surface litter was carefully removed before sampling. In each plot, five soil cores (2.5 cm in diameter) at 0–10 cm depth from random sampling points were combined as one composite sample. After carefully removing the surface organic materials and visible roots, each composite sample was passed through a 2 mm sieve and divided into two subsamples; one subsample was stored at 4 °C for soil physical–chemical analysis, and the other was stored at −20 °C for PLFA analysis.

Soil pH was measured every month during the experimental period. Three soil cores at 0–10 cm depth were sampled in each plot and combined as a composite sample. Soil pH was determined in a 1:2.5 soil-water slurry using a combination glass electrode. Soil temperature was measured with the DS1922l temperature logger iButton® (Dallas Semiconductor Corp., Dallas, TX, USA) every 2 h at 5 cm soil depth in situ from December 2011 to December 2012. Soil moisture contents (SMC% g of water per 100 g of dry soil) was determined gravimetrically after the soil was dried for 48 h at 105 °C.

The soil microbial community was characterized by the phospholipid fatty acid (PLFA) analysis [27,28]. Briefly, field-moist, sieved soil samples equivalent to 8 g dry mass were extracted in

a chloroform–methanol–phosphate buffer (1:2:0.8 *v*/*v*/*v*), and the extracted lipids were fractionated into neutral lipids, glycolipids, and polar lipids on silica acid columns by successive elution with chloroform, acetone, and methanol. The methanol fraction (containing phospholipids) was subjected to mild alkaline methanolysis to transform the fatty acids into free methyl esters and was analyzed on a gas chromatograph (Hewlett-Packard 6890, Agilent, Santa Clara, CA, USA). Peaks were identified using bacterial fatty acid standards and MIDI peak identification software (MIDI, Inc., Newark, DE, USA). Concentrations of each PLFA were standardized relative to 19:0 internal reference concentrations. Fungal biomass was considered to be indicated by 18:2ω6,9 [28]; bacterial biomass was considered to be indicated by the following 10 PLFAs: i15:0, a15:0, 15:0, i16:0, 16:1w7, i17:0, a17:0, 17:0, cy17:0, and cy19:0; total PLFAs and the ratio of fungal to bacterial PLFAs were also calculated. All PLFAs were calculated as ng g–1 dry soil.

#### *2.4. Litter Decomposition Experiment*

Litter bags (26 cm × 15 cm with a mesh size of 1 mm × 1 mm) were used to measure *Eucalyptus* litter decomposition rates. Fresh *Eucalyptus* litter was collected from litter traps (1 m × 1 m nylon net) in the experimental plantation. Each litter bag contained 10 g oven-dried litter, and six litter bags were placed on the soil surface in each subplot of the four treatments in December 2011. Afterward, one litter bag in each subplot was harvested bimonthly until all the six bags were harvested. The retrieved litter decomposition residuals were collected and oven-dried at 65 °C to a constant mass, and then weighed for the mass. The litter decomposition rate was calculated as the bimonthly mass loss rate: Litter decomposition rate (%) = (mass loss bimonthly/initial weight of litter) × 100%, where the initial weight of litter in each litter bag was 10 g.

#### *2.5. Statistical Analysis*

Repeated-measure ANOVAs by the general linear model were employed to examine the effects of time and treatments throughout the experimental period. Soil physico-chemical properties and microbial biomass measured during the experiment were analyzed using linear mixed-effect models, with lime application and understory removal as categorical fixed effects, and replicated (*n* = 3) as a random effect. Redundancy analysis (RDA) was performed to determine the relationships between microbial communities (PLFA profile) and soil physico-chemical properties. The most discriminating soil property variables were selected by the "forward selection" procedure. Data were transformed (natural log, square root, or rank) when necessary to meet assumptions of normality and homogeneity of variance. Statistical significance was determined at *P* ≤ 0.05. All univariate analyses were performed in SPSS 18.0 (SPSS Inc., Chicago, IL, USA). RDA was performed with the CANOCO 4.5 software (Ithaca, NY, USA). The forward selection was implemented based on the Monte Carlo permutation (*n* = 499).

#### **3. Results**

#### *3.1. Soil Physico-Chemical Properties*

Contents of soil organic carbon (SOC), total nitrogen (TN), dissolved organic carbon (DOC), as well as soil pH, did not differ among plots before the treatments were applied (Table 1), indicating no spatial heterogeneity between the experimental plots. Soil NH4 <sup>+</sup>-N contents varied with time but did not differ among the four treatments during the experimental period (Table 2). Soil pH, however, was significantly higher in the LA treatment (4.19–4.56 units) than in the control (3.78–3.83 units) or the UR treatment (3.57–3.87 units) across the sampling events (Figure 1). Soil NO3 –-N contents were increased by LA, but decreased by UR. Soil moisture contents (SMC) were significantly lower in both LA and UR treatments than in the control. For DOC, there was a slight decline in the UR treatment, but the trend was minimal (*P* = 0.06). In addition, no interaction between LA and UR on any soil physico-chemical property measured in this study was observed (Table 2).


**Table 1.** Soil characteristics before treatments were applied.

SOC, soil organic carbon; TN, total nitrogen; DOC, dissolved organic carbon; Control, no lime application and no understory removal; UR, understory removal; LA, lime application; LUR, lime application with understory removal. Values are means (with SE in parentheses; *n* = 3). Statistical significance was determined at *P* < 0.05. ns, no significant difference between treatments using ANOVA.

'HF -DQ )HE0DU\$SU0D\-XQ -XO \$XJ6HS 2FW 1RY

**Figure 1.** Soil pH under different treatments during the experimental period. CK, control without lime application and with intact understory; UR, understory removal; LA, lime application; LUR, lime application with understory removal. Values are means ± SE, *n* = 3.

**Table 2.** Effects of time, lime application, understory removal, and the interaction of lime application and understory removal on soil moisture content, soil temperature at 0–5cm depth, dissolved organic carbon, and inorganic nitrogen content, *n* = 3.


Time, sampling time; LA, lime application; UR, understory removal; LA × UR, interactions between lime application and understory removal; SMC, soil moisture content; ST, soil temperature; DOC, dissolved organic carbon. Results are from a three-way factorial ANOVA, with factors of time (levels: Dec 2011, Feb 2012, May 2012, Sep 2012, and Dec 2012), lime application (levels: limed, not limed), and understory removal (levels: understory removed, understory not removed).

#### *3.2. Soil Microbial Community*

Soil microbial biomass (as indicated by the total PLFAs, i.e., the sum of fungal and bacterial PLFAs) ranged from 453 to 2701 ng g–1 under different treatment during the experimental period. Sampling time had a significant effect on total PLFAs, bacterial PLFAs, fungal PLFAs, and the F:B ratio (ratio of fungal to bacterial PLFAs) (Figure 2, Table 3). Across sampling time, the LA treatment significantly decreased both fungal and bacterial PLFAs and hence the total PLFAs, but had no significant effect on the F:B ratio (Figure 2, Table 3). The UR treatment, however, significantly decreased the fungal and

total PLFAs without affecting bacterial PLFAs, and thus reduced the F:B ratio (Figure 2, Table 3). There were no significant interaction effects of LA and UR on the composition of soil microbial communities during the entire experimental period. In addition, no interaction effect between sampling time and LA was observed, while the interaction between sampling time and UR on the total PLFAs was remarkable (Table 3).

**Figure 2.** Soil microbial phospholipid fatty acids (PLFAs) under different treatments. (**a**) Soil Total PLFAs; (**b**) Soil Fungal PLFAs; (**c**) Soil Bacterial PLFAs; (**d**) the Ratio of Fungal PLFAs to Bacterial PLFAs. UR, understory removal; LA, lime application; LUR, lime application with understory removal. F:B indicates the ratio of fungal to bacterial PLFAs. Values are means ± SE, *n* = 3.



Time, sampling time; LA, lime application; UR, understory removal; LA × UR, interactions between lime application and understory removal; F:B indicates the ratio of fungal to bacterial PLFAs. Results are from a three-way factorial ANOVA, with factors of time (levels: Dec 2011, Feb 2012, May 2012, Sep 2012, and Dec 2012), lime application (levels: limed, not limed), and understory removal (levels: understory removed, understory not removed).

#### *3.3. Relationships between Soil Physico-Chemical Properties and Microbial Community*

Redundancy analysis showed that soil microclimate (*P* < 0.01), soil DOC (*P* < 0.01), soil pH (*P* = 0.02), and soil NO3 –-N content (*P* = 0.02) were significantly correlated with whole soil microbial communities (Figure 3). Specifically, soil pH and SMC were correlated negatively to 18:2ω6,9, i16:0, 15:0, and a17:0 (Figure 3). Soil temperature was correlated positively to i17:0 and negatively to a15:0, i15:0, and 16:1ω7 (Figure 3). Soil NO3 –-N content and soil DOC, however, were correlated negatively to i17:0, but correlated positively to a15:0, i15:0, and 16:1ω7 (Figure 3). All the environmental data explained 92.4% of the total variance, with the first axis (Axis 1) explaining 66.3% of the variance and the second axis (Axis 2) explaining another 26.1% (Figure 3).

**Figure 3.** Redundancy analysis of soil microbial PLFA biomarkers. Ordination diagrams presenting species scores and environmental factor scores (vectors). ST, soil temperature at 0–5 cm depth; DOC, dissolved organic carbon; SMC, soil moisture content.

#### *3.4. Litter Decomposition*

One year after the treatments were applied, both LA and UR treatments significantly reduced litter mass loss (both *P* < 0.01). There was no significant interaction of LA and UR on the litter mass loss during the experimental period (Figure 4).

**Figure 4.** Mass of *Eucalyptus* litter remaining per litter bag under different treatments. Vertical bars represent SEs. LA, lime application; UR, understory removal. The inserted *P* values were from repeated-measure ANOVA.

#### **4. Discussion**

#### *4.1. The Impacts of Lime Application on Soil Microbial Community*

Liming is commonly used to improve environmental conditions for the development of acid-intolerant microbes, resulting in increased microbial biomass and activities [4,29]. In this study, however, we observed a reduction of soil fungal and bacterial PLFAs after lime application. Consistent with our results, some other studies have also found that lime application inhibited the soil microbial biomass in both field- and lab-based conditions. For instance, Pawlett et al. (2009) reported that liming decreased soil microbial PLFAs and bacterial PLFAs [11]. Lin et al. (2018) also found that some soil bacterial groups were consistently decreased by a 27-year lime application in South China, and lime application may potentially affect the N-fixation ability by altering the soil community structure of soil diazotrophs [30]. Lime application did not affect soil the F:B ratio in this study, due to the decreases in both fungal and bacterial PLFAs (Figure 2, Table 3). Likewise, Xue et al. (2010) observed that liming had minimal effects on the F:B ratio, though soil microbial community structure had been significantly affected [31].

The most likely reason for lime application inducing alterations in soil microbial biomass may be associated directly with the increased soil pH, as shown by the fungal biomarker 18:2ω6,9, and the bacterial biomarkers *a*17:0, *i*16:0, and 15:0 were correlating negatively with soil pH (Figure 3). It is well recognized that most fungi are moderately acidophilic [7], and therefore, the fungal activities would be inhibited and the microbial communities would typically shift from fungal to bacterial as soil pH increased after lime application. Furthermore, the signature of *a*17:0, *i*16:0, and 15:0 had dominated the decreased bacterial biomass with higher pH in liming soils in the present study. In addition, our results showed that lime application reduced the soil water content, which might in turn suppress the NO3 –-N leaching from soil. The increased inorganic N might have "toxicity effects or salt effects" on soil microbes and consequently suppress the microbial biomass [32]. In line with our results, studies from Wu et al. (2011) and Cox et al. (2010) also demonstrated that the microbial diversity and biomass were reduced by the increased N availability [19,33].

#### *4.2. The Impacts of Understory Removal on Soil Microbial Community*

Soil microorganisms are greatly affected by both the quantity and quality of resources, as well as the soil habitat. Soil sampling time had significant effects on soil microbial biomass and composition in this study (Figure 2, Table 3), probably due to the seasonal variations in soil temperature and/or moisture. The present study region is within a typical subtropical monsoon climate, with distinct wet–hot and dry–cold seasons. Over 83% and 66% of annual precipitation and radiation inputs to this region occur during the wet–hot season (from April to September). The temporal changes in soil microbial biomass and composition could be ascribed to the changes in moisture, soil temperature, and the substrate availability, as indicated by other studies [24].

In addition, soil sampling time had a significant interactive effect with UR on soil total PLFAs, which should be considered when conducting "before-after-control-impact" (BACI) experiments in the monsoon climate regime region. Our findings showed that understory removal significantly reduced total PLFAs, fungal PLFA, and the F: B ratio. The effects of understory plant removal on soil microbial communities have also been observed in other forest ecosystems. For instance, Winsome et al. (2017) reported that understory removal significantly reduced soil microbial biomass in a ponderosa pine plantation [34]. Tong et al. (2015) found that understory removal apparently decreased soil fungal PLFAs and the F:B ratio in a Chinese fir plantation [35]. Yin et al. (2016) found that the biomass of all microbial community groups (i.e., fungi, bacteria, arbuscular mycorrhizal fungi, and actinomycetes) decreased with decreasing understory herbs biomass in subtropical plantations, which suggested that understory vegetation would exert strong controls on soil microbial communities [36]. Zhao et al. (2012) and Giuggiola et al. (2018) also suggested that the indirect effects of understory removal on soil microclimate wound be the most likely factor affecting soil microbial community composition [20,37]. Furthermore, RDA showed that the increased soil temperature and decreased soil moisture derived from the understory removal treatment resulted in changes in soil microbial composition. This may be attributed to the fact that the soil microbial community was significantly correlated with soil temperature and moisture, and all of the environmental variables explained up to 92.4% of the total variance in soil microbial communities in this study (Figure 3).

It is generally accepted that soil microbial community is controlled by the bottom-up force [38,39]. Normally, understory vegetation is considered as a minor portion of biomass in a forest ecosystem and contributes less to the bottom-up force on soil decomposers than overstory plants. Wu et al. (2011) reported that the DOC remained unaffected by understory removal in both two-year and 24-year *Eucalyptus* plantations [19], which indicated that the bottom-up control induced by understory vegetation was minor. However, there was an apparent trend that understory removal decreased DOC (*P* = 0.06) (Table 2), and thus resulted in the decrease in fungal PLFAs in the six-year *Eucalyptus* plantation in this study, which was proven by the positive correlation between the fungal PLFA biomarker (i.e., 18:2ω6,9) and the DOC (Figure 3). Additionally, understory removal decreased soil moisture (*P* < 0.01), and subsequently increased the fungal PLFAs, as indicated by the negative relationship between the fungal PLFA biomarker 18:2ω6,9 and soil moisture (Figure 3). Consequently, the suppressing effect of understory removal on soil fungal PLFAs could be ascribed to the decrease in DOC. Furthermore, our previous study in the same experiment site found that understory removal significantly reduced the soil CO2 flux, and the biomass of understory fern occupied up to 20% of the *Eucalyptus* plantation [17], which indicated that the understory dominated by *D. dichotoma* might have a bottom-up control on soil microbial communities. Since the ecological functions of understory vegetation in a forest ecosystem largely depend on their species [21] and biomass [6,40], the differentiated results from different ecosystems indicate that the effects of understory vegetation on soil microbial communities are context dependent and require further studies.

#### **5. Conclusions**

In conclusion, soil microbial communities in this study were affected by both lime application and understory removal, but without any interaction. Lime application increased soil pH and soil NO3 –-N content, which decreased the soil fungal and bacterial PLFAs. Understory presence is favorable for sustaining soil microclimates and providing resources for soil microorganisms. Understory removal inhibited soil total PLFAs and fungal PLFAs, and hence altered the soil microbial communities. In addition, both lime application and understory removal significantly reduced the litter decomposition rates (Figure 4), which indicates the lime application and understory removal may impact the microbe-mediated soil ecological process (i.e., soil nutrient cycling). We propose that lime application may not be suitable for the management of subtropical *Eucalyptus* plantations, while the understory fern *D. dichotoma* functions as a facilitator in *Eucalyptus* plantations and should not be removed. In general, our findings in this study would be helpful for the development of management practices that optimize the yield and sustainability of subtropical plantations.

**Author Contributions:** S.W. and Z.L. designed the experiment. S.W., Y.C., and J.Z. performed the study data collection. J.L., S.W, and Q.Y. analyzed the results. J.L. and Z.L. reviewed the article. S.W. wrote the final article. All authors approved the final version of the manuscript.

**Funding:** This study was funded by the National Science Foundation of China (41867007 and 31500341), Guangdong Provincial Key Laboratory of Applied Botany, South China Botanical Garden, Chinese Academy of Science (AB2016003), the Open Foundation of the State Key Laboratory of Urban and Regional Ecology of China (SKLURE2016-2-5), and the Youth Innovation Promotion Association CAS.

**Acknowledgments:** We thank Shenglei Fu and Weixing Zhang from Henan University for their comments and discussions on the early version of the manuscript. We also thank Xiaoli Wang and Xiaolin Zhu for their assistance in samples analysis and data collection.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

### **E**ff**ects of Phosphate Solubilizing Bacteria on the Growth, Photosynthesis, and Nutrient Uptake of** *Camellia oleifera* **Abel.**

**Fei Wu 1,2, Jianrong Li 3, Yanliu Chen 1, Linping Zhang 1,2,\*, Yang Zhang 1, Shu Wang 1, Xin Shi 4, Lei Li <sup>5</sup> and Junsheng Liang <sup>5</sup>**


Received: 14 March 2019; Accepted: 15 April 2019; Published: 20 April 2019

**Abstract:** Phosphorus (P) is a necessary nutrient for plant growth and plays an important role in plant metabolisms; however, the majority of P in soil is in insoluble forms. Phosphate solubilizing bacteria (PSB) can convert the insoluble phosphates into plant-available forms and may have the potential for use in sustainable agricultural practices. This study examined the effects of two native PSB, namely *Bacillus aryabhattai* (JX285) and *Pseudomonas auricularis* (HN038), and a mixture of both strains (1:1) on the growth of *Camellia oleifera* Abel. seedlings. The results showed a significant promotion of the growth of *C. oleifera* plants by three inoculation treatments. All the PSB inoculation treatments could improve the leaf nitrogen (N) and P content and had positive effects on the available N, P, and potassium (K) content of rhizosphere soil. A co-inoculation of the two native PSB strains caused a synergistic effect and achieved the best benefit. In conclusion, *B. aryabhattai* and *P. auricularis* could be used as biological agents instead of chemical fertilizers for agricultural production to reduce environmental pollution and increase the yield of tea oil.

**Keywords:** phosphate solubilizing bacteria; nutrition; oil tea

#### **1. Introduction**

Phosphorus (P) is an essential nutrient for plant growth and development. Although a large amount of P is present in soil, the majority of it is unavailable to plants [1]. In agro-forestry practices, it is a fact that many problems arise in the application of phosphate fertilizer [2]. On one hand, phosphate rock is a nonrenewable resource and may run out in 50–100 years due to anthropogenic exploitation [3]. On the other hand, about 70% of the phosphate fertilizer in soluble forms that is applied to the soil is rapidly combined with cations such as Ca2<sup>+</sup>, Mg2<sup>+</sup>, Al3<sup>+</sup>, and Fe3<sup>+</sup> and converted to insoluble forms [4]. Therefore, improving the absorption and use of P by crops is of great significance from both the ecological and economical perspective [5].

Phosphate-dissolving microorganisms are microbial resources closely associated with plant nutrition and account for 10% of all soil microorganisms [6,7]. Phosphate solubilizing bacteria (PSB) are a group of these microorganisms that can transform insoluble P compounds into available forms

by secreting organic acids, and they may be used as inoculants to enhance P availability for plants [8]. In addition, they can also promote plant growth via producing hormones, such as cytokinin and indole acetic acid [9,10]. Despite some PSB present in plant rhizospheres and in soil, the amount of P released by these microorganisms is usually insufficient to meet the demand of growing plants [11]. High-efficiency PSB have the potential for making a great contribution to the decrease of environmental pollution and promoting ecological balance by replacing chemical fertilizers [12]. Consequently, there is an urgent need to investigate the effects of selected high-efficiency PSB on plant nutrition and growth.

*Camellia oleifera* Abel. (Theaceae), a unique edible oil tree species to China, is one of the world s four famous woody oil plants [13]. Tea oil obtained from seeds has an unsaturated fatty acid content of up to 90%, much higher than vegetable oil, peanut oil, and soybean oil, and contains special physiologically active substances, such as camellia, which have extremely high nutritional values [14]. *Camellia oleifera* generally grows in the mountains and hills of subtropical regions in southern China and is also of great value in soil and water conservation and maintaining ecological balance [15]. However, the acid soil in southern China has a low available P content and limits the growth and productivity of *C. oleifera* [16]. In our previous study, we screened two high-efficiency PSB strains (JX285 and HN038) from the rhizosphere soil of *C. oleifera* [17,18]. The present study's aim was to determine the impacts of JX285 and HN038 on the growth, photosynthesis, and nutrient uptake of *C. oleifera*. The results may provide a theoretical basis for the development of microbial fertilizers for use in agroforestry.

#### **2. Materials and Methods**

#### *2.1. Strains, Plant Material, and Growth Medium*

The two PSB strains, JX285 (*Bacillus aryabhattai*) and HN38 (*Pseudomonas auricularis*), were isolated from the rhizosphere of *C. oleifera* in our previous experiments [17,18]. The PSB inoculum was prepared as follows: the two PSB strains were individually grown in liquid Luria-Bertani (LB) medium (containing tryptone 10 g/L, yeast extract 5 g/L, and NaCl 10 g/L) with shaking at 180 rpm for 24 h at 28 ◦C, and then the broth was centrifuged at 10,000 r/min for 10 min. The supernatant was discarded, and the pellet was re-suspended and washed with sterilized water three times. The final concentration of PSB inoculum was adjusted to 108 colony-forming units (CFU) mL<sup>−</sup>1.

The seeds of *C. oleifera* were provided by the Jiangxi Academy of Forestry Research, China. The seeds were first surface sterilized with 0.3% potassium permanganate solution for 1 h and rinsed with sterilized water three times. They were then germinated on wet gauze at 30 ◦C. The germinated seeds were transplanted into 1 kg pots with sterilized sand. Finally, uniform seedlings were selected and transferred to a new pot containing 3 kg of growth medium. Each pot was planted with one seedling.

The growth medium was a mixture of krasnozem soil, sand, and peat soil (6:3:1, *v*/*v*). The soil used in the experiment was collected from the campus of Jiangxi Agricultural University (China) and was sieved after being air dried (1 mm). The growth medium was autoclaved at 0.11 MPa and 121 ◦C for 2 h. The basic physicochemical properties of the growth medium were as follows: pH 5.25 (soil:water—1:5), 43.71 mg/kg organic matter, 94.68 mg/kg available nitrogen (N), 2.80 mg/kg available P, and 15.80 mg/kg available potassium (K).

#### *2.2. Experimental Design and Growth Condition*

A pot experiment was conducted using a two-factor completely randomized block design containing five inoculation treatments and three P treatments [19]. The inoculation treatments were as follows: inoculated with *P. auricularis* HN038 (T1), inoculated with *B. aryabhattai* JX285 (T2), inoculated with the mixture of HN038 and JX285 (1:1) (T3), inoculated with LB medium (CK1), and inoculated with sterile water (CK2). The P treatments were as follows: no phosphate fertilizer (0 g/kg), (P1), added 3 g calcium superphosphate (1 g/kg) (P2), added 6 g calcium superphosphate (2 g/kg) (P3). Each treatment contained 5 replicates. Furthermore, in order to prevent the plant grow from stopping due to a lack of nutrient requirements, urea and potassium nitrate were added to the soil to reach

300 mg/kg N and 200 mg/kg K, respectively. The pot was placed in a greenhouse with 12 h of light per day. After five months of growth, the growth and photosynthetic parameters were measured and the seedlings were harvested.

#### *2.3. Plant Growth Measurement*

Three healthy seedlings with consistent growth were selected for each treatment to measure their growth parameters. A measuring tape (Swordfish, China) was used to measure the plant height. The leaves, stems, and roots were harvested and dried at 75 ◦C to a constant weight and weighed to calculate the total biomass of each plant.

#### *2.4. Measurement of Gas Exchange Parameters and Relative Chlorophyll Content*

A LI-6400 portable photosynthesis system (Li-Cor, Lincoln, NE, USA) was used to determine the photosynthetic parameters. When the seedlings were growing vigorously (mid-September) and the temperature was relatively stable (09:00–11:00 am), healthy leaves in the upper part of the seedlings were chosen to determine the net photosynthetic rate (Pn) and the transpiration rate (Tr). The water use efficiency (WUE) was calculated as WUE = Pn/Tr. A chlorophyll meter (SPAD-502) was used to measure the relative chlorophyll content of the plants. A random measurement method was used between each treatment to eliminate the inevitable errors caused by different measurement times.

#### *2.5. Measurement of Plant Nutrient and Soil Nutrient*

The dry leaves were finely ground and homogenized to determine the N and P concentrations. The P content was assayed by the dry ash digestion method [20]. The N content was determined by the Kjeldahl digestion method [21].

Samples of rhizosphere soil were collected to test the soil's nutrient status. The N content was measured using the alkaline hydrolysis method [22]. The P content was determined by the molybdenum ruthenium anti-colorimetric method [23]. The K content was measured by atomic absorption spectrometry [24].

#### *2.6. Data Analysis*

Statistical tests were performed using SPSS 13.0 (SPSS Inc., Chicago, IL, USA). A two-way ANOVA was performed to analysis the significance of P application, the PSB inoculation treatment and the interaction of P × PSB treatment. A one-way ANOVA was employed to examine the differences among the different inoculation treatments under each P level. Duncan's multiple range test was performed at *p* = 0.05 in case of significant impact by factor.

#### **3. Results**

#### *3.1. Plant Height and Biomass*

The two-way ANOVA results showed that plant height and biomass were significantly (*p* ≤ 0.01) influenced by the interaction between P application and PSB inoculation (Table 1). The plant height and biomass were significantly (*p* ≤ 0.05) improved by the single and mixed inoculation of two PSB strains compared to the non-inoculation control under each P level (Table 1). Co-inoculated plants showed a greater plant height and dry weight than those of plants inoculated with single PSB strains and non-inoculated plants. There was no significant difference in plant height and biomass between the CK1 and CK2 treatments. Moreover, the degree of the increase by PSB inoculation differed to some extent under different P levels. The maximum positive effects of the PSB inoculation on plant height and biomass were observed at the P3 level.


**Table 1.** Effect of PSB on plant height and biomass of *C. oleifera* seedlings.

Data are means ± SD (*n* = 3). P: phosphorus; PSB: phosphate solubilizing bacteria; P1: no phosphate fertilizer (0 g/kg); P2: added 3 g calcium superphosphate (1 g/kg); P3: added 6 g calcium superphosphate (2 g/kg); T1: inoculated with *P. auricularis* HN038; T2: inoculated with *B. aryabhattai* JX285; T3: inoculated with the mixture of HN038 and JX285 (1:1); CK1: inoculated with LB medium; CK2: inoculated with sterile water. Different lowercase letters within the same column indicate significant differences (*p* ≤ 0.05). \*\*, significant effect at *p* ≤ 0.01.

#### *3.2. Gas Exchange Parameters*

The Pn was significantly (*p* ≤ 0.05) influenced by P treatment and PSB treatment, while the Tr and WUE were significantly (*p* ≤ 0.05) influenced by the interaction between P application and PSB inoculation (Table 2). At the same P application level, the Pn, Tr, and WUE were significantly (*p* ≤ 0.05) enhanced by the single and mixed inoculation of PSB strains (Table 2). The Pn, Tr, and WUE of plants co-inoculated with *B. aryabhattai* and *P. auricularis* were higher than that of other treatments. The inoculation of single or mixed PSB strains promoted Pn and Tr more effectively at a high P level (P3 treatment), while the PSB's best improvement effect on WUE appeared at the intermediate level (P2 treatment).

**P Application Levels PSB Inoculation Pn [**μ**mol**/**(m2**·**s)] Tr [**μ**mol**/**(m2**·**s)] WUE (**μ**mol**/**mmol) SPAD** P1 T1 3.41 ± 0.07 c 2.47 ± 0.07 a 1.38 ± 0.03 b 48.40 ± 1.64c T2 3.64 ± 0.05 b 2.59 ± 0.06 a 1.40 ± 0.03 b 51.37 ± 0.64 b T3 4.48 ± 0.20 a 2.69 ± 0.09 a 1.66 ± 0.07 a 64.83 ± 2.70 a CK1 2.46 ± 0.25 d 1.58 ± 0.09 b 1.74 ± 0.10 a 43.63 ± 2.21 d CK2 2.79 ± 0.43 d 1.73 ± 0.03 b 1.62 ± 0.23 a 44.83 ± 1.72 d P2 T1 4.14 ± 0.35 c 2.20 ± 0.09 ab 1.88 ± 0.08 a 62.93 ± 2.40 b T2 4.73 ± 0.40 b 2.66 ± 0.07 a 1.78 ± 0.14 a 68.87 ± 1.66 a T3 5.69 ± 0.36 a 2.80 ± 0.04 a 2.03 ± 0.21 a 71.30 ± 1.64 a CK1 2.93 ± 0.45 d 1.87 ± 0.04 b 1.56 ± 0.21 b 47.83 ± 1.50 c CK2 3.07 ± 0.36 d 1.83 ± 0.04 b 1.68 ± 0.16 b 48.31 ± 1.14 c P3 T1 3.76 ± 0.10 c 2.37 ± 0.03 ab 1.58 ± 0.02 a 52.73 ± 2.01 c T2 3.93 ± 0.29 b 2.42 ± 0.05 a 1.62 ± 0.10 a 57.87 ± 0.40 b T3 4.82 ± 0.20 a 2.87 ± 0.06 a 1.68 ± 0.04 a 66.17 ± 2.15 a CK1 2.89 ± 0.23 d 1.79 ± 0.07 b 1.61 ± 0.10 a 46.23 ± 1.00 d CK2 3.12 ± 0.16 d 1.90 ± 0.10 b 1.59 ± 0.02 a 47.36 ± 0.92 d Two-way ANOVA results P 0.00 \*\* 0.03 \* 0.00 \*\* 0.00 \*\* PSB 0.01 \* 0.00 \*\* 0.01 \*\* 0.00 \*\* P × PSB 0.17 NS 0.00 \*\* 0.17 NS 0.00 \*\*

**Table 2.** Effect of PSB on the photosynthetic characteristics of *C. oleifera* leaves.

Data are means ± SD (*n* = 3). Pn: net photosynthetic rate, Tr: transpiration rate, WUE: water use efficiency, SPAD: relative chlorophyll content. Different lowercase letters within the same column indicate significant differences (*p* ≤ 0.05). \*, significant effect at 0.01 < *p* ≤ 0.05; \*\*, significant effect at *p* ≤ 0.01; NS, no significant effect.

#### *3.3. Chlorophyll Content*

The chlorophyll content was significantly (*p* ≤ 0.01) influenced by the interaction between P application and PSB inoculation (Table 2). The single and mixed inoculations of PSB strains could significantly (*p* ≤ 0.05) increase the chlorophyll content of leaves regardless of the amount of P applied in the medium (Table 2). The co-inoculated plants had a higher chlorophyll content than that of either plants inoculated with single PSB strains or non-inoculated plants. The best promoted effect of the PSB strains on chlorophyll content was found at the intermediate P level. No significant difference in chlorophyll content was observed between the CK1 and CK2 treatments.

#### *3.4. Nitrogen and Phosphorus Content in Leaves*

The N content of the leaves was significantly (*p* ≤ 0.05) influenced by P treatment and PSB treatment, while the P content of the leaves was significantly (*p* ≤ 0.01) influenced by the interaction between P treatment and PSB treatment (Table 3). The trend of the changes in the N content of leaves was similar to that of the P content of leaves. The N and P content of the leaves increased to varying degrees by the single and mixed inoculations of PSB strains under different P application levels (Table 3). At the low (P1) and intermediate (P2) P levels, the N and P content of leaves were significantly increased by the PSB inoculation. The mixed inoculation of PSB strains increased the leaf N and P content more effectively than the single inoculation. However, there were no significant differences in the N and P content with the inoculation of PSB strains at a high P level (P3).


**Table 3.** Effect of PSB on the nitrogen content of *C. oleifera* leaves.

Data are means ± SD (*n* = 3). Different lowercase letters within the same column indicate significant differences (*p* ≤ 0.05). \*, significant effect at 0.01 < *p* ≤ 0.05; \*\*, significant effect at *p* ≤ 0.01; NS, no significant effect.

#### *3.5. Nitrogen, Phosphorus, and Potassium Content of Rhizosphere Soil*

The soil available N content was significantly (*p* ≤ 0.01) influenced by P treatment and PSB treatment (Table 4). A significant increase in the available N content was found in the treatment of the inoculated PSB strains compared to the non-inoculated control (*p* ≤ 0.05) (Figure 1). The maximum available N in the soil was found in the co-inoculated treatment under the intermediate P level. Moreover, no significant difference was observed in the available N content between the CK1 and CK2 treatments.


**Table 4.** Results of a two-way ANOVA of the effects of phosphate solubilizing bacteria (PSB) inoculation,

phosphorus (P) application and their interactions on the soil parameters.

\*\*, significant effect at *p* ≤ 0.01; NS, no significant effect.

**Figure 1.** Effects of PSB on the available nitrogen content of *C. oleifera*. Note: Data are means ± SD (*n* = 3). P1: no phosphate fertilizer (0 g/kg); P2: added 3 g calcium superphosphate (1 g/kg); P3: added 6 g calcium superphosphate (2 g/kg); T1: inoculated with *P. auricularis* HN038; T2: inoculated with *B. aryabhattai* JX285; T3: inoculated with the mixture of HN038 and JX285 (1:1); CK1: inoculated with LB medium; CK2: inoculated with sterile water. Different lowercase letters within the same column indicate significant differences (*p* ≤ 0.05).

The results showed that the soil available P content was significantly (*p* ≤ 0.01) influenced by the interaction between P treatment and PSB treatment (Table 4). Different inoculation treatments caused different degrees of increase in the available P content in the soil (Figure 2). Under the intermediate and high P levels, the inoculation of single and mixed PSB strains significantly increased the content of available P in the soil (*p* ≤ 0.05). Conversely, there were no significant differences in the available P content among the inoculation treatments under a low P level (*p* > 0.05).

**Figure 2.** Effects of PSB on available phosphorus content of *C. oleifera*. Note: Data are means ± SD (*n* = 3). Different lowercase letters within the same column indicate significant differences (*p* ≤ 0.05).

The soil available K content was significantly (*p* ≤ 0.01) influenced by the interaction between P treatment and PSB treatment and (Table 4). An inoculation of PSB strains increased the available K content in the soil (Figure 3). The co-inoculated treatment had a higher available K content than single and non-inoculation treatments under each P level.

**Figure 3.** Effects of PSB on the available potassium content of *C. oleifera*. Note: Data are means ± SD (*n* = 3). Different lowercase letters within the same column indicate significant differences (*p* ≤ 0.05).

#### **4. Discussion**

Phosphorus deficiency in soil is an important limiting factor in global agro-forestry production [25,26]. Some microorganisms can increase the concentration of available P by secreting organic acids and various degrading enzymes (phytase, nuclease, phosphatase, etc.) to decompose insoluble phosphate in the soil [27,28]. However, these microorganisms release little P in their natural state [11]. Therefore, it is necessary to artificially inoculate some high-efficiency P solubilizing microorganisms to increase the amount of available P. In our previous experiments, we isolated native PSB strains from the rhizosphere soil of *C. oleifera* and screened two strains (JX285 and HN38) with a high phosphorus solubilizing ability [17,18]. The results of the present study showed that the single and mixed inoculations of JX285 and HN38 had positive effects on the growth, photosynthetic capacity, N and P content in the leaves of *C. oleifera*, and available N, P, and K content in the soil.

Plant growth is the most obvious characteristic for evaluating the effects of PSB [29]. In this study, the single and mixed inoculations of two native PSB strains increased the plant height and biomass of *C. oleifera*, supporting the findings reported by previous studies [10,30]. This might be due to the PSB strains of JX285 and HN38 dissolving the insoluble phosphate in the soil and enhancing the available P content by producing organic acid and extracellular phosphatases [27,28]. Another possibility might be related to the metabolism of PSB, producing a variety of plant hormones, acids, and vitamins [30].

Plants use the light energy absorbed by chlorophyll molecules to drive photosynthesis [31]. This study showed the beneficial effect of PSB inoculation on the chlorophyll content and photosynthetic capacity of *C. oleifera* that is consistent with a previous observation in rice [32]. The higher Tr and WUE found in plants inoculated with PSB strains indicated that PSB improved the water status of *C. oleifera* [33]. An enhancement of Tr suggested an increased water uptake capacity, providing additional water for transpiration and improving soluble nutrient uptake [34].

Phosphorus is an essential nutrient, being a component of vital molecules in plants, and is involved in many metabolic processes [28]. PSB may convert insoluble P compounds into soluble forms by the processes of chelation, acidification, and exchange reactions [30]. In the current study, the increased available P content in the rhizosphere soil by inoculating PSB strains was only found in the intermediate and high insoluble P levels, while no significant difference was observed between the

inoculated and non-inoculated plants under a low P level. This may result from the uptake of P by *C. oleifera* for growth under a low P level, as supported by the improved leaf P content and the growth performance of PSB under a low P level. Under a high P level, the available P released by PSB was sufficient for the plants' growth, which led to the PSB with no effect on the leaf P concentration. PSB not only solubilize and mineralize P from insoluble compounds but also release other nutrients [1,35]. In this study, the available N and K content of rhizosphere soil and the leaf N content were promoted by the PSB inoculation, demonstrating that PSB elevated the amounts of available N, P, and K in the soil and subsequently provided better nutrition for plant growth [36,37].

This study also found that the inoculation effect of mixed PSB strains was better than that of single strains. According to previous studies, there may be a synergistic effect between different microorganisms [38,39]. On one hand, the combination of different PSB strains may have a higher and more stable cell activity than a single strain [40]. On the other hand, the combination could better exploit the limited P sources in soil [30]. However, the mechanisms of a synergistic interaction remain to be explored.

#### **5. Conclusions**

In this study, we probed the effects of single and mixed inoculations of two native PSB strains, JX285 and HN38, on the growth of *C. oleifera* seedlings. The results showed the positive influence of JX285 and HN38 on plant growth, photosynthetic capacity, the N and P content of the leaves, and the available N, P, and K content of rhizosphere soil. The two PSB strains in the co-inoculation acted synergistically with each other and strengthened the beneficial effects on plant growth performance. The use of PSB as inoculants may provide an alternative to chemical fertilizers and promote sustainable agroforestry.

**Author Contributions:** L.Z. conceived the idea and established the direction; F.W. wrote the first draft of the manuscript; J.L. (Jianrong Li), Y.C., S.W., Y.Z., L.L. and X.S. carried out the experiment. S.W., J.L. (Junsheng Liang), and Y.C. analyzed the data. All authors contributed with suggestions and corrections, and approved the final manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (31660189, 31570594), and Hunan Provincial Natural Science Foundation of China (2018JJ2217, 2018JJ3281).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Forests* Editorial Office E-mail: forests@mdpi.com www.mdpi.com/journal/forests

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18