Spatial and seasonal variations of leaf area index ( LAI ) in subtropical secondary forests related to floristic composition and stand characters

1 Faculty of Life Science and Technology, Central South University of Forestry and Technology, Changsha 410004, Hunan Province, China 2 Huitong National Field Station for Scientific Observation and Research of Chinese Fir Plantation Ecosystem in Hunan Province, Huitong 438107, China 3 National Engineering Laboratory of Applied Technology for Forestry & Ecology in Southern China, Changsha 410004, China 4 Changsha Environmental Protection College, Changsha 410004, China 5 Institute of Environment Sciences, Department of Biological Sciences, University of Quebec at Montreal, Montreal, QCH3C 3P8, Canada


Introduction
Many fundamental ecological processes in forest ecosystems, such as carbon (C) flux as well as water and energy exchanges, take place between the canopy layer and atmosphere (GCOS, 2006;Brut et al., 2009;Alonzo et al., 2015;Liu et al., 2015b).At a finer scale, leaves within the canopy are the primary organ to perform a series of physiological activities (i.e.photosynthesis, respiration, and evapotranspi-Published by Copernicus Publications on behalf of the European Geosciences Union.ration) (Aragão et al., 2005) and physical reactions (i.e.rainfall and radiation interception) (Aston, 1979;Smith, 1981;Crockford and Richardson, 2000).Therefore, the amount of leaves in a forest is the determinant of above-ground ecological processes and ecosystem functions.Leaf area index (LAI), defined as total one-sided leaf area per unit ground surface area (Biudes et al., 2014), is a widely used parameter (Kross et al., 2015) to quantitatively describe the vegetation canopy structure (Woodgate et al., 2015), to simulate ecological process models (Brooks et al., 2006;Sprintsin et al., 2007;Facchi et al., 2010;Gonsamo and Chen, 2014), and to reveal tree growth and productivity in forests at stand scale and landscape level (Lee et al., 2004;Liu et al., 2015b).In addition, LAI is listed as one of the essential variables for observation of global climate (Mason et al., 2003;Manninen et al., 2009) and for remote sensing data validation (Asner et al., 2003;Clark et al., 2008).Thus, accurate estimates of LAI values are important to understand ecological processes in forest ecosystems.
At present, various direct and indirect methods have been developed to measure LAI in forests.Direct estimation methods including leaf harvest (Clark et al., 2008), allometric equations, and litter collection (Ryu et al., 2010;Liu et al., 2015a) are recognised as the most accurate.However, leaf harvest and allometric equations methods need timeconsuming, labour-intensive, and destructive sampling processes, while litter collection is more feasible for temperate deciduous forests.Obviously, the direct methods are less applicable to large-scale and long-term LAI monitoring (Bequet et al., 2012;Biudes et al., 2014).Indirect methods include using a plant canopy analyser (Licor LAI-2000), hemispherical or fisheye photography (Macfarlane et al., 2007), and remote sensing (Biudes et al., 2014).The indirect methods retrieve LAI value from light transmittance through canopies or from canopy image analysis.For large-scale LAI estimates, remote sensing is the most effective method but requires validation with ground-based LAI data.LAI estimates on the ground at small scales are still a challenge due to the problems of sampling strategies associated with accepted level of accuracy, time, and cost considerations (Richardson et al., 2009).Hemispherical photography is a relatively simple and easily operated method among many indirect methods to retrieve LAI value at small scales (Demarez et al., 2008).Correction of the effects of woody materials, clumping and zenith angels or exposure is critical to improve the accuracy of LAI estimation (Liu et al., 2015b).Analysis software development and portable and timely characteristics allow hemispherical photography to measure spatial heterogeneity and seasonal variations of LAI in forests.
Forest canopy structure is highly complex so LAI values show great temporal and spatial variations at scales ranging from stand to global scale.For example, LAI values in the 7.9 ha plot of an old humid temperate forest tended to increase spatially as elevation increased and showed a temporal variation with plant phenology (Naithani et al., 2013).
The spatial patterns of LAI values at stand scale were significantly influenced by spatial distribution of tree species, which was dependent on topography and soil types (Naithani et al., 2013).The coefficient of variation (CV) in LAI decreased as the scale increased and LAI values did not have any relationship to biome type and climate patterns, but they were influenced by land use and land cover, terrain features, and soil properties at stand scale (Aragão et al., 2005).The CV of LAI of three species (i.e.beech, oak, and pine) had different degrees of spatial variation in a 1 ha plot at stand level (Bequet et al., 2012).LAI values in sagebrush displayed strong spatial patterns with time after disturbance and increased with stand age and total plant cover (Ewers and Pendall, 2007).The LAI values derived from MODIS data (Myneni et al., 2002;Huang et al., 2008) revealed strong spatial variations at global scale, which were correlated with latitude (Tian et al., 2004).At the global scale, temperature is the limiting factor for LAI under cool conditions while water plays a predominant role under other conditions, and this pattern differed among plant functional types (Iio et al., 2014).The factors that govern the spatial variations in LAI values at stand level include forest types, stand structure (Bequet et al., 2012), climate (Shao and Zeng, 2011), topography, soil moisture condition (Breshears and Barnes, 1999), and human disturbance and management activities (Huang and Ji, 2010).Although effects of topography, soil properties (Aragão et al., 2005;Naithani et al., 2013), and stand characters (Bequet et al., 2012;Yao et al., 2015) on LAI values have been investigated in detail, the effect of forest type, stand structural diversity, and stand structure on spatial heterogeneity and seasonal variations of LAI has yet to be fully understood.
Chinese subtropical forests contain a diversity of tree species with complex canopy structure that mostly grow on heterogeneous topography and soil conditions.As a result, LAI in subtropical forests may exhibit great spatial and seasonal variations, which is worthy of further investigation.However, LAI data of subtropical forests are relatively deficient in the global database (see Asner et al., 2003).In this study, we selected three different forests: Pinus massoniana-Lithocarpus glaber coniferous and evergreen broadleaved mixed forests, Choerospondias axillaris deciduous broadleaved forests, and L. glaber-Cyclobalanopsis glauca evergreen broadleaved forests, which were measured by using hemispherical photography to measure LAI values.Spatial heterogeneity of LAI was investigated through geostatistical analysis, and generalised additive models (GAMs) were used to examine how stand structural diversity and stand characters affect LAI variations in the three forests.Specifically, the objectives of this study were (1) to examine differences and seasonal variations in LAI among three forests in subtropical China; (2) to analyse spatial heterogeneity of LAI values within a specific forest; and (3) to identify how forest types, stand structural diversity, and stand characters control the spatial heterogeneity and seasonal variations of LAI values in three forests.

