Geographic variation in growth and wood traits of Neolamarckia cadamba in China

Forest trees are important because they provide not only timber but also raw materials for various purposes, including pulp, fuel, bioenergy feedstock, plywood, charcoal, medicines, and food. In addition to their usefulness, forest trees also contribute to terrestrial ecosystems. For example, they deliver important ecological services such as water conservation, climate regulation, biodiversity protection, and carbon recycling, and they have helped to shape the diversity of human culture^[1]. As the human population grows exponentially, the demand for wood is increasing sharply, leading to a worldwide acceleration in deforestation. Estimates suggest that about 2 billion hectares of forests must be restored worldwide to meet these needs^[2], and forest ecosystems are in urgent need of protection for sustainable development. Many countries have implemented natural forest protection projects, and these are very important for protecting biodiversity and improving the ecology of many regions. Such projects can have a remarkable effect on improving our environment and alleviating possible ecological disasters; for example, the implementation of strict forest protection policy in China led to an increase in vegetation that accounted for more than 25% of the total global increase from 2000 to 2017^[3]. However, such projects will inevitably intensify the conflict between timber demand and ecological/environmental considerations. Such a challenging problem can be addressed, to some degree, by conventional tree breeding, which can maximize the natural production potential of forest trees through genetic selection. However, the growth traits and wood properties of forest trees are polygenetic or complex traits regulated by multiple genes, which may also interact with the environment^[4]. Climatic factors such as precipitation, temperature and light levels all play important roles in the growth of tree species. Larsen et al. found that precipitation was significantly positively correlated with the annual ring-width growth of white spruce (Picea glauca) and jack pine (Pinus banksiana)^[5]. Nagamitsu et al. found that provenances of Japanese larch (Larix kaempferi) from eastern Japan, where seasonal variation in temperature is lowest and that of precipitation is highest, showed the highest performance in stem and branch growth^[6]. To improve the yield and quality of forest products and optimize the economic and ecological benefits of forest trees, we can use strategies for the selection of germplasm resources, such as provenance tests, which enable gene pyramiding and selection of the best provenances harboring genes that positively regulate trait(s) of interest. This also allows us to take advantage of geographic genetic variation^[7] and match the most appropriate trees to specific sites.

Neolamarckia cadamba (Roxb.) Bosser (Rubiaceae)^[8] is an evergreen broadleaved tree species found from South Asia to Australia. It has been introduced to Surinam, South Africa, Costa Rica, and other tropical and subtropical areas^[9] because of its tremendous economic and ecological value^[10]. N. cadamba is used not only for furniture-making but also for pulp and paper production^[11]. Its leaves, bark, flowers, and fruit are widely used in traditional and modern medicine^[12]. The species has also attracted significant attention for other reasons. For example, its pollen is used for apiculture^[13], its leaves for silage^[14], and the whole tree for its value in the landscape^[15]. To date, growth traits and wood properties have been the focus of the N. cadamba breeding program^[16−19]. At the World Forestry Congress in 1972, N. cadamba was described as 'a miraculous tree'. Under normal conditions, it can attain a height of 17.67 m and a trunk diameter of 25.3 cm at breast height within 9 years^[20]. Therefore, it has potential, in suitable regions, to meet the increasing demand for wood products.

Progeny tests of N. cadamba from Southern China have been conducted in Guangdong. Preliminary results from the field trials indicate that there are significant differences between provenances and families within provenances with respect to height and diameter at breast height, and these characteristics exhibited moderate heritability (0.53–0.79)^[21,22]. However, to date, the geographic variation among genetic breeding resources remains unexplored, and corresponding wood properties have not yet been reported. In the present study, we quantified differences in the growth and wood quality of different provenances and characterized associated patterns of geographic variation. In addition, we identified provenances that exhibited rapid growth and desirable wood properties, and we collected materials to preserve as germplasm resources and for use in future breeding activities.

