Association distribution characteristics of natural and artificial Larix gmelinii forests along environmental gradients in northeast China

Background: Larix gmelinii forest is one of coniferous forests in cold-temperate zone, which is a vital part of national strategic landscape of ecological security of China. Plant association distribution is different in natural and artificial Larix gmelinii forests, meanwhile, determining mechanisms in typical associations of Larix gmelinii forests still need to be explored. The study focused on which environmental factors actuated association distribution of natural and artificial Larix gmelinii forests in northeast China. Two-way indicator species analysis (TWINSPAN) and canonical correspondence analysis (CCA) were used to classify plant associations and explored the relationship between species and environment. Results: All the plots (n=175, size=30 m×30 m) of Larix gmelinii forests were classified into 6 plant associations by TWINSPAN. Species diversity in natural forests were higher than that in artificial forests. Ass.III and Ass.IV only appeared in natural forests, meanwhile, Ass.VI only appeared in artificial forests. The primary environmental drivers of species diversity patterns in natural forests were annual mean temperature, followed by annual precipitation, elevation, slope aspect, and canopy density. However, elevation and annual precipitation had strong effects in determining association distribution in artificial forests. Conclusions: Plant association distribution showed habitat preferences, besides, natural forests had higher species diversity and more resistance than artificial forests. The study can be used as a reference for Larix gmelinii forest ecosystem protection in northeast China and a theoretical basis for scientific management in similar areas.


Background
Exploring plant community assembly has long been a central issue for understanding species coexistence and biodiversity maintenance in community ecology [1]. At biogeographical scales, patterns of spatial species turnover provided important insights into historical and regional processes underlying the composition and dynamics of regional biotas [2]. At local scales, spatial community dissimilarity was driven by several non-mutually exclusive processes, including dispersal limitation, habitat filtering and species interactions [3]. Dispersal, biotic interactions, and gap dynamics were likely to produce spatial structure most evident at relatively fine scales depending on underlying geomorphology [4]. Environmental heterogeneity was inherently spatially auto-correlated, and thus 3 the use of spatially variables helped to account for its variation, and also the spatial scale at which they influenced β-diversity [5]. Comparisons among ecological or taxonomic groups are necessary to gain a more comprehensive understanding of diversity patterns [6].
Larix gmelinii forests with a deciduous needle leaf conifer are adapted to grow in harsh climates and widely distribute over a range from 40°N to 72°N and 110°E to 130°E [7]. They are dominant timber tree species in northern China. To meet the needs of the timber production, most of the larch plantations are developed by replacing the secondary forests. The total area of larch plantations reaches 2.61 × 106 ha [8]. However, the soil fertility of larch plantations declined gradually due to single species composition and mono ecosystem structure [9]. Understanding community composition difference is the foundation for sustainable development of Larix gmelinii forests with the main purpose of timber production and ecological function [10]. In tropical forests, it was found that nearly four-fifth of the species examined were associated with topographically defined habitats [11]. Recent research focused on quantitative analysis to quantify the relative importance of environmental drivers on community composition and structure [12]. How environmental variation limits plant association distribution is of considerable current interest.
In the study, climate, topography and canopy density were considered as environmental variables.
TWINSPAN and CCA were used to analyze the relationship between species distribution and environmental variables. The study was designed to answer the following questions: (1) What is the present situation in terms of species composition and structure for natural and artificial Larix gmelinii forests? (2) Which environmental factors affected different association distribution patterns in natural and artificial Larix gmelinii forests? Based on the questions, we expected to find out the limiting environmental factors and provide important guidance for managing and protecting Larix gmelinii forests.

Study area
The study was conducted from July to September each year from 2016 to 2018. Based on typical and comprehensive sampling methods, a total of 175 Larix gmelinii forest plots (130 natural Larix gmelinii 4 forest plots and 45 artificial Larix gmelinii forest plots) were investigated in the Khingan Mountains (118.83°E -135.09°E, 38.72°N -53.56°N) (Fig.1). The region belongs to a cold temperate continental monsoon climate zone with a mean annual temperature and precipitation of -4℃-11.5℃ and 400-1100 mm, respectively. The elevation ranges from 200 to 1400 m a.s.l. The annual sunshine duration is 2300-3000 h and frost-free period is 100-160 d. The zonal vegetation types in the area are forest, shrub, grassland, wetland meadow, marsh, and aquatic vegetation. The most common tree species are Larix gmelinii, Pinus pumila, Betula platyphylla, Populus davidiana. Pinus pumila shrubs, Larix gmelinii forest and wetland distribute vertically from high to low along altitude. The soil type is brown coniferous forest soil, meadow soil, boggy soil [13].