Study site description
The study was carried out at Dashanchong Forest Farm (latitude 28  , 2006).Evergreen broadleaved forest is the climax vegetation of the region.As a result of human disturbance and management activities, the farm has no primary forest, and possesses a range of secondary forests in different stages of succession (based on species composition) dominated by different tree species, including (1) early-stage P. massoniana-L.glaber coniferous and evergreen broadleaved mixed forests dominated by the shade-intolerant coniferous species typical of early succession, (2) middle-stage C. axillaris deciduous broadleaved forests dominated by shadeintolerant deciduous broadleaf species, and (3) late-stage L. glaber-C.glauca evergreen broadleaved forests dominated by the shade-tolerant evergreen broadleaved species commonly observed in the late stage of succession in this farm (Xiang et al., 2015;Ouyang et al., 2016).

Determination of stand characteristics
We established a permanent plot for each of three forests (i.e. 90 m × 190 m irregular plot for P. massoniana-L.glaber mixed forests, 100 m × 100 m plot for C. axillaris deciduous forests, and 100 m × 100 m plot for L. glaber-C.glauca evergreen broadleaved forests).Each plot was divided into 10 m × 10 m subplots, where tree species, diameter at breast height (DBH, in centimetres), tree height (H , in metres), height under the lowest live branch (in metres) and crown width (in metres) were measured for the individual stem with DBH larger than 1 cm.Stand characteristics for the trees with DBH > 4 cm of the three forests are presented in Table S1 in the Supplement.
To identify the factors that control spatial heterogeneity of LAI values in the forests, we selected individual trees with H larger than average height of each stand (see Table S1) and calculated their stem number, average DBH, H , total basal area at breast height (BA), crown width, crown coverage (calculated from crown diameter measured for individual trees within a stand), tree species diversity, tree size diversity, the proportion of BA of three functional groups (coniferous, deciduous, and evergreen broadleaved species) to total stand BA within a subplot.Tree species diversity (biodiversity index, BDI) was determined using the Shannon-Wiener index as follows: where P i is important value of ith species and is calculated by dividing the sum of relative abundance degree and relative dominance degree of ith species within a subplot by two.
Based on the Shannon-Wiener index, 2 cm was used for the DBH class, so tree size diversity (H ) was determined using the formula of Lei et al. (2009): where P i is the proportion of basal area for the ith diameter class.