Ten geographic provenances were used in this study, covering the entire natural distribution of N. cadamba in China (Fig. 1). The variation in all growth traits is shown in Table 1. For all traits, the degree of variation between provenances was significant (multiple comparison analysis at the P < 0.05 level). The provenances with the highest mean diameter at breast height (DBH), height (H), and tree volume (V) were GXFCG (13.31 cm), GXNN (13.09 m), and GXNN (0.1022 m³), respectively. Those with the lowest mean values of these traits were GDGZ (11.41 cm), GDGZ (10.91 m), and YNBS (0.0695 m³). The largest individual values of DBH, H, and V were all recorded in provenance GXFCG (20.90 cm, 19.5 m, and 0.3143 m³, respectively), and the smallest individual values of these traits were recorded in YNDH (3.9 cm), GXLZ (1.00 m), and YNDH (0.0022 m³). The ranges and coefficients of variation (CVs) are also presented in Table 1. The provenances with the greatest DBH, H, and V ranges were GXLZ (4.10–20.80 cm), GXLZ (1.00–19.00 m), and GXFCG (0.0127–0.3134 m³), respectively. Of all growth traits, V exhibited the highest ratio of the maximum to minimum value, reaching 24.68 fold. The YNJH provenance had the smallest range of DBH (6.50–18.20 cm), H (7.00–18.00 m), and V (0.0109–0.2128 m³), and the YNJH provenance had the lowest CV for all growth traits (21.94%, 23.64%, and 58.79% for DBH, H, and V, respectively), indicating low variation in growth traits for this provenance. The provenances with the highest CVs for DBH, H, and V were YNDH (31.54%), GXLZ (32.22%), and YNMN (81.63%), respectively.

Figure 1. 

Geographic origins of the 10 N. cadamba provenances (black dots) and location of the test site (red dot).

The mean coefficients of variation for DBH, H, and V were 27.79%, 29.05%, and 72.33%, respectively, indicating substantial potential for genetic improvement of growth traits among provenances. Out of all traits, DBH had the smallest CV, and V had the largest, far greater than those of the other traits, indicating a high potential for genetic improvement in N. cadamba.

Phenotypic variation in wood properties

Table 2 shows descriptive statistics for the wood properties of all provenances. Each trait exhibited a different degree of variation, and all traits except Cr differed significantly among the provenances. The GXLZ provenance had the largest mean values of FL, FD, FL/FD, VL, and VD, whereas the GDYF provenance had the smallest mean values of FL, FD, FL/FD, and VD. The largest and smallest mean values of WBD were for YNBS (0.3600 g/cm³) and YNDH (0.3253 g/cm³), respectively. The YNBS provenance had the largest range of FL and FD, and the corresponding CVs were also the largest. The largest and smallest ranges of other traits were from different provenances.

The mean CVs for FL, FD, FL/FD, VL, VD, VL/VD, WBD, and Cr were 8.50%, 3.46%, 6.63%, 13.17%, 15.55%, 12.59%, 8.15%, and 4.75%, respectively. Among these traits, the largest and smallest CVs were for VD and FD, respectively. The mean CVs for vessel traits were all greater than 10%, indicating that the potential for genetic improvement among provenances was greater for these traits than for the other wood properties investigated.

Differences in growth and wood property traits among provenances

The variance components, provenance heritabilities (h²), and genetic variation coefficients (CV_G) among the provenances were estimated (Table 3). The provenance variance components (V_P) for DBH, H, and V accounted for 4.97%, 3.87%, and 4.33% of the corresponding total variance, respectively. The provenance by block interaction variance components (V_PB) for DBH, H, and V accounted for 0.84%, 4.37%, and 1.47% of the corresponding total variance, respectively. The V_PB percentages for DBH and V were lower than the corresponding V_P percentages, whereas the opposite was true for H. This suggests that interactions between provenance and block were lower for DBH and V than for H. The variance components of provenance were always lower than the variance components associated with random error (V_e), suggesting that random environmental effects could be a major cause of variation in the studied traits, whereas the genetic effects were more limited. The provenance heritabilities of DBH, H, and V were 0.67, 0.59, and 0.64, respectively. This indicates that growth traits of N. cadamba are under moderate to high genetic control. By contrast, the provenance heritabilities for wood properties were in the range of 0.02–0.45, suggesting that these properties were under low to moderate genetic control. The differences in variance for FL, FD, FL/FD, VL, VD, and Cr among the provenances were not significant (the standard error was much larger than the corresponding estimated value).

For traits with a significant V_P, the genetic variation coefficient was in the range 0.59%–11.83%. Only the CV_G of V was greater than 10%, indicating that the provenance selection potential of V is greater than that of the other traits.

Analysis of phenotypic and genetic correlations