Data collection
The plots basically covered the distribution range of Larix gmelinii forests in northeast China. The size of forest plots was 900 m 2 (30 m×30 m), and consisted of nine 10 m×10 m subplots [14]. Plant species in tree, shrub, and herbaceous layers were recorded in each plot. Species name, diameter at breast height (DBH), and height of all individual trees with a DBH 2.5 cm were investigated. Two of the nine subplots were selected for a survey of the shrub layer at the cross corners. Species name, abundance, and average height of all woody individuals with DBH 2.5 cm (including sapling trees) were recorded. Five 1 m×1 m herbaceous quadrats within each subplot were set up. Species name, average coverage, and average height of each vascular species occurring in the herbaceous quadrats were recorded.
Climate data that contained with 19 bioclimatic variables with 30s (ca.1 km) spatial resolution obtained was from the World Climate Database (www.worldclim.org). Detrended canonical correspondence analysis suggested that mean annual temperature (bio1) and annual precipitation (bio12) as climatic variables were more appropriate than others. Habitat factors were recorded, including altitude, slope, aspect, canopy density and disturbance. Ultimately, 6 environmental factors that mainly influenced species distribution of Larix gmelinii forests were added to model operation.

Data analysis
Species importance value (IV) was used as the basic data in multivariate analysis. The species data 5 matrix functions were consisted of importance values of 399 species in 130 natural Larix gmelinii forests plots (399×130) and 309 species in 45 artificial Larix gmelinii forests plots (309×45).
TWINSPAN is one of common methods in vegetation classification [15]. We named associations after the dominant species in the tree layer, followed by indicators species of shrub layers. CCA ordination can represent the variation of a data matrix in a reduced number of dimensions [16].  Table 2).

The response of different associations to environmental gradients
In the CCA analysis, the Monte Carlo permutation test indicated significant eigenvalues for all canonical axes (P 0.01). The species-environment correlations of CCA axes respectively were 0.944, 0.777, 0.806, 0.708. The cumulative percentage variance of the species-environment relationship for CCA axes was 83.5% in natural forests ( Table 3). The species-environment correlations of CCA axes respectively were 0.950, 0.922, 0.862, 0.861. The cumulative percentage variance of the speciesenvironment relationship for four CCA axes was 80.4% in artificial forests (Table 4). CCA ordination 6 performed well in plant association distribution along environmental gradients. In addition, natural forests had better fitting effect than artificial forests.
In two-dimensional CCA ordination diagram of natural forests, Asso.I and Asso.III usually distributed in higher elevation and slope, shady slope, lower annual mean temperature. Asso.II was immensely limited by annual precipitation, and mainly distributed in moist sunny slope with higher canopy density. Asso.IV usually distributed in cold shady slope with lower annual precipitation. Asso.V usually distributed in lower elevation and more flat areas (Fig.2). In artificial forests, Asso.I distribution trends was similar to the trends in natural forests. Asso.II and Asso.V were usually distributed in lower elevation with higher canopy density. Asso.VI was immensely limited by annual mean temperature and annual precipitation in lower elevation (Fig.3).

The main limiting environmental factors and species diversity patterns in different associations
In natural forests, environmental variable ranking of marginal effects in order was annual mean temperature, annual precipitation, elevation, slope aspect, canopy density, slope. Five environmental factors had significant marginal effects on community distribution (P < 0.05). Removing the collinearity of leading factors, only 5 environmental factors past through significance test. The conditional effect order of environmental factors didn't change (Table 5).
In artificial forests, environmental variables ranking of marginal effects in order was elevation, annual mean temperature, annual precipitation, slope aspect, canopy density, slope. Five environmental factors had significant marginal effects on community distribution (P < 0.05). Removing the collinearity of leading factors, only 2 environmental factors past through significance test. The conditional effect order of environmental factors changed, and ranking of environmental variables in order was elevation and annual precipitation ( Table 6).
Species richness ranged from 16.07 to 28.18, the Shannon-Wiener index ranged from 1.98 to 2.44, and evenness ranged from 0.69 to 0.81 in natural forests (Fig.4). Species richness ranged from 14.25 to 26.18, the Shannon-Wiener index ranged from 1.56 to and 2.46 evenness ranged from 0.58 to 0.71 in artificial forests (Fig.5). Ass.III and Ass.IV only appeared in natural forests and Ass.VI only appeared in artificial forests. Species diversity in Ass.I, Ass.II, Ass.V of natural forests were higher than that of artificial forests.