Sampling design for LAI measurement
At the centre of each subplot of the three forests, hemispherical photographs were taken using a LAI measuring instrument (SY-S01A, Shiya Scientific and Technical Cooperation, Hebei, China) throughout four measurement seasons, i.e. in April (spring), July (summer) and October (autumn) in 2014, and January (winter) in 2015.The operation was carried out below canopy with the fisheye lens (Pentax TS2V114E, Japan) 1.0 m above the ground (Manninen et al., 2009) with a viewing angle of 180 • .The picture exposure is automatic exposure set by the manufacturer, and we took the photographs (768 × 494 pixels, BMP) in the morning, at dusk, or when cloudy in order to minimise influence of direct sunshine (Rich, 1990;Bequet et al., 2012).The images were processed and effective LAI values (L e ) were recorded using plant canopy analysis software developed by the manufacturer, for which appropriate pixel classification (thresholding) was chosen (752(H ) × 494(V ), where V is vertical resolution), viewing angle was considered (150 • ), and the hemispherical photography was divided into five rings to obtain results.To obtain accurate LAI (L), the correction was made to L e based on previous theory (Chen, 1996): where α is the ratio of woody to total area and reflects the contribution of woody materials to L e , and E is the clumping index that quantifies the effect of foliage clumping beyond shoots level.In the method getting accurate E values, the hemispherical photography was divided into 10 sectors.γ E is the needle-to-shoot-area ratio and quantifies the effect of foliage clumping within shoots.Photoshop software (Adobe Photoshop CS5, Adobe Systems Incorporated, North America) was used to calculate α.Photoshop software, we used the clone stamp tool to select the image of the woody materials (e.g.stems) and excluded the pixels, leaving only leaves on the photos, recorded as LAI of leaves (LAI leaf ).The value of α was calculated accordingly: The logarithm averaging method proposed by Lang and Xiang (1986) was applied to calculate E : where P (θ ) is the average gap fraction (expressed without the bar in the text), ln[P (θ )] is the logarithm average of the gap fraction, and P k (θ ) is the gap fraction of segment k.For deciduous and evergreen broadleaved species, γ E = 1.0; for coniferous species γ E is always > 1.0, but we ignored the effect of needle to shoot area on LAI in this study.