The genetic and phenotypic correlations between all studied traits are presented in Table 4 (the genetic correlations for FD and VD were not estimated because the corresponding provenance variances were zero). The genetic and phenotypic correlations between growth traits (DBH, H, and V) were large, positive, and significant (0.97–0.99 for genetic correlations, 0.87–0.98 for phenotypic correlations), and the genetic correlations were always larger than the corresponding phenotypic correlations. Low to moderate genetic correlations between wood properties were observed, except between FL and FL/FD. In the analysis of phenotypic correlations between wood properties, the correlations between Cr and other wood properties were always low (absolute values of the correlation coefficients were in the range 0.02–0.18). The correlations between WBD and other wood properties were also low (absolute values 0.02–0.28), except for those with FD (−0.52) and VL (0.99). There were significant moderate to high positive correlations between fiber traits and vessel traits (0.41–0.99), except for the correlation between VL/VD and FL (0.28) and the correlation between VL/VD and FL/FD (0.20). The genetic correlations between growth traits and wood properties were always low (absolute values 0.02–0.22) and negative (except for those between H and FL and between H and FL/FD), and the corresponding standard errors were always larger than the estimated correlation coefficients (except for those of correlation coefficients between WBD and growth traits). The phenotypic correlations between growth traits and wood properties were low to moderate (absolute values 0.07–0.41). The phenotypic correlations between growth traits and fiber/vessel traits were significant and positive (except for those between growth traits and VL/VD). The phenotypic correlations between growth traits and WBD or Cr were weak and negative.

Geographic patterns

To better understand possible trends associated with geographic patterns, we performed a binary quadratic trend surface analysis (Table 5). The regression equations for DBH, H, V, VL/VD, and WBD were significant, whereas those for FL, FD, FL/FD, VL, VD, and Cr were not, and subsequent discussions therefore focus on the former traits.

The trend surface diagrams (Fig. 2) revealed that the patterns of geographic variation in growth traits were basically opposite to those in wood properties. All three growth traits (DBH, H, and V) displayed the same pattern of geographic variation: a gradual increase from the periphery to the central region. Thus, as latitude and longitude increased, the growth traits first increased then decreased. Variation was lower in the central region than at the periphery, and latitude had a greater impact on growth traits than longitude. Unlike the growth traits, WBD showed a gradual decrease from the north and south to the central region; as latitude increased, WBD first decreased then increased, and as longitude increased, WBD exhibited a slight increase. The trend for VL/VD was also different, showing a gradual decrease from the southeast and northwest to the central region. Geographic variation in VL/VD was affected by latitude more than by longitude.

Figure 2. 

Power contour-trend surfaces for (a) diameter breast height (DBH), (b) height (H), (c) volume (V), (d) wood basic density (WBD), and (e) vessel length/vessel diameter (VL/VD) of N. cadamba. The surface represents geographic variation, and the lines on the surface represent contours.

Assessment of provenance performance

The main goal of forest tree breeding is to select superior provenances with rapid growth and high wood density; however, forest trees intended for different purposes should be assessed on the basis of different indicators. In the present study, we used V and WBD as indicators of wood product value and V and FL as indicators of good pulpwood sources. Taking the overall means of the relevant traits as thresholds for the selection of superior provenances for wood products, we identified five provenances for V (mean 0.0844 m³) and four for WBD (mean 0.3390 g/cm³). However, only one provenance—YNMS—had values that surpassed both thresholds (Fig. 3a). The realized gains (G) of this provenance for all studied traits are shown in Table 6. The largest realized gain was for V (12.50%), and the smallest was for VD (−3.15%). The realized gains for growth traits and VL/VD, WBD, and Cr were positive, whereas those for other fiber traits and vessel traits were negative.

Figure 3. 

Selection of superior provenances for wood products (a) and pulpwood (b).

For pulpwood, selection yielded five provenances with mean V values greater than the overall V mean (0.0844 m³) and five provenances with mean FL values greater than the overall FL mean (1,422.97 μm). Three provenances—GXLZ, GXFCG, and GXNN—had values that surpassed both thresholds (Fig. 3b). The realized gains (G) of these provenances for all studied traits are shown in Table 6. The largest realized gain was in V (25.00%), and other growth traits also had relatively high gains (5.11% and 5.93% for DBH and H, respectively), whereas the smallest gain was in VL/VD (−5.12%). The traits that are beneficial for pulpwood had moderate realized gains (3.47% and 2.73% for FL and FL/FD, respectively).