Discussion
Larix gmelinii forests are major components of temperature coniferous forests in northern China, and play an important role in climate change mitigation, water and soil conservation, and resource conservation [17]. Based on a field investigation of Larix gmelinii forests, we found that the composition of plant flora was rich, but significantly different among different mountain areas due to the complex topographic conditions and high spatial heterogeneity. In all plots, Compositae and Ranunculaceae of the angiosperms held a dominant position in the species community composition.
The classification demonstrated that Ass.III and Ass. IV only appeared in natural forests and Ass.VI only appeared in artificial forests. The results were in accordance with facts that Ledum palustre,Vaccinium vitis-idaea and Rhododendron dauricum usually appeared in stable community of Larix gmelinii forests [18]. However, Sorbaria sorbifolia and Geranium wilfordii were partial to distribute in lower altitude e and often appeared in artificial forests with strong disturbance [19].
Compared with the natural forests, some shrub layer was absent in the artificial forests. A tendency of increasing clustering for the woody plants along the temperature gradient, and indicated that environmental filtering drove the closely related species to occur in the same community [20]. The results of our study were consistent with other studies, which indicated environmental filtering might drive the community structure at a broad scale [21]. The Monte Carlo permutation test showed significant species-environment correlations with the CCA axes and natural Larix gmelinii forests showed better imitative effect.
Species diversity is an important feature of forest composition and structure, and its change is frequently used as an indicator of forest dynamics. Many researchers had pointed out that species richness was often strongly correlated with climate. However, most of the previous studies at a large scale were based on species richness within geographic grids [22]. In the analysis, the data at plot scale was used to examine geographic diversity patterns in relation to climate. It was found that species diversity in natural forests were higher than that in artificial forests due to natural 8 environment and association type. Higher canopy cover might cause a decrease in the redundancy of species with similar characteristics. The limited light conditions, microhabitat and interspecific competition played important roles in the patterns of species diversity for the artificial forests [23].
Stands at early successional stages were typically dominated by shade intolerant, nutrient demanding, and fast growing species, which resulted in phylogenetic and functional clustering, due to the availability of abundant resources subsequent to a disturbance [24]. Late-successional communities were more often dominated by shade tolerant, slowed growing and distantly related species with dissimilar functional traits, which were characterized by phylogenetic and functional over dispersion caused by competitive exclusion as resources become limited [25]. Such transformation has been generally attributed to shifts in community assembly processes, from environmental filtering to competitive exclusion [26]. Species diversity in the natural forests was more stable, however, higher coverage might limit light availability and enhance interspecific competition.
Therefore, plant community would like to regulate themselves to maintain appropriate biodiversity as succession.
Marginal effect and conditional effect provided ranking of environmental variables in order of importance in Larix gmelinii forests. The differences in forest type and structure caused the variations in light condition, soil property, interference intensity, which had strong effects in determining species composition and interspecific relationships [27]. CCA performed well in describing association distribution along environmental gradients. Due to the complex topographical conditions and high spatial heterogeneity in natural forests, the primary environmental drivers of species diversity patterns were mean annual temperature, followed by annual precipitation, elevation. Most studies have shown that the patterns and underlying mechanisms of woody plants and herbaceous plants were different due to the different response patterns [28]. As to natural forests, more environmental factors took part in determining species distribution, such as annual mean temperature, annual precipitation, elevation, slope aspect, canopy density. The results suggested that natural forests were relatively well developed and the system stability was affected by climatic, topographic, structural factors [29][30][31][32]. Nutrient resorption of plants can respond to the variation of soil nutrients by reducing the dependence of plants on soil nutrient supply to ensure their maximum growth [33]. Elevation and annual precipitation had strong effects in determining association distribution of artificial forests.
Overall, habitat preferences of association distribution in Larix gmelinii forests were obvious, besides, natural forests had higher species diversity and more resistance than artificial forests.

Conclusions
Plant association distribution showed habitat preferences in natural and artificial Larix gmelinii forests. Objective guidelines was established to reflect actually association spatial situation of Larix gmelinii forests by CCA. Natural forests were more resistant to disturbance than artificial forests. The results demonstrated that association distribution was affected by multiple environmental and ecological variables acting at different spatial scales overtime. The primary environmental drivers of species distribution in natural forests were annual mean temperature, followed by annual precipitation, elevation, slope aspect, canopy density. In addition, elevation and annual precipitation had strong effects in determining association distribution in artificial forests. Our results provided a theoretical basis for scientific conservation and restoration of Larix gmelinii forests in northeast China.