Data analysis
The minimum, maximum, mean value, standard deviation, and CV were calculated for the LAI data measured in 100 plots within each forest.Two-way analysis of variance (ANOVA) was used to detect effect of forest type and measurement season on LAI value.The LAI data in the three forests were tested for normal distribution using the K-S test (P < 0.05).We followed Chiang et al. (2003) in regarding LAI values as normal when they fell within the mean value ±3 standard deviations.Otherwise, the LAI values were regarded as outliers and replaced with the maximum or the minimum of normal values.Because the geostatistical analysis requires that the data meet normal distribution, the transformation was applied when the data did not meet normal distribution (Dai et al., 2014).Most values required natural logarithm transformation to meet assumptions of normality.The exception is for L. glaber-C.glauca in April and in November, which were ARTAN-transformed.
To investigate spatial heterogeneity of LAI values over four seasons measured in the three forests, the semivariance function was calculated as follows: where γ (h) is the semivariance value of lag distance h, N (h) is the number of pair data for lag distance h, and Z(x i ) and Z(x i + h) represent LAI values at coordinate x i and (x i + h) (Rossi et al., 1992).Based on the semivariogram plotting γ (h) values against hvariable, the appropriate models were fitted and we obtained the values of nugget (C 0 ), sill (C 0 +C), range (A 0 ) (Ewers and Pendall, 2007), and the ratio [C / (C 0 + C)] that reflected the degree of spatial autocorrelation of LAI values in a forest.Because spatial autocorrelation and semivariogram theory make unbiased optimal estimation for regional variables in a limited area (Bivand et al., 2013), the Kriging interpolation method, an unbiased estimation of the regional variables of the sampling points using the structure of the data and semivariogram function, was used to predict unknown LAI values in the forests from the data measured and to produce spatial distribution maps of LAI values for the three forests and four seasons.Compared with other methods, the Kriging method can overcome the difficulty in analysing error of interpolation, does not produce the boundary effect of regression analysis, and estimates the spatial variability distribution of measured parameters.Ordinary Kriging -one of the Kriging methods -is a least-squares method of spatial prediction based on the assumption of an unknown mean.It is the most common type of Kriging in practice (Dai et al., 2014) and is widely used in soil spatial heterogeneity studies (Elbasiouny et al., 2014).In our study, we also used the ordinary Kriging interpolation method to investigate spatial heterogeneity of LAI values.
Because the largest amount of defoliated leaves occurs in January and leaves fully expand in July in subtropical forests, we chose LAI values measured in January and July in three forests as response variables.The explanatory variables include forest types, stand structural diversity (species richness, tree species diversity, and tree size diversity) and stand characters (stem number, average DBH, H , BA, crown width, crown coverage, and the proportion of two functional groups (deciduous and evergreen conifer species) to total stand BA).GAMs are able to analyse complex and nonlinear relationships (Guisan et al., 2002;Austin, 2002;Wood, 2006).Therefore, we used GAMs to examine how the factors affect LAI values.The function of GAMs is the addition of many smooth functions and each smooth function has an explanatory variable.In our study, we chose smooth spline with two splines as the smooth method for GAMs.The variance inflation factor (VIF) -the ratio of the regression coefficient variance for a variable when fit with all variables to that for the variable if fit on its own -was used to test the multi-collinearity of explanatory variables (James et al., 2013).When the VIF of an explanatory variable is between 0 and 10, the variable was retained to the model; otherwise, we discarded the variable (Shen et al., 2015).The Akaike information criterion (AIC) or generalised cross-validation (GCV) was used to determine whether the model was good or bad (Clark, 2013).The factors selected after the multicollinearity test were used for multi-factor analysis.After all the possible models in multi-factor analysis, we determined the optimal model based on the significant influence of all explanatory variables in the model with the smallest AIC or GCV (Dong et al., 2012).Geostatistical analysis was performed with GS+ software (Gamma Design Software).Statistical analysis and GAMs analysis were operated in R 3.2.1 (R Development Core Team, 2015).The car packages were used to test multi-collinearity and the GAM packages were used to select the optimal model.3 Results

Variation in LAI values
The LAI values varied with forest type and measurement season (Table 1).Generally, LAI differed significantly between measurement seasons (P < 0.001), but LAI difference was not significant among forest types (P > 0.05).Interactive effects of measurement seasons and forest types on LAI were significant (P < 0.01).Among three forests, LAI in the P .massoniana-L.glaber forest had relatively low variation, while LAI in the L. glaber-C.glauca forest had the highest variation.In the P .massoniana-L.glaber forest, LAI showed the largest variation (the highest CVs) in October and the lowest variation (the smallest CVs) in January.In the C. axillaris forest, the largest variation in LAI was found in April and the lowest was found in January.In the L. glaber-C.glauca forest, LAI showed the largest variation in April and had the lowest variation in July.
Mean LAI values in the three forests showed different seasonal variation patterns (Fig. 1).The C. axillaris forest exhibited a unimodal pattern of seasonal variation, with the maximum mean LAI value (3.11 ± 1.18) occurring in July and the minimum mean LAI value (1.28 ± 0.44) in January.In the P. massoniana-L.glaber forest and L. glaber-C.glauca forest, the maximum mean LAI values occurred in October and the minimum mean LAI values appeared in January.During the growing season (April and July), the C. axillaris forest had the highest mean LAI value and the L. glaber-C.glauca forest had the lowest mean LAI value.During the non-growing season (October and January), the L. glaber-C.glauca forest had the highest mean LAI value in January, while the P. massoniana-L.glaber forest had the highest mean LAI value in October, and the C. axillaris forest had the lowest mean LAI values.
Mean α values in the three forests showed different seasonal variation patterns (Table 2).The C. axillaris forest exhibited a unimodal pattern of seasonal variations in mean α value, with the maximum mean α value occurring in January and the minimum mean α value in July.No obvious seasonal variations were found for the mean α value in the P. massoniana-L.glaber forest and in the L. glaber-C.glauca forest.Mean E values in the three forests were between 0.84 and 0.92, but they did not show clear seasonal variations, and the standard deviations were small.