Yield per unit volume is a key indicator in forestry production, and individual tree volume is the most important indicator for predicting plantation yield potential. In the present study, all growth traits differed significantly between provenances, and the coefficient of variation for V was significantly greater (p < 0.001) than those for DBH and H, with the genetic variation coefficient exhibiting the same pattern. Similar results have been found in other studies that examined provenances of N. cadamba^[23] and other forest tree species such as Betula alnoides^[24], black spruce^[25], and Chinese Fir^[26]. These results suggest that, in most forest tree species, the variation in tree volume between provenances is greater than that of other growth traits, thus indicating that volume has greater potential for selection in provenance trials. The provenance heritability of different growth traits was in the range 0.59–0.67, similar to the results of other provenance experiments involving N. cadamba^[22,23]. This result suggests that the growth traits of N. cadamba are under moderate to high genetic control.

As one of the most important indicators of the physical properties of wood, wood density has a great influence on the final yield and quality of forest trees. We measured the wood properties of N. cadamba, including basic density, fiber properties, and vessel properties. The range of WBD was 0.26–0.41 g/cm³, similar to that found by Krisnawati et al.^[27] in Malaysia, who reported wood density (with 15% moisture content) in the range 0.29–0.56 g/cm³. Another study from Indonesia^[28] produced similar results (0.259–0.606 g/cm³). Fiber length and the ratio of fiber length to fiber diameter are the two most important indices for pulp quality^[29]. In the present work, FL and FL/FD were in the ranges 1,072.86–1,692.28 μm and 36.98–50.99, respectively. Rahayu et al.^[28] reported similar results with respect to fiber length (800–1,700 μm). There were some differences in wood properties (except for Cr) between provenances, but only those for WBD and VL/VD were significant. Wood properties not only vary between tree species, provenances, and families but also radially within a single plant. At the same time, they tend to be correlated with the growth rate of DBH, although this occurs mainly in fast-growing tree species. The wood properties of Scots pine (Pinus sylvestris) differ significantly between families^[30], the wood density of poplars (Populus spp.) differs significantly between clones^[31], and the wood density of Larix kaempferi differs significantly from the pith to the bark^[32]. In N. cadamba, radial growth rate can also affect wood properties^[28]. In the present study, however, the data we used were average values for provenances, and this may be why the provenance variance of some wood traits was not significant.

Genetic and phenotypic correlations between traits

In this study, the genetic and phenotypic correlations between growth traits were always strong and positive, as reported in other tree species^[26,31,33]. Such results may provide a theoretical basis for the indirect selection of growth traits. In previous literature, the correlation between growth traits and wood quality in conifers has always been found to be negative^[34]. For example, the genetic correlation between diameter and density was −0.61 for Picea abies^[35], −0.48 for Pinus radiata^[36], and −0.34 for Chinese fir^[26]. However, in broadleaved tree species, some researchers have found the genetic correlation between growth traits and wood quality to be weak. For example, the genetic correlation between tree height and density of Eucalyptus grandis declined towards zero with age^[37], and the genetic correlation between diameter and density of Populus tomentosa was −0.05^[38]. Similarly, the genetic correlations between growth traits and wood properties were weak in our study, suggesting that there is weak correlation between the genes controlling growth traits and wood quality. This result indicates that growth traits and wood traits can be selected independently when identifying superior provenances of N. cadamba.

Patterns of geographic variation in N. cadamba provenances from southern China

The geographic variation associated with the provenance of forest trees is one of the main components of variation at different levels within a tree species. To make use of this variation, provenance experiments can be used to study the effects of geographic variation on growth traits, wood quality, stress resistance, and so forth, thus identifying high-yielding, high-quality, and stable germplasm resources for forest production. Several studies from Southeast Asia have shown that there is significant variation among provenances of N. cadamba^[23,39,40]. In the present study, significant geographic variation was found in DBH, H, V, VL/VD, and WBD. However, the patterns of geographic variation associated with growth traits were the opposite of those for wood quality. This result may be related to the negative correlation between them, although this correlation was weak. Latitude had a greater influence on the traits than longitude, and the growth traits first increased and then decreased with increasing latitude. Similar results were reported in a study of Chinese fir^[41]. This pattern is closely related to the annual average temperature and precipitation in the original distribution area of the trees. Compared with the periphery of the natural range of N. cadamba, the central region has a relatively high annual average temperature and a more suitable precipitation regime. Because N. cadamba is very sensitive to low temperatures, the frostless period and minimum temperature are also major causes of geographic variation, and the central region has a relatively longer frostless period and higher minimum temperature than the periphery. On the other hand, Diaphania glauculalis is an important pest that attacks N. cadamba, and its damage period occurs mainly during the hot and humid season^[42]; such damage may influence growth indirectly. This may also be one of the main reasons for the pattern of geographic variation in N. cadamba growth traits. Thus, in addition to target traits such as growth and wood traits, resistance traits such as cold resistance and insect resistance should also be considered in the selection of superior N. cadamba provenances.