Spatial heterogeneity in LAI values
The semivariogram results for LAI across the three forests during different measurement seasons are summarised in Table 3.The spatially dependent variance [C] accounted for 88.9-98.4% of the total variance [C + C 0 ] for LAI values measured in January in the three forests and also in April, July, and October in the L. glaber-C.glauca forest.This indicated the strong spatial autocorrelations of LAI values over short distances.These LAI data were best fitted with a Gaussian model or exponential model (r 2 > 0.50).
Spatial autocorrelation ranges of LAI values differed among forests and measurement seasons (Table 3).In January, the largest spatial autocorrelation range was found in the P. massoniana-L.glaber forest, and the smallest was found in the C. axillaris forest.In April, the largest spatial autocorrelation range of LAI was found in the C. axillaris forest, and the smallest was found in the P. massoniana-L.glaber forest.In July, the largest spatial autocorrelation range of LAI was in the P. massoniana-L.glaber forest, while the smallest was in the C. axillaris forest.In October, the largest spatial autocorrelation range of LAI was in the L. glaber-C.glauca forest, while the smallest was in the P. massoniana-L.glaber forest.Seasonal changes of range showed one peak pattern for C. axillaris forest and L. glaber-C.glauca forest, where the large range appeared in the growing season (April and July) and the small range appeared in the non-growing season (October and January).
Spatial distribution pattern of LAI values also varied with forest type and measurement season (Fig. 2).For example, LAI values in January across the three forests exhibited obvious patch and heterogeneous spatial distribution.In April and July, less spatial heterogeneity was found for LAI values especially in the P. massoniana-L.glaber forest.In October, heterogeneous and patch spatial distributions of LAI values appeared in the

Factors affecting LAI variation
The multi-collinearity test indicated that the explanatory variables in January and July did not have multi-collinearity.Thus, forest type, species richness, tree species diversity, tree size diversity, stem number, average DBH, H , BA, crown width, crown coverage, and the proportion of two functional groups (deciduous and evergreen conifer species) to total stand BA were included as explanatory variables in multi-factor analysis for LAI values measured in January in the three forests.After comparing all possible models, the best-fitted GAMs for LAI values in January were expressed as LAI ∼ s (stem number, 2) + s (crown coverage, 2) + s (PESB, 2) + s (PDSB, 2) + factor (forest types) (Table 4).
For LAI values measured in July, all these factors selected by the multi-collinearity test were included as explanatory variables in multi-factor analysis.The best-fitted GAMs for LAI values in July were expressed as LAI ∼ s (stem number, 2) + s (PDSB, 2) (Table 4).
The explanatory variables included in GAMs reflected their effects on or relationship with LAI variations.Given that other variables were fixed, LAI measured in January tended to decrease as stem number increased.LAI showed a positive nonlinear relationship with crown coverage up to ∼ 200 m 2 and then decreased with increasing crown coverage.The LAI values tended to increase as the proportion of evergreen conifer species to total stand BA increased and tended to decrease as the proportion of deciduous species to www.biogeosciences.net/13/3819/2016/Biogeosciences, 13, 3819-3831, 2016 The significance of the regressions (P ) are * , * * , and * * * for P < 0.05, 0.01, and 0.001 respectively.
total stand BA increased (Fig. 3).Given that other variables were fixed, LAI measured in July tended to increase as stem number increased up to ∼ 7 and then decreased at higher values.The effect of the proportion of deciduous species to total stand BA on LAI appeared more complicated, in that LAI increased as the proportion of deciduous species to total stand BA increased up to ∼ 0.7 and then decreased at higher values (Fig. 4).

Discussion
4.1 Seasonal variation in LAI value among forest type LAI data in subtropical forests in southern China are lacking compared to other global regions (Asner et al., 2003).This study provided seasonal LAI data in three subtropical forests that consist of contrasting functional types of species.Their mean LAI values varied from 1.28 ± 0.44 to 3.28 ± 1.26 (Table 1).This result is close to the LAI range (from 1.0 in winter to 4.0 in summer) retrieved by remote sensing techniques from the subtropical area of China from 2000 to 2010 (Liu et al., 2012).Compared with the LAI values estimated from al-lometric equations (Xiang et al., 2016) and specific leaf area values in 40 m × 40 m plots in this study (5.29-9.19), the LAI values measured by hemispherical photography are low but significantly correlated (r 2 = 0.40 and P = 0.035).Previous studies (see Lopes et al., 2015) have proved the underestimation of LAI using hemispherical photography.However, the method is feasible to obtain forest LAI data and to investigate spatial and seasonal variation in such values (Coops et al., 2004;Dovey and Toit, 2006).
The ratio of woody to total area (α) and the clumping index ( E ) have been recognised as the error sources in LAI measurement by optical methods (Chen et al., 1997;Bréda, 2003;Liu et al., 2015a).So far these two parameters have been measured in northeastern China (Liu et al., 2015a, b), which showed that the α values ranged from 0.04 ± 0.01 to 0.69 ± 0.12 and E values ranged from 0.88 ± 0.04 to 0.96 ± 0.01.These values were measured in temperate forest in northeastern China and differed from our study (mean α values varied from 0.04 ± 0.03 to 0.15 ± 0.09 and mean E values varied from 0.84 ± 0.09 to 0.92 ± 0.08) (Table 2), so they are not suitable for LAI correction in subtropical forests.Also literature on α and E values in subtropical forests is scarce.The variations in α are probably due to the seasonal variations and spatial heterogeneity of canopy structure in the three forests.In general, the α values are consistent with the amount of leaf litter.Our results showed that the large mean αvalues occurred in autumn for the P. massoniana-L.glaber forest and the C. axillaris forest but in spring and autumn for the L. glaber-C.glauca forest (Table 2).This seasonal change in mean α value in three forests was generally consistent with the amount of leaf litter collected by a litter trap installed in each forest type (Guo et al., 2015).The average E value (0.87) in this study was smaller than the values of mixed broadleaved Korean pine forest in northeastern China (Liu et al. 2015b) and this could be attributed to the different region and forests.The values of α and E obtained in this study fill the gap of calibration for optical measurement of LAI in subtropical forests.
Mean LAI values differed among the three forests and the differences were significant between the C. axillaris forest and the other two forests at a given measurement season.The C. axillaris forest had a relatively high mean LAI value dur- ing the growing season but changed to the lowest mean LAI value during the non-growing season.The change in mean LAI values in the C. axillaris forest was consistent with the study of a deciduous species-dominated forest reported by Naithani et al. (2013).It has been reported that the forests consisting of different plant functional types showed different LAI values (Asner et al., 2003;Iio et al., 2014).The differences and seasonal variations of LAI values in the three forests could be attributed to floristic composition and phenological defoliation patterns of tree species especially the deciduous species.The C. axillaris forest consisted of 74.15 % deciduous species, 25.80 % evergreen broadleaved species, and 0.05 % evergreen coniferous species, while the proportions of deciduous species were 10.05 and 25.70 % in the P. massoniana-L.glaber and L. glaber-C.glauca forests respectively.Seasonal growth and defoliation of different functional types of species lead to the change in leaf lifespan and foliage area (Niinemets, 2010) during different seasons related to temperature and water availability, which are responsible for the unimodal pattern of seasonal variation in mean LAI values.This agrees with the results of Liu et al. (2012), where the highest LAI was found in summer (July), followed by autumn (October) and spring (April), and the lowest was found in winter (January).