Provenance selection for wood products and pulpwood

Our analysis showed that the provenances with high growth traits did not necessarily have the best wood properties. Moreover, there was a negative genetic correlation between growth traits and wood properties, although this correlation was weak. This finding highlights the need for careful consideration when selecting superior provenances, rather than simply looking at the performance of a single trait. Taking account of the different intended uses, we selected WBD and V as selection indicators for wood products and FL and V for pulpwood. YNMS was selected as an excellent provenance for the cultivation of wood products. Compared with the overall average, V had the highest realized gain (12.50%), followed by DBH (3.91%). Among the wood properties, WBD, Cr, and VL/VD also showed moderate realized gain (2.94%, 2.47%, and 2.93%, respectively), whereas the realized gains for other wood properties were negative. Provenances GXLZ, GXFCG, and GXNN were selected as being superior for pulpwood: the realized gain for V was 25.00%, and other growth traits also had relatively high gains (5.11% and 5.93% for DBH and H, respectively). Traits that are beneficial for pulpwood production all showed moderate realized gains (3.47% and 2.73% for FL and FL/FD, respectively). Similar results were found in selection of Chinese fir^[26]; when constraining selection to DBH and WBD, growth traits and WBD increased dramatically, and hygroscopicity and tracheid diameter exhibited negative realized gains. In a breeding program for Pinus kesiya^[43], the gain in volume was found to be −9.78% when selecting only for the single trait of wood density. However, when selection was constrained to both DBH and WBD, there were simultaneous increases in strength and volume. Similar research results were also found in breeding programs for broadleaved tree species. In poplar^[44], when wood density or volume was selected alone, the genetic gain of another character was zero or even negative. In Eucalyptus dunnii, when a single trait was used as the selection criterion for pulp production^[45], the elite clones selected were different each time. In the long term, the best selection strategy may be to develop a multiple-trait selection index for each application. The selection of superior provenances will inevitably lead to a narrowing of genetic diversity^[6]. Therefore, it is necessary to continuously supplement new superior provenances in the breeding population.

Our results revealed that growth traits and most wood properties of N. cadamba differed among provenances in southern China, and the growth traits were under moderate genetic control. There were strong positive correlations between growth traits, whereas the correlations between growth traits and wood properties were weak. The patterns of geographic variation differed between growth traits and wood properties, although latitude had the greatest influence on both. The YNMS provenance was suitable for selecting superior timber-wood trees, and the GXLZ, GXFCG, and GXNN provenances showed promise for pulpwood trees. However, in the long term, a multiple-trait selection index for N. cadamba should be considered for different applications, including the production of medicine, juice, nectar, and silage, as well as the species' value as a landscape tree. Multi-site provenance experiments are also required to further refine our understanding of tree–site interactions and better match provenances with sites.

The data were collected from a half rotation-aged progeny trial in Leizhou (21°10′06″ N, 110°21′34″ E), Guangdong province. The site has an annual mean temperature of 22 °C and annual rainfall of 1711.6 mm. The minimum and maximum temperature in this region are 15.5 and 28.4 °C, respectively. The experimental plantation was established in the spring of 2014. The experimental design in the field comprised randomized complete blocks, with 10 blocks and 5-tree plots in a 3 m × 3 m square spacing. Ten geographic provenances were planted in the trial, covering the entire natural distribution of N. cadamba in Southern China (Table 7). The sampled trees were chosen because they were considered phenotypically average or above average with respect to stem DBH and total height compared with neighboring trees in the population. The distance between mother trees within the population was a minimum of 100 m to reduce genetic relatedness between seed lots. Seeds were collected by climbing the trees. The seed lots from each tree were kept separate, and their location and number of samples were recorded. Dr. Mingxuan Zheng undertook the formal identification of the plant material used in our study, and the plant specimens are housed in the herbarium of South China Agricultural University (CANT32205).