Within-forest spatial heterogeneity and factors controlling LAI
Semivariograms of LAI values in the three forests were fitted with spherical, Gaussian, exponential, or linear models (Table 3).Based on the fitted models, the degree of spatial autocorrelation could be evaluated.Spatial autocorrelation is weak when the determination coefficient (r 2 ) of the bestfitted semivariogram model is less than 0.5 (Duffera et al., 2007).The ratio [C / (C 0 +C)] is also used to describe the degree of spatial autocorrelation.A ratio of between 0 and 0.25 indicates a weak spatial autocorrelation, of between 0.26 and 0.75 indicates moderate autocorrelation, and of more than 0.75 indicates strong autocorrelation (Lopez-Granados et al., 2004).Spatial autocorrelation of LAI in this study varied with forest and measurement season (Table 3).Strong spatial autocorrelation in LAI values at a short range measured in January in all three forests indicated the sampling distance is reasonable for LAI variables within the spatial range (Liu et al., 2008).On the contrary, weak autocorrelation indicated that more samples and smaller sampling intervals should be taken to determine spatial dependency of LAI, such as for LAI measured in April in the P. massoniana-L.glaber forest.Spatial heterogeneity in LAI values was different for forest type and measurement season.Our study described spatial variations in LAI value by CV and geostatistical analysis, and the results were largely consistent with each other.In general, the CVs of LAI values in the three forest types (in particular C. axillaris forest) were higher for the period of leaf onset (April) and senescence (October) than for the period of leaf maturity (July) (Table 1).This reflects changes in leaves due to plant phenology and is consistent with the study of Naithani (2013) where LAI became increasingly homogenous from leaf onset to maturity but became more heterogeneous from maturity to senescence.As a result, degree of heterogeneity in LAI value for all three forests tended to dwindle from leaf non-growing season to growing season (Fig. 2).
The complex hydrothermal environment results in complex vertical and horizontal variation in canopy layer and formed unique spatial heterogeneity in LAI values.The effects of stand characters on LAI have been examined and positive and negative effects have been reported (Tobin et al., 2006;Bequet et al., 2012;Yao et al., 2015).In our study, results from GAMs showed that forest types, stand structural diversity, and stand characters affected spatial heterogeneity of LAI values significantly in the three forests.This finding that floristic composition and stand characters affected LAI values measured in July is consistent with the study of Yao et al. (2015); LAI values increased with stem number but when the stem number was larger than 7, LAI values decreased with stem number mainly due to the floristic composition in these study areas.Because July is the period of leaf maturity for deciduous species and leaves fully expand in this season, LAI values tended to increase as ratio of deciduous species increased, but when the ratio was higher than ∼ 0.7, its negative relationship with LAI probably could be explained by the strong competition among tree species, with diverse species composition and the canopy overlap among tree species (Fig. 4).Our results indicated that LAI values did not exhibit a significant relationship with stand BA, consistent with the findings of Mcdowell (2007); total LAI did not exhibit a clear pattern in relation to stand BA.
Until now, the non-growing season relationship of LAI variation with forest type and stand characters has been seldom reported.In this study, forest type, stem number, crown coverage, proportion of evergreen conifer species to total stand BA, and proportion of deciduous species to total stand BA and forest type were the factors significantly affecting LAI variation in January.As January is mainly the leaf senescence period of deciduous species, LAI values in January decreased with stem number and decreased with deciduous species ratio.The relationship between LAI value and the evergreen species ratio was generally the reverse of that between LAI and the deciduous species ratio.The fact that LAI values in January decreased with increasing crown coverage when crown coverage was larger than ∼ 200 m 2 could be explained by large crown coverage resulting in more defoliation (in particular for deciduous species) in the forest in January (Fig. 3).The proportion of deciduous species to total stand BA significantly affected LAI variations in January and July, and the relationship between LAI and the deciduous species proportion was reversed when the ratio was smaller than 0.7 in these two seasons, which is consistent with the growth law of deciduous species.Thus, deciduous species play an important role in LAI variations across seasons.Also, the seasons have a significant effect on LAI variation by affecting leaf growth.The partial effects of stem number and crown coverage on the LAI values observed in January showed that these smooth functions were large at both ends of the 95 % confidence interval.This was due to the small sample number in this range, and most were concentrated in the middle parts, the same as the partial effects of stem number on the LAI values observed in January (Figs. 3,4).
Although the factors selected by regression could explain a small proportion (4 %) of spatial heterogeneity of LAI measured in July, the factors selected in January could explain 35 % of the LAI spatial heterogeneity (Table 4).The LAI heterogeneity also could be affected by several other factors, such as the topography (Naithani et al., 2012), soil feature (Choler et al., 2010), soil temperature (Vitasse et al., 2009;Hardwick et al., 2015), microclimate, human activity, and other physicochemical properties.However, full leaf expansion of all tree species, which covers up the effect of other physicochemical properties on LAI, leads to a small difference in LAI in July.The effects of environmental factors (e.g.temperate and rainfall) on LAI in the forests at the fine scale should be taken into account in future studies.
Spatial heterogeneity of LAI in the three forests can yield some useful information for sampling strategy to accurately estimate of LAI using indirect measurement.An optimal sampling strategy should consider appropriate sampling plot size and the lowest sampling number that, as far as possible, obtains a high sampling accuracy and a low sampling error (Bequet et al., 2012).Our study found that strong spatial autocorrelations range were ∼ 13-27 m (the minimal range was 13.80 m, and the maximal range was 27.00 m) (Table 3), indicating that the range from 13 to 27 m might serve as the reference for sampling plot size to estimate LAI in subtropical forests.In addition, LAI heterogeneity was closely related to floristic composition and stand characters; thus stand structural variables (BA or DBH) are important for sampling strategy to measure LAI in forests (Bequet et al., 2012).

Conclusions
This study measured LAI in three subtropical forests using a hemispherical photography method over four seasons, and offered reliable data to analyse spatial and seasonal variations in LAI.Our results indicated that LAI differed greatly with forest type and measurement season.Seasonal variation in LAI across the three forests reflects defoliation due to plant phenology.LAI values for all three forests exhibited different spatial autocorrelation in the four seasons.A clear patch distribution pattern in LAI value was found during the nongrowing seasons and this pattern gradually dwindled in the growing seasons.While stem number, crown coverage, proportion of evergreen conifer species to total stand BA, the proportion of deciduous species to total stand BA, and forest type significantly affected spatial variations in LAI values in January, stem number and proportion of deciduous species to total stand BA significantly affected spatial variations in LAI values in July.These findings supplement LAI data for global synthesis and will provide valuable information for sampling strategies to enable more accurate estimates of LAI for simulated models of production and hydrological cycles in subtropical forests.

Figure 3 .
Figure 3. Partial effects of stem number, crown coverage (m 2 ), the proportion of evergreen conifer species to total stand BA (PESB), the proportion of deciduous species to total stand BA (PDSB), and forest types (calculated for overstorey trees with height larger than average stand height) on the LAI values observed in January in P. massoniana-L.glaber, C. axillaris, and L. glaber-C.glauca forests.

Figure 4 .
Figure 4. Partial effects of stem number and the proportion of deciduous species to total stand BA (PDSB) (calculated for overstorey trees with height larger than average stand height) on the LAI values observed in July in P. massoniana-L.glaber, C. axillaris, and L. glaber-C.glauca forests.

Table 1 .
L. glaber-C.glauca forest, and banded spatial distributions of LAI values obviously appeared in the C. axillaris forest.Descriptive statistical characteristics of LAI values measured from April 2014 to January 2015 in P .massoniana-L.glaber, C. axillaris, and L. glaber-C.glauca forests (n = 100).

Table 2 .
Average woody to total leaf ration (α) and clumping index ( E ) values in P .massoniana-L.glaber, C. axillaris, and L. glaber-C.glauca forests.Values in parenthesis are the standard deviation of α and E values (n = 100).

Table 4 .
Estimated coefficients of the generalised additive models (GAMs) for the factors with effects on LAI values measured in P .massoniana-L.glaber, C. axillaris, and L. glaber-C.glauca forests.