Data collection

Diameter at breast height (DBH in cm, 1.3 m above ground level) and height (H in m) of all trees were recorded. Individual tree volume (V in m³) was calculated using the following formula^[20]:

Five average trees were selected from each provenance based on DBH. From each, a 5.02-mm core was extracted at breast height using a tree growth cone. Wood basic density (WBD in g·cm⁻³) was determined using the water displacement method^[46] based on two measurements: volume of water displaced by immersion of the wedge (w₁) and oven-dry weight (w₂). Wood basic density was then calculated using the formula:

Vessel length (VL), vessel diameter (VD), the ratio of VL to VD (VL/VD), fiber length (FL), fiber diameter (FD), and the ratio of FL to FD (FL/FD) were determined following Chen and Xie^[47]. Each sample was placed in a test tube, and chromic acid-nitric acid separation solution equivalent to 10–20 times the sample volume was added; the mixture was then boiled over an alcohol lamp for 5–8 min. A small sample of the material was removed and placed on a glass slide, then pressed lightly with tweezers to determine whether it disintegrated. If not, heating was continued for another 2–3 minutes. After rinsing thoroughly with clean water, a small amount of 0.5% safranin aqueous solution was added, and the material was mashed with a glass rod. Samples were then examined under an optical microscope. Each sample was measured three times, and a total of > 30 values were obtained.

The degree of crystallinity (Cr) was measured by X-ray diffraction following Segal et al.^[48]. The diffraction intensity of wood fiber is at its maximum value at 2θ = 22°, and its integrated intensity is defined as I_u. The wave trough appears near 2θ = 18°, which is the scattering intensity of the diffraction in the amorphous area of the wood fiber, and its integrated intensity is defined as I_a. Crystallinity is then calculated using the formula:

Statistical analysis

R version 4.0.2^[49] and Microsoft Excel 2013 were used to analyze variation in growth traits and wood properties, including the mean value, standard error, amplitude, and CV. To determine the differences between phenotypic variables of the provenances, Duncan’s multiple range tests were performed using the agricolae package^[50] in R version 4.0.2. Variance and covariance components for genetic analyses were estimated using the sommer package^[51] in R version 4.0.2 based on a mixed linear model:

where y is a vector of trait phenotypes, β is a vector of fixed effects (block), u is a vector of random effects (provenance, provenance by block interaction), and ε is the residual. X and Z are incidence matrices for fixed and random effects, respectively.

To estimate the degree of genetic control for each trait, provenance heritability (h²) was calculated for all traits in the provenances overall using the following formula, based on the variance component estimates from the model analyses:

where V_P is the provenance variance, V_PB is the provenance by block interaction variance, V_e is the random error variance, n_h is the adjusted number of trees, and b is the number of blocks.

Because some of the poorly adapted saplings died, the number of trees in each block was not constant, and the data were unbalanced. Thus it was necessary to replace the original number of trees with the actual number of trees in each plot^[52]:

Where p is the number of provenances and n_ij represents the number of individual trees of the j^th provenance within the i^th block.

The genetic variation coefficient (CV_G) was calculated using the following formula^[53]:

Where $ \overline X$ is the average phenotypic mean of the trait, and V_P is the provenance variance.

The realized gain (G) was estimated by:

Where $ {\overline X_i} $ and $ \overline X $ are the mean values of the selected superior provenance and the overall trait, respectively.

The genetic and phenotypic correlations between traits were calculated as follows:

Where r_g and r_P are genetic and phenotypic correlations, respectively, and Cov_a1,a2 is the provenance covariance between traits a₁ and a₂. Cov_e1,e2 is the error covariance between traits a₁ and a₂, V_a1 and V_a2 are the provenance variances for trait a₁ and a₂, and V_e1 and V_e2 are residual error variances for traits a₁ and a₂.

Surfer 13.0 (Golden Software, Golden, CO, USA) software was used for trend surface mapping. The regression equation for the trend surface analysis was as follows:

Where β is the regression coefficient, x is the latitude, y is the longitude, and ε_ij is the random error.