Associations between Road Density , Urban Forest Landscapes , and Structural-Taxonomic Attributes in Northeastern China : Decoupling and Implications

A better understanding on the associations between road density (RD), urban forest structural-taxonomic attributes, and landscape metrics is vital for forest ecological service evaluations and suitable management in sprawling urban areas with increasing road networks. We chose Harbin, a fast growing provincial capital city in northeast China, as a case study to address this issue. We utilized ArcGIS software (Esri, version 10.0; Redlands, CA, USA) and FRAGSTATS (V4.2.589) to digitize GF-1 images (Gaofen No.1 remote sensing images) to acquire road net characteristic information and landscape metrics of urban forests in Harbin. Together with forest structural-taxonomic attributes from a stratified random sampling survey, statistical methods such as an analysis of variance, a regression analysis, and a redundancy analysis were used to determine the road-dependent differences and to decouple the associations between them. The results indicated that road area percentages, road length/imperious surface area (ISA) ratios, road area/ISA ratios, and road cross-points sharply increased from low to heavy RD areas. This road intensification was strongly associated with increased urban forest area, patch density, and diverse patch shapes; smaller tree sizes, lower tree densities, and diverse tree species compositions were generally observed. Redundancy-based variation partitioning showed that part of the variations in structural-taxonomic attributes of forests could be explained by road intensity characteristics. In low RD (0–1.5 km/km2) regions, the road characteristics significantly affected forest characteristics (Shannon Wiener diversity index, species richness, and evenness index); however, such associations weakened with increasing forest landscape-related associations in medium to heavy RD (1.5–6 km/km2) regions. Our findings highlighted that road development is strongly associated with forest characteristics in Harbin city, and RD-dependent forest landscape regulating management could favor the maximization of forest ecological services that are related to structural and species identities.


Introduction
Road density (RD), commonly used to measure the impact of roads on landscapes [1,2], is strongly associated with the ecological risks for local ecosystems because road networks run through various landscapes [3].Areas with high RD are often characterized by large land areas under construction and those with low forest coverage [4], leading to a significant decline in landscape structure and ecosystem health [5].Linear transportation infrastructure (e.g., roads) and vehicles affect the structure of ecosystems and the dynamics of ecosystem functions.In addition, they have several direct or indirect ecological impacts on ecosystem components, including animal and plant species composition [6], original habitat fragmentation, changes in the physical-chemical environment and microclimates [1], invasion of weeds and pest animals, and elevated rates of poaching and wildfires [7].During the process of road development in urban regions, the number of original trees decrease and landscape fragmentation increases [1].To date, many studies have investigated the impact of road development on various landscape patterns [4,8,9], but research on the impact of road development on forest structural or taxonomic attributes is scarce.In China, road-affected areas were 18.37% of the total terrestrial area [10] and the government has invested 8.91 billion RMB into road greening practices.The total greening road mileage was over 50 thousand km in 2017 (http://www.ce.cn/xwzx/gnsz/ gdxw/201803/12/t20180312_28436372.shtml).The policy implementation of roadside afforestation has significantly increased the total forest area in China [11].These authority regulations on the greening of road networks might strongly interact with forest characteristics; nevertheless, there are very limited well-defined researches on their associations to date.
As an important part of urban green infrastructures, forests and trees have certain effects on the water, heat, carbon, and pollution cycles of a city [12][13][14][15].The fractal dimension index, patch numbers, and average patch area of forests may decrease with the development of the road network density [16].The most common conclusion regarding the influence of RD on landscape patterns is that RD could reduce the average original vegetation patch area and increase the number and distance between these patches [4,17].Moreover, there is a negative correlation between the plant species richness and the intensity of road construction [17].At the city level within the buildup region, a well-designed configuration of roads and urban forests could potentially maximize ecological services, and a basis for such configuration is an exact understanding of their complex associations [18].Landscape metrics derived from remote sensing classification results have been used to evaluate temporal and spatial changes of forests as well as management effectiveness [19][20][21][22].Together with remote sensing methods and field surveys, urbanization-induced variations of forest structural traits, such as tree height, diameter at breast height, canopy size, and compositional traits, and their association with carbon sequestration have been studied [13,15,23].Furthermore, it is reported that urban forest landscape metrics (e.g., forest aggregation and patch sizes) are strongly associated with structural-taxonomic attributes, and these metrics are possible indicators of urban forest characteristics [24].To date, no study has simultaneously studied road characteristics (both representation at urbanization levels and a specific ecological meaning), landscape metrics of forest, and forest characteristics with regards to species composition and individual sizes.However, the association and decoupling of these metrics will possibly favor urban greenspace management based on road ecological considerations [6,7,25].
Recent methodological advances dealing with complex associations in the natural environment gave hints for decoupling relationships between road development, urban forest landscapes, and structural-taxonomic attributes.The associations between road development and forest landscape pattern has been explored through correlation analysis and curve fitting [26], as well as association coefficient and the trend analysis [27].Principal component analysis (PCA) could quantitatively analyze the impact of road development on forest landscape patterns [17].Redundancy ordination (RDA) and variation partitioning have been used to decouple complex associations among various ecological factors [23,28,29], such as association between glomalin-related soil carbon sequestration and climate and soil physiochemical properties [23,30], tree microclimate regulation's association with background conditions [28], plant species diversity's associations with forest community features [29], and the association between carbon sequestration and forest landscape patterns [31].Co-utilization of regression analysis, RDA, and variation partitioning will favor the decoupling of the complex association among road development, forest landscape, and structural-taxonomic attributes.
Using the city of Harbin (the north-most provincial capital city in China, with a population of over 9.6 million in the administrative regions) as an example, in this study, we hypothesized that road development negatively affects forest structural-taxonomic attributes and forest landscape patterns (e.g., road development would reduce forest structural and taxonomic attributes and drive landscape fragmentation) and that the decoupling of their associations may favor urban greenspace management.The main questions addressed in this study are as follows: (1) What are the changing patterns of urban forest landscape metrics, structural-taxonomic attributes, and road characteristics at different RDs?
(2) How should the associations among RD, landscape metrics, and forest structural-taxonomic attributes be ordinated, and which landscape metrics and road characteristics are indicative of forest structural-taxonomic attributes?
(3) Are there any suggestions for the reasonable planning and management of urban forests that will optimize the ecological service benefits of urban forests?
To verify our hypothesis, we first obtained urban forest field-survey data and extracted and calculated road and landscape metrics data.Then, the road-dependent changes of the forest structural-taxonomic attributes and landscape patterns were studied by heavy, medium, and low road density classification.Finally, the associations and decoupling between road characteristics, forest structural-taxonomic attributes, and landscape patterns were explored by ordination and variation partitioning analysis for possible improvements in road-related landscape planning and management of urban forests.

Study Area
This study was conducted in Harbin (45 • 45 N, 126 • 38 E), the capital city of Heilongjiang Province in Northeastern China.Harbin is a provincial metropolis with the highest latitude and lowest temperature in China.The climate is typical of a temperate continental monsoon climate with an average annual rainfall of 569.1 mm and an average annual temperature of 5.2 • C [32].The surveyed region examined in this study was an area within the four-ring road area of Harbin, with an area of about 345.31 km 2 (Figure 1a).The historical natural vegetation in Harbin is mainly grassland with sparse elm trees (equivalent to savannahs in temperate regions).The other tree species are willow and poplar [33].The urban forest types we investigated included road forests, institutional forest, landscape forests, and farmland shelterbelts.Harbin has improved its green infrastructure and planted many trees.The urban forest types we investigated included road forests, institutional forests, landscape forests, and farmland shelterbelts.In a recent survey, the total forest area in the Harbin buildup area (within the 4th ring road) is 39.7 km 2 , a sharp increase from 6.19% to nearly 12% of the total land area [34].
The reform and opening-up policy is a national top-down economic and governing system reform in China that has been recognized as the only way for China to achieve modernization and become a developed country.Since the reform and opening-up policy began in 1978, the pace of construction in Harbin, the economic development, and the urban population growth increased yearly [35].In the past 30 years, the buildup area in Harbin increased by 65% (from 200 km 2 in the 1980s to 330 km 2 in 2010).During this period, Harbin's ring road was expanded from a two-ring road to a four-ring road [15].As the road development of Harbin represents a typical case in China, it has been used to study the effects of urbanization [23].Many studies in Harbin have focused on CO 2 dynamics [36], remote sensing tree measurements [37], habitat variability [38], mycorrhizal changes [23,39], soil fertility mapping [40], forest carbon sequestrations [31], the biodiversity of trees and birds [32,41], microclimate regulations from urban trees [14,42], and the tree census method development [43].The associations between road development and forest characteristics were not reported to date, and new studies will provide support for urban green infrastructure development.

Road Density Identification and Road-Related Traits Computation
The technical route of the urban road and forest landscape extraction road density grading is shown in Figure 1b.By referencing the Google Roadnet map, the road identification was based on GF-1 images (full color resolution 2 m and multispectral resolution 8 m) and extracted manually in the ArcGIS software (Esri, version 10.0; Redlands, CA, USA) (Figure 1b).We mainly extracted urban first-class roads, second-class roads, and parts of third-class roads within the four-ring road of Harbin (Figure 2a).RD was calculated as the length of roads divided by the terrestrial area [44].The RD fraction ranged from 0 to 6 km/km 2 (Figure 2a), and it was used as a proxy for the road network intensity.We split the study area into multiple 3 km × 3 km grid cells (a total of 80 grid cells), and the average RD value in each grid was assigned to one of the three levels of RD in each grid by using a partition statistics method of ArcGIS (Esri, version 10.0; Redlands, CA, USA) (Figure 2b).The three levels of RD were as follows: low RD region (RD < 1.5 km/km 2 ), medium RD region (1.5 km/km 2 ≤ RD ≤ 3 km/km 2 ), and heavy RD region (RD > 3 km/km 2 ).By using Google Roadnet map data as a reference, a precision check was carried out and we found that the remote sensing extraction data matched well with the reference data (R 2 = 0.72, p < 0.001) (Figure 2b).
We also calculated several road characteristics associated with RD, including road area density (RAD) and cross point density (CPD), which are equal to road area and the number of road cross points per unit terrestrial area (3 km × 3 km grid), respectively.GF-1 images were also used to digitize and acquire imperious surface area (ISA).ISA was acquired using the object-oriented feature extraction method in the ENVI 5.2 (Exelis Visual Information Solutions, USA) and by manual modification using the ArcGIS software (Esri, version 10.0; Redlands, CA, USA) to increase accuracy.The ISA data were subsequently used to compute road characteristic-related parameters, including road length/ISA (RLI) and road area/ISA (RAI), which are equal to the length and area of roads per unit ISA-impervious surface area within a 3 km × 3 km grid, respectively.

Urban Forest Characteristics Survey: Structural-Taxonomic Attributes
We conducted field investigations from August to October 2014 using a stratified random sampling method for obtaining the structural and compositional attributes of forests in Harbin city (Figure 1b).Each plot covered at least an area of 400 m 2 of urban forests.Regular shaped plots in this study was 20 m × 20 m.For irregular shaped plots, we adjusted the width and length to ensure an area coverage of at least 400 m 2 .A total of 194 plots, distributed among 61 grid cells, were field surveyed, including 23 grid cells in low RD regions, 19 grid cells in medium RD regions, and 19 grid cells in heavy RD regions.We identified every tree species in each sample point during the field survey.We measured six structural attributes, including tree height, crown size, diameter at breast height (DBH), under branch height (UBH), tree density (TD), and section area at breast height (DBH section area).Crown size was measured as the projected size for four directions.A mean value of each tree was used as the radius of the projection canopy size of the measured trees.DBH was measured at 1.3 m at the aboveground level, and measurements were recorded only if DBH was over 2 cm.
Three taxonomic attributes were measured in each plot, including woody plant species richness (SR), Shannon Wiener species diversity index (H'), and evenness index (J').The relevant formulae are given below [45][46][47] where s is the total number of species in each plot, RA is relative abundance, and RA i represents the number of trees of each species i/total number of trees [47].
All the field-surveyed data (1-13 plots in each grid) were averaged to single values for structural attributes (height, DBH, crown size, UBH, DBH section area, and TD) and taxonomic attributes (SR, H', and J') for future data analysis with a good match of road characteristics and landscape patterns.All the field-surveyed data (1-13 plots in each grid) were averaged to single values for structural attributes (height, DBH, crown size, UBH, DBH section area, and TD) and taxonomic attributes (SR, H', and J') for future data analysis with a good match of road characteristics and landscape patterns.

Analysis of Urban Forest Landscape Patterns
We used an object-oriented-based classification method to extract the urban forest of Harbin according to texture information, spectral information, and spatial attributes of remote sensing

Analysis of Urban Forest Landscape Patterns
We used an object-oriented-based classification method to extract the urban forest of Harbin according to texture information, spectral information, and spatial attributes of remote sensing images (Figure 1b).The attribute assignment and manual modifications were performed using the ArcGIS software (Esri, version 10.0; Redlands, CA, USA).An object-based validation method was used to evaluate the accuracy of the road extraction by referencing manually digitalized results from Google Earth satellite images (resolution, 0.59 m).The overall accuracy of the tree coverage extraction was >97%, with a miss factor of 0.06 and detection rate and quality percentage of >94%.The details can be obtained from our previous publication [31].We used the ArcGIS software (Esri, version 10.0; Redlands, CA, USA) to cut out the image of Harbin urban forest with a grid of the same size as the four-ring road area of Harbin, and then imported the forest image of each 3 km × 3 km grid cell into the FRAGSTATS (V4.2.589) software (University of Massachusetts, Amherst, MA, USA) to calculate the landscape metrics (a total of 80 grid cells).
Landscape metrics are often used for quantitative analysis of landscape patterns [22,24].In our study, the dispersed homogeneous forest region in each grid cell was considered a forest patch.Forest patch-related landscape metrics were then calculated as 10 frequently used metrics for quantifying forest landscape variations based on their ecological meanings and research purposes (Table 1).These 10 metrics include the following: two area-edge metrics (edge density (ED) and total area (TA)); three shape metrics (mean fractal dimension index (FRAC-MN), mean perimeter area ratio (PARA-MN), and area-weighted mean contiguity index (CONTIG-MN)); and five aggregation metrics (patch density (PD), mean Euclidian nearest neighbor distance (ENN-MN), interspersion and juxtaposition index (IJI), landscape shape index (LSI), and patch cohesion index (COHESION)).In this study, raster data and FRAGSTATS (V4.2.589) software version 4.1 were used to calculate the aforementioned metrics.
Table 1.List of landscape metrics and their descriptions.
TA equals the total area (m 2 ) of the landscape, divided by 10,000 (to convert to hectares).

Edge Density (ED)
ED = E A (10000) E = total length (m) of edge in landscape.ED ≥ 0; Units, meters per hectare.
ED equals the sum of the lengths (m) of all edge segments in the landscape, divided by the total landscape area (m 2 ), and multiplied by 10,000 (to convert to hectares).
PD equals the number of patches in the landscape, divided by total landscape area (m 2 ), and multiplied by 10,000 and 100 (to convert to 100 hectares).PD is one of the most important indices to describe landscape heterogeneity.
Landscape Shape Index (LSI) LSI equals 0.25 (adjustment for raster format) times the sum of the landscape boundary and all edge segments (m) within the landscape boundary and divided by the square root of the total landscape area (m 2 ).LSI measures the aggregation of patches.A large LSI value indicates a more irregular landscape.

Mean Fractal Dimension Index
FRAC-MN equals the mean value of the sum of 2 times the logarithm of the patch perimeter (m) divided by the logarithm of the patch area (m 2 ) for each patch of the corresponding landscape.It reflects shape complexity.FRAC approaches 1 for shapes with very simple perimeters such as squares and approaches 2 for shapes with highly convoluted, plane-filling perimeters.
Mean Perimeter Area Ratio (PARA-MN) PARA-MN equals the mean value of the ratio of the patch perimeter (m) to area (m 2 ), and PARA is a simple measure of shape complexity.

Metrics Calculation and Range Description
Mean Contiguity Index (CONTIG_MN)

N
c ijr = contiguity value for pixel r in patch ij; v = sum of the values in a 3-by-3 cell template; a ij = area of patch ij in terms of number of cells.
CONTIG equals the average contiguity value for the cells in a patch (i.e., sum of the cell values divided by the total number of pixels in the patch) minus 1 and divided by the sum of the template values (13 in this case) minus 1.The contiguity index assesses the spatial connectedness, or contiguity, of cells within a grid-cell patch to provide an index on the patch boundary configuration and thus the patch shape.

Mean Euclidian
Nearest Neighbor Distance (ENN_MN) from patch ij to nearest neighboring patch of the same type (class).ENN > 0; Units, meters.
ENN_MN is the mean value of the distance (m) to the nearest neighboring patch of the landscape, based on the shortest edge-to-edge distance.Note that the edge-to-edge distances are from cell center to cell center.ENN_MN measures the isolation of the patches of the landscape.The straight line distance is the shortest path between the patches.

Interspersion and Juxtaposition Index (IJI)
e ik = total length of edge in landscape between patch types i and k; m = number of patch types present in the landscape.1 < IJI ≤ 100; Units, percent.
IJI equals a negative sum of the length (m) of each unique edge type divided by the total landscape edge (m), multiplied by the logarithm of the same quantity, and summed over each unique edge type; divided by the logarithm of the number of patch types, times the number of patch types, minus 1, and divided by 2; and multiplied by 100 (to convert to a percentage).IJI provides a measure of isolating the interspersion or intermixing of patch types.

Patch Cohesion Index (COHESION) COHESION
p ij = perimeter of patch ij in terms of number of cell surfaces; a ij = area of patch ij in terms of number of cells; Z = total number of cells in the landscape.1 < COHESION < 100 COHESION equals 1 minus the sum of the patch perimeter (in terms of number of cells), divided by the sum of the patch perimeter, times the square root of the patch area (in terms of number of cells) for all patches in the landscape, divided by 1 minus 1 over the square root of the total number of cells in the landscape, and multiplied by 100 to convert to a percentage.The patch cohesion increases as the patch type becomes more clumped or aggregated in its distribution.

Statistical Analysis
For finding differences in different RD regions, the analysis of variance (ANOVA) together with multiple comparisons were performed using the SPSS software (version 19.0, 2010, IBM, USA).Before ANOVA, the normal distribution and homogeneity of the variance of the experimental data (all variables) were statistically analyzed.We use QQ-plots to test whether the raw data and the transformed data conform to the normal distribution and to verify by frequency distribution histograms (Figures S1-S3).For data that did not conform to the normal distribution or the approximate normal distribution, we performed logarithmic transformation, reciprocal transformation, and so on, until the data were in line with the normal distribution (p > 0.05) or approximate normal distribution (p > 0.01) (Figure S3).Then, ANOVA and multiple comparison tests were used to determine differences in road characteristics, forest landscape patterns, and forest structural-taxonomic attributes across three RD regions.When variances were homogeneous, we used Duncan's tests; Games-Howell tests were used otherwise.
The associations between road development and forest landscape pattern were generally explored through correlation analysis and curve fitting [26,27], as well as principal component analysis [17].These analyses together with redundancy analysis (RDA) were used in this study to explore the associations between road intensity and urban forest landscapes as well as structural-taxonomic attributes.
RDA ordination and variation partitioning have been used to decouple complex associations among various ecological factors [23,28,29].Moreover, these methods were used in this study to Forests 2019, 10, 58 9 of 23 identify the relative contributions of road characteristics and urban forest landscape metrics to forest structural and taxonomic variations, and a two-group analysis ("Var-part-2groups-simpleeffects-tested" in Canoco (v.5.0; Biometrics, The Netherlands) was conducted.Using structural attributes and taxonomic attributes as dependent factors, two groups of explanatory factors (RD attributes and landscape patterns) were used to find their relative contribution to variations in these dependent factors.In the partial permutation test, the candidate variable was used as the only explanatory variable, and variables already selected were used as covariates.The effect of the variable tested in such a context is referred to as its conditional (or partial) term effect.At the start of the forward selection process, when no explanatory variables had been entered into the selected subset, each variable was tested separately to estimate its independent, simple (marginal) term effect.During RDA ordination processes, both simple term effects and conditional term effects were included in the study to assess the parameters that were most likely responsible for urban forest (structural and taxonomic) variations in different RD regions (low, medium, and heavy).Canoco (v.5.0; Biometrics, The Netherlands) was used for RDA analyses.
Regression analyses were performed using landscape metrics, structural attributes, taxonomic attributes (dependent variables), and RD (independent variables), and linear or polynomial relationships were used in the regression.Raw data for each grid is listed in Table S1, and their quality was checked by the standard error of the mean, standard deviation, and total replicating measurements in each grid.Stepwise regression was also performed to identify associations between structural attributes, taxonomic attributes, and landscape metrics in different RDs.The entering criterion for the selected variables in the regression equations was defined as <0.05 for the F value.The coefficients of determination (R 2 ), standard coefficient for each entered parameter, and p-values (p) were determined for each model to compare the relative contribution and best-fitness of the model.After the equation was established by regression analysis, the model was diagnosed (heteroscedasticity, independence of the residuals, multicollinearity test, and so on, which has shown in Table S2 and Figure S5).All statistical analyses were carried out using SPSS version 19.0 (SPSS Inc., Chicago, IL, USA).

Characterization of the Study Object in Road Characteristics
Characterization of road development at 3 RD regions are shown in Table 2.The mean values of four road characteristics (RAD, RAI, RLI, and CPD) exhibited an increasing trend from low road RD to heavy RD areas, and these increasing trends were linearly and statistically significant (p < 0.05).Furthermore, the mean values of the five road characteristics significantly differed at different RDs (p < 0.05) (Table 2).In contrast with the other five parameters, the ISA percentage in the three road regions were not statistically significant, and no linear changes were observed (Table 2).In summary, the roads in Harbin were characterized as the linear increase in road area and road length for both the total terrestrial land area and the total impervious surface area (Table 2).

Landscape Metrics of Urban Forests
Landscape heterogeneities across different RDs are depicted in Figures 3 and 4, and regression analysis results between road density and landscape metrics are showed in Figure S4.In the case of the area-edge and shape metrics of urban forests, only TA showed significant changes at low, medium, and heavy RD regions, with the exception of PARA-MN (Figure 3).Mean values of TA in low RD regions were approximately 34.85 ha lower than TA in heavy RD regions (98% decline, p < 0.05), and regression analysis indicated a logarithmic increase in RDs (p < 0.05).This indicates that the high RD regions in Harbin have a larger total forest landscape area than the low RD regions.Besides this significant change, an increasing tendency in ED with RDs was also found in the mean values, while a decreasing tendency was observed in CONTIG-MN (Figure 3).This indicates that the fragmentation degree of forest landscape edge in the high RD regions in Harbin was higher than that in the low RD regions, but the complexity degree of the forest patch in the high RD regions was lower than that in the low RD regions.

Characterization of the Study Object in Road Characteristics
Characterization of road development at 3 RD regions are shown in Table 2.The mean values of four road characteristics (RAD, RAI, RLI, and CPD) exhibited an increasing trend from low road RD to heavy RD areas, and these increasing trends were linearly and statistically significant (p < 0.05).Furthermore, the mean values of the five road characteristics significantly differed at different RDs (p < 0.05) (Table 2).In contrast with the other five parameters, the ISA percentage in the three road regions were not statistically significant, and no linear changes were observed (Table 2).In summary, the roads in Harbin were characterized as the linear increase in road area and road length for both the total terrestrial land area and the total impervious surface area (Table 2).p < 0.05), and regression analysis indicated a logarithmic increase in RDs (p < 0.05).This indicates that the high RD regions in Harbin have a larger total forest landscape area than the low RD regions.Besides this significant change, an increasing tendency in ED with RDs was also found in the mean values, while a decreasing tendency was observed in CONTIG-MN (Figure 3).This indicates that the fragmentation degree of forest landscape edge in the high RD regions in Harbin was higher than that in the low RD regions, but the complexity degree of the forest patch in the high RD regions was lower than that in the low RD regions.Compared to the area-edge and shape metrics, the landscape aggregation metrics of urban forests showed more significant changes among different RDs (Figure 4).LSI and PD significantly differed across low to heavy densities (Figure 4B,D).Regression analysis results showed that LSI and PD logarithmically increased along RDs (p < 0.01) (Figure 4B,C).Regression equations indicated that a rapid increase in LSI and PD appeared in the lower RD region (<1.0 km/km 2 ) while the increase became much slower at higher RD regions (Figure 4, Table 3).The mean values of ENN-MN showed a steadily decreasing trend along the RDs.However, this variation did not differ significantly across different RDs (p > 0.05) (Figure 4).
Overall, road intensification in Harbin urban regions aligned with that of increasing urban forest areas (TA), and more patch numbers and more diversified patch shapes (LSI) with closer patch distances (ENN-MN) were observed (Figures 3 and 4; Table 3).

Structural and Taxonomic Attributes of Urban Trees
Forest structural attributes at different RDs are depicted in Figure 5, and the regression analysis of the road density with all the forest taxonomic attributes are showed in Figure S4.Among these structural attributes, only TD significantly differed across different RDs (Figure 5).Tree height logarithmically decreased with RDs, and the decreased rate was equivalent to 1.10 m per 1 km/km 2 increase in RD (Figure 5C and Table 3).DBH section areas linearly decreased from low RD to heavy RD, and logarithmical regression equations in Table 3 show that a unit increase in RD was accompanied by a 8.14 cm 2 /m 2 logarithmical decrease in the section area (p < 0.01) (Figure 5D).In addition, DBH and TD decreased with increased RD (Figure 5E,F); however, regression analyses of the raw data between RD and these two parameters showed non-significant changing trends (p > 0.05).Overall, the road intensification in Harbin was usually accompanied by smaller trees (lower tree height, DBH, and TD), resulting in smaller wood section areas (DBH section area).
of the road density with all the forest taxonomic attributes are showed in Figure S4.Among these structural attributes, only TD significantly differed across different RDs (Figure 5).Tree height logarithmically decreased with RDs, and the decreased rate was equivalent to 1.10 m per 1 km/km 2 increase in RD (Figure 5C and Table 3).DBH section areas linearly decreased from low RD to heavy RD, and logarithmical regression equations in Table 3 show that a unit increase in RD was accompanied by a 8.14 cm 2 /m 2 logarithmical decrease in the section area (p < 0.01) (Figure 5D).In addition, DBH and TD decreased with increased RD (Figure 5E,F); however, regression analyses of the raw data between RD and these two parameters showed non-significant changing trends (p > 0.05).Overall, the road intensification in Harbin was usually accompanied by smaller trees (lower tree height, DBH, and TD), resulting in smaller wood section areas (DBH section area).Abbreviations: patch density (PD), landscape shape index (LSI), Shannon Wiener diversity index (H'), species richness (SR), species evenness index (J'), section area at breast height (DBH section area), and total area in the grid (TA).
Forest taxonomic attributes across different RDs are displayed in Figure 6, and the regression analysis of the road density with all the forest structural attributes are showed in Figure S4.H' and SR values differed significantly across different RDs (p < 0.05).The regression analyses indicated that SR, H', and J' logarithmically increased with RD.The increased logarithmic rates were 2.07, 0.32, and 0.1, respectively, at one unit increases in the road length density fraction (Figure 6 and Table 3).Overall, the road intensification in Harbin was accompanied by smaller tree sizes (particularly in height), lower TD (both in section area and bole number), and more diversified tree species compositions (SR, H', and J') (Figures 5 and 6; Table 3).

Association Decoupling
To determine the relative contributions of the road characteristics and forest landscape metrics to the variation of the forest structural-taxonomic characteristics, RDA ordination-based variation partitioning analyses were performed (Figure 7).Regarding structural attribute variations, road characteristics alone could explain 51.2% of the structural attribute variations, and this was 2.7-fold higher than those of the unique characteristics that could explain the influence of landscape metrics alone.The interactions between the RD characteristics and landscape metrics could explain 30.1% of the variation (Figure 7).In the case of the variations in forest taxonomic attributes, the explanatory power of the RD characteristics (38.1%) was 1.6-fold higher than that of the unique landscape metrics (24.2%).Moreover, the interactions between the RD characteristics and forest landscape metrics could explain 37.6% of the variation (Figure 7).
Overall, RD characteristics could powerfully explain the variation in the forest structural attributes and taxonomic attributes compared to the forest landscape metrics.Moreover, the influences of urban forest landscapes on forest characteristics were more likely achieved via their interactions with RD characteristics.To identify the importance of landscape metrics with regards to the regulation of forest characteristics, different RD regions should be classified to emphasize the importance of landscape metrics (Figure 7).
Overall, RD characteristics could powerfully explain the variation in the forest structural attributes and taxonomic attributes compared to the forest landscape metrics.Moreover, the influences of urban forest landscapes on forest characteristics were more likely achieved via their interactions with RD characteristics.To identify the importance of landscape metrics with regards to the regulation of forest characteristics, different RD regions should be classified to emphasize the importance of landscape metrics (Figure 7).

Significant Parameters that Explain Forest Variation
Based on all the forest characteristics, the simple term effects and conditional effects from the RAD ordination were used to identify the significant factors responsible for the variations in the forest characteristics in different RD regions (Figure 8).In low RD regions, the taxonomic attributes of forests were negatively related to the forest structure parameters, i.e., a higher diversity of the forest was associated with smaller tree sizes (TH) and lower wood density (section area).Moreover, several parameters of landscape metrics (ED, IJI, etc.) and road characteristics (RAD, RLD, and CPD) were closely and positively related to the species diversity indices (H', J', and SR) (Figure 8a).The RDA axis extracted 50% of the information from raw data.In the case of the simple term effects, three road characteristics (RAD, RLD, and CPD) were identified as three significant explaining factors for variations in the forest characteristics (p < 0.05), followed by significant factors related to the landscape patterns of urban forests (IJI and ED); however, these factors were not statistically significant, (p = 0.076 and 0.06, respectively) (Figure 8b).When excluding the collinear effects (conditional effects), one road characteristic (RAD) and one forest landscape metrics (PARA_MN) exhibited significant associations with the forest structural and taxonomic attributes (p < 0.05) (Figure 8).
In the medium and heavy RD regions, much weaker road-related associations with the forest structural-taxonomic attributes were observed than those in low RD regions.The RDA axis extracted 27.4% of the information from raw data, which is approximately half of that identified in the low RD region.The length of the arrows for road attributes (RDA, RLD, and CPD) was similar to those from forest landscape metrics (IJI, PD, and LSI).However, in the low RD region, the lengths were much longer than those from the latter, indicating much stronger influences from the former (Figure 8).

Significant Parameters that Explain Forest Variation
Based on all the forest characteristics, the simple term effects and conditional effects from the RAD ordination were used to identify the significant factors responsible for the variations in the forest characteristics in different RD regions (Figure 8).In low RD regions, the taxonomic attributes of forests were negatively related to the forest structure parameters, i.e., a higher diversity of the forest was associated with smaller tree sizes (TH) and lower wood density (section area).Moreover, several parameters of landscape metrics (ED, IJI, etc.) and road characteristics (RAD, RLD, and CPD) were closely and positively related to the species diversity indices (H', J', and SR) (Figure 8a).The RDA axis extracted 50% of the information from raw data.In the case of the simple term effects, three road characteristics (RAD, RLD, and CPD) were identified as three significant explaining factors for variations in the forest characteristics (p < 0.05), followed by significant factors related to the landscape patterns of urban forests (IJI and ED); however, these factors were not statistically significant, (p = 0.076 and 0.06, respectively) (Figure 8b).When excluding the collinear effects (conditional effects), one road characteristic (RAD) and one forest landscape metrics (PARA_MN) exhibited significant associations with the forest structural and taxonomic attributes (p < 0.05) (Figure 8).
In the medium and heavy RD regions, much weaker road-related associations with the forest structural-taxonomic attributes were observed than those in low RD regions.The RDA axis extracted 27.4% of the information from raw data, which is approximately half of that identified in the low RD region.The length of the arrows for road attributes (RDA, RLD, and CPD) was similar to those from forest landscape metrics (IJI, PD, and LSI).However, in the low RD region, the lengths were much longer than those from the latter, indicating much stronger influences from the former (Figure 8).Both the simple term effects and the conditional term effects indicated that no parameters had statistically significant explanatory power with regard to the forest structural-taxonomic variations in medium and heavy RD regions (Figure 8b).
Overall, road characteristics, particularly RAD, exhibited the most significant associations with the forest characteristics in low RD regions.However, such associations became weaker with more road development in medium and heavy RD regions, and the forest landscape-related associations become stronger based on the close-to-significant explanatory power of PARA-MN (p = 0.062) (Figure 8b).Stepwise regression analysis showed that landscape metrics and road characteristics had different associations with regards to the variable structural-taxonomic attributes, and low RD regions usually had much stronger associations than medium-heavy RD regions (Table 4), indicating a cross-method validation with the RDA ordination (Figure 8).In low RD regions, seven out of the eight forest structural-taxonomic attributes (Height, DBH section area, TD, UBH, canopy size, H', and Stepwise regression analysis showed that landscape metrics and road characteristics had different associations with regards to the variable structural-taxonomic attributes, and low RD regions usually had much stronger associations than medium-heavy RD regions (Table 4), indicating a cross-method validation with the RDA ordination (Figure 8).In low RD regions, seven out of the eight forest structural-taxonomic attributes (Height, DBH section area, TD, UBH, canopy size, H', and SR) were significantly associated with landscape metrics and road characteristics, and r 2 ranged from 0.28 to 0.87.However, in medium to heavy RD regions, five out of the eight parameters showed significant associations, and R 2 was generally less than 0.29 (Table 4).
As shown in the multicollinearity diagnosis in Table S2, we used the Eigenvalue and Condition Index as the indexes to judge multicollinearity.The multidimensional Eigenvalue at about 0 and the condition index > 10 indicates the existence of multicollinearity (SPSS user manual), and in most cases of our studies, the endpoints of the stepwise regression is the existence of multicollinearity.For example, after 6 dimensions in the DBH sections area-related models, both the Eigenvalue and Condition Index showed an existence of collinearity among the dependent variables, and then the stepwise regression model ended in the 7 dimensions of the regression.All other variables showed much lower collinearity than the DBH at the low RD region (Table S2), indicating that our stepwise model makes sense of excluding multicollinearity.
We also used the Durbin-Watson value as a measure of the residual independent diagnosis of the stepwise regression models (SPSS user manual).In general, the closer the Durbin-Watson value is to 2, the less biased residual terms are.The closer the Durbin-Watson value is to 0, the stronger positive biased residual terms will be.The closer the Durbin-Watson value is to 4, the stronger the negative biased residual terms will be.As shown in Table S2, the Durbin-Watson values for all the stepwise models are all around 2, which basically determines that the residuals are independent, and an unbiased estimation was achieved in our stepwise models.This kind of unbiased estimation can also be found in Figure S5.

Table 4.
Stepwise regressions between forest attributes and road characteristics (structural attributes and taxonomic attributes) and landscape metrics in the different road density (RD) regions.The selected variables were selected at a probability of less than 0.05 for the F value.The coefficients of determination (R 2 ) and p-value (p) are also given for each model, and the data are listed in this table.Contrary to previous studies, our findings highlighted that road intensification increased the forest area, mainly through increased patch numbers and patch shape complexity.Previous studies found that road networks have both positive and negative ecological effects, but the majority of these effects are negative [2].Increases in RD and urbanization will lead to decreased forest patch areas, increased quantities, landscape fragmentation, and intensification of landscape heterogeneity [32,48,49].RD is a measure of urbanization and connection of various landscapes, which is affected by both forest landscapes and forest attributes.Moreover, RD impacts most forest structural-taxonomic attributes compared to landscape metrics, and this result indicated a much stronger explanatory power on the variations of forests from road attributes.This is consistent with the conclusion of previous studies and our hypothesis.The fragmented urban forest could affect landscape beautification, recreational services [50], tree and bird species conservation [32,51], tree and soil carbon stock functions [15], etc.Our results indicated that long-term plantation practices could increase forest areas during road development in urban regions (although fragmentation might occur).

Region
Trees play a large role in urban environment regulations because of their large body size-related vertical regulation in multiple environmental functions [14,28].In this study, among forest structure attributes, the tree height and DBH section areas in low RD regions were much larger than those in heavy RD regions, and the section areas were significantly negatively correlated with RD.This result is consistent with our hypotheses and the findings of Fahey and Casali [48] who examined forest ecosystem distributions for over two centuries in heavy urbanization regions.The results of their study also indicated that remnant forests had higher canopy cover and basal areas than recently established forests [48].In areas with convenient transportation and economic development, more exotic tree species were often afforested for environmental beatification, so the tree height and basal area were smaller than the forests in areas with fewer transportation facilities [52], where legacy trees could grow for a long time.Based on our study in Harbin, heavy RD areas should enhance the protection of old trees (larger tree size), and low RD areas should appropriately increase the planting of trees to achieve a reasonable urban forest density (basal DBH area per ha).
Regarding forest taxonomic attributes, the H' and SR values were significantly higher in heavy RD regions than in low RD regions.However, this is contrary to our hypothesis that road development leads to a reduction in the urban forest taxonomic attributes.Previous studies in western countries showed that people prefer to live in medium developed areas [53][54][55]; therefore, tree species diversity was highest in these areas.Currently, people in China tend to live in the downtown regions to obtain convenient daily lives [24,56].This preference possibly results in the planting of more diverse species in highly developed areas with more road networks, corresponding to high species diversity, species richness, and evenness as observed in this study.The species diversity and richness of trees can be enhanced by strengthening intensive protection and suitable management practices of various urban forests [24,57,58].
The association between road development and forest characteristics is possibly due to China's governance policy.For example, afforestation is highly encouraged because of the shortage of forests in China, and the harvesting or cutting of trees is strictly inhibited without official permission, which is regulated by forest laws (https://baike.baidu.com/item/中华人民共和国森林法/719190).The long-term afforestation policy in China has been implemented for over 40 years, and forest coverage has increased from 8% to 22% presently.Along-road afforestation is an important part of afforestation practices [11] and is a must-do action in the process of road construction in China.In China, a detailed road-side afforestation plan should be proposed before road construction, be monitored during road development, and be fully checked by administrative authorities after road construction.The national-guided standards (Code for Planting Planning and Design on Urban Road (CJJ75-97) (https://baike.baidu.com/item/《城市道路绿化规划与设计规范》GJJ75-97/14703300))have been implemented for over 20 years.These regulations and laws in China have potentially increased the close association between forest attributes, including compositional traits and road development, as observed in this study.Thus, forest-related studies, evaluation, and management in China should fully consider road development.Urbanization has become a more worldwide homogenization [59], and road development are the important infrastructure of urbanization.Other cities in the world had similar historical backgrounds, geographical locations, climatic conditions, and economic development statuses.In particular, similar afforestation practices implemented for the long term may also show similar road-development-related changes and forest (both landscape and tree community) characteristics.Their applicability in other parts of the world should be investigated in future studies.

Road-Dependent Landscape Regulation Could Improve Forest Ecological Services
The existence and expansion of road networks could affect the surrounding landscape ecology, patterns, and processes, potentially impacting regional ecology security [60].Landscape metrics could be used as indicators of forest management effects [19,20].Moreover, the indicator functions of landscape metrics in association with ecosystem services are increasingly favored by ecologists [61,62].Landscape heterogeneity metrics, which are crucial to the configuration of landscape patterns and ecosystem regulation, are particularly favored [14].In this study, we found that unique effects of road attributes as indicators of RD characteristics in Harbin could explain over half (51.2%) of the variation observed in forest structural attributes as well as 38.1% of the variation in forest taxonomic attributes (Figure 7).Furthermore, forest landscape metrics and their interactions with road attributes were responsible for the remaining 48.8%-61.9% of structural attributes and forest taxonomic attributes.The interaction between urban forest landscapes and RD, rather than landscape alone, provided better explanatory power with regards to variations in forest structural and taxonomic attributes.In our study, TA, PD, and LSI had maximum values in the heavy and medium RD areas, indicating a tendency for landscape fragmentation in urban forests in the most urbanized regions of Harbin.Any forest-landscape regulation measures should take RD into consideration.Thus, different methods in different RD regions are potentially the most reliable way to implement landscape-related measures that can be used to optimize forest structure and composition.
The structural attributes of trees are the basis of the multiple ecological services provided by urban forests, such as increasing carbon sinks [63,64], reducing air temperature [28,42,65], and decreasing storm water runoff [66].In the low RD regions of Harbin, taller and bigger trees were more likely to be associated with lower patch aggregation (IJI) and lower patch areas (TA), as shown by the stepwise regression (Table 4) and RDA ordination results (Figure 8).However, road intensification (in medium and heavy RD regions) largely weakened these associations (Table 4).In the heavy RD regions, larger trees were not observed, and this was likely the result of road construction, which usually occurs at the expense of remnant old forests.Alternatively, young trees were usually afforested at roadsides after road construction, in accordance with local regulations and national laws.In the low RD regions, the decreasing aggregation of forest landscapes possibly favored the conservation of larger urban trees, and this coincided with the general observations of tall and large sparse trees in the rural regions and old parks in Harbin (https://heilongjiang.dbw.cn/system/2017/06/09/057670066.shtml).Therefore, regarding the management of urban forests, we should pay more attention to RD, and different RD regions and landscape metrics should be used to regulate and predict urban tree size-related vertical structures.
Biodiversity had a moderating effect, and it affected all levels of the ecosystem [29,67].Previous studies showed that the shapes and sizes of the patches were strongly associated with the ecosystem service functions, and the degree of landscape fragmentation was negatively correlated with the species diversity indices [24,68].Landscape metrics as a significant indicator function of woody plants and overall species richness [69] were once used in tropical forest research, which concluded that patch shape and similarity metrics can be used as indicators to predict species diversity [67].In this study, lower species diversity, richness, and evenness values were associated with lower road densities, patch shape complexities, patch aggregation, and edge densities (p < 0.05); however, these associations were much weaker in medium to heavy RD regions (Figure 8).Similarly, stepwise regression results indicated that the RD attributes exhibited a much stronger association with the diversity variation in low RD regions, while the forest landscape metrics of the interspersion features (IJI) and shape features (LSI) were included in the regression models of other regions (Table 4).Therefore, in low RD regions, the forest species compositional traits were most likely regulated by RD.In contrast, in heavy RD regions, the road-related regulations became weaker, and the forest patch attributes were associated with the regulation of the tree composition.These results provided hints for the suitable road-dependent landscape regulation of urban forest species compositions.

Technical Implications and Uncertainty
Our findings in this study indicate that road intensification parameters (i.e., road intensity related to road length and road area) in urban regions were closely associated with urban forest characteristics.In previous studies, urbanization intensity, classified using impervious surface percentage, was used to study associations between forest structural-taxonomic attributes and forest landscape metrics [24].In contrast to impervious surface data (extracted using remote sensing), road data are much more approachable via internet providers, including Google (http://www.google.cn/maps),Baidu (https://map.baidu.com/),and Tencent (https://map.qq.com/) and via professional road map providers such as Gaode map (https://www.amap.com/).With these free and available road data and our findings in this study (i.e., identification of the most indicatable road parameters that predict urban forest structural and taxonomic traits, RLD, and RAD), it is possible to make a more precise evaluation of urban forest changes during urbanization processes.Our results strongly suggested that in the future, this kind of bigdata on the internet could be technically used in ecological studies of road ecology.Moreover, some stepwise regression models (Table 4) could possibly be used in these evaluations, and the most probable indicator for each forest characteristic could also be found in the ordination map (Figure 8).
Decoupling complex associations have challenges in field ecological studies.For instance, in urban regions, both natural conditions (soil nutrients, climatic conditions, and species competition) and social development (e.g., road construction and economic development) should be fully considered to obtain a scientific understanding of the underlying process and mechanisms.The ordination and stepwise regression analyses have been used in various studies to decouple associations [42,70].Although statistical analyses were performed in this study, conclusions regarding possible causal relationships require additional caution (Table 4, Figures 7 and 8).In the future, more advanced statistical methods (e.g., structural equation models) will require the identification of direct and indirect effects [23,30,71,72].Moreover, long-term detailed field-controlled experiments (such as large plots in natural forest studies [73]) will be useful when testing the hypotheses proposed in this study [74].
External factors and experimental design should also be considered for the uncertainty of this study.Besides road development tested in this paper, other factors possibly affecting forest characteristics include environmental factors (pollution, climate change, soil compaction, and impervious surface in urban regions) [24] and human disturbances (planting, cutting, and cultivation) [75].In the data analysis and discussion, most of the environmental factors, such as climate and pollution, were assumed as random factors.Furthermore, the experimental design in this paper utilized a sample size of 400 m 2 for the tree census.The larger plot will increase tree species richness and change species diversity indices as well as tree size parameters.In this paper, this kind of scale dependence was not investigated.In the future, more detailed experiments are warranted to clarity this uncertainty.

Figure 1 .
Figure 1.Location of the studied Harbin city (a) and main processes of road and forest landscape data extraction (b).Abbreviations : Gaofen No.1 remote sensing image (GF-1 image), Redundancy analysis (RDA), Analysis of variance (ANOVA).

Figure 1 .
Figure 1.Location of the studied Harbin city (a) and main processes of road and forest landscape data extraction (b).Abbreviations: Gaofen No.1 remote sensing image (GF-1 image), Redundancy analysis (RDA), Analysis of variance (ANOVA).

Figure 2 .
Figure 2. Urban road distribution in Harbin (a); classification of low, medium, and heavy road density (RD) regions (b); and road extraction precision with reference to Google Roadnet map (c).

Figure 2 .
Figure 2. Urban road distribution in Harbin (a); classification of low, medium, and heavy road density (RD) regions (b); and road extraction precision with reference to Google Roadnet map (c).

Figure 3 .
Figure 3. Area-edge and shape metric changes in urban forests across different road densities (RD).(A) area-weighted mean contiguity index (CONTIG-MN), (B) mean fractal dimension index (FRAC-MN), (C) mean perimeter area ratio (PARA-MN), (D) edge density (ED), and (E) total area (TA).Boxplots show the mean (black square), median (black line in the box), inter-quartile range (25%-75% in the box), and the range of normal values (the whisker).Mean values sharing different letters are significantly different (p < 0.05), and mean values sharing the same letters are not significantly different (p > 0.05).Figure insets depict significant regression lines from the raw data.The red line represents the regression line, and the green lines represent the moving mean value.The dual-shaded shadows represent the best-fitting and predicted interval ranges.

Figure 4 .Figure 4 .
Figure 4. Urban forest landscape aggregation changes from low to heavy road density regions.(A) patch cohesion index (COHESION), (B) landscape shape index (LSI), (C) patch density (PD), (D) interspersion and juxtaposition index (IJI), and (E) mean Euclidian nearest neighbor distance (ENN-MN).The boxplots show the arithmetic mean (black square), median (black line in the box), interquartile range (25%-75% at the box upper and lower edge), and the normal value range (the whisker).The mean values sharing different letters are significantly different (p < 0.05), and different letters at different RD regions are significantly different in their arithmetic means (p < 0.05).The figure insets depict significant regression lines from the raw data.The red line is the regression line, and the green

Forests 2019 ,Figure 6 .
Figure 6.Taxonomic attributes at different road densities (RD).(A) Shannon Wiener diversity index (H'), (B) species richness (SR), and (C) evenness index (J').The boxplots show the mean (black square), median (black line in the box), inter-quartile range (25%-75% in the box), and the range of normal values (the whisker).The mean values sharing different letters are significantly different (p < 0.05), and the different letters at different RD regions are significantly different in their arithmetic means (p < 0.05).The figure insets depict significant regression lines from the raw data.The red line is the regression line, and the green lines represent the moving mean value.The dual-shaded shadows represent the best-fitting and predicted interval ranges.

Figure 6 .
Figure 6.Taxonomic attributes at different road densities (RD).(A) Shannon Wiener diversity index (H'), (B) species richness (SR), and (C) evenness index (J').The boxplots show the mean (black square), median (black line in the box), inter-quartile range (25%-75% in the box), and the range of normal values (the whisker).The mean values sharing different letters are significantly different (p < 0.05), and the different letters at different RD regions are significantly different in their arithmetic means (p < 0.05).The figure insets depict significant regression lines from the raw data.The red line is the regression line, and the green lines represent the moving mean value.The dual-shaded shadows represent the best-fitting and predicted interval ranges.

Forests 2019 , 25 Figure 7 .
Figure 7. Partitioning of the RDA ordination-based forest attribute variation for the entire urban region.(A) Variation partitioning for structural attributes and (B) variation partitioning for taxonomic attributes.

Figure 7 .
Figure 7. Partitioning of the RDA ordination-based forest attribute variation for the entire urban region.(A) Variation partitioning for structural attributes and (B) variation partitioning for taxonomic attributes.

Forests 2019, 10 ,Figure 8 .
Figure 8. RDA ordination on forest structural-taxonomic attributes from landscape metrics and road characteristics in different road density (RD) regions.(a) low RD region and (b) medium and heavy RD regions.

Table 2 .
Characterization of the study object at different density levels.

Table 2 .
Characterization of the study object at different density levels.

Table 3 .
The relationships between landscape metrics, structural-taxonomic attributes, and road densities (RDs) measured as road density fractions.

Table 3 .
The relationships between landscape metrics, structural-taxonomic attributes, and road densities (RDs) measured as road density fractions.

Items Equations with RDs as X (km/1000 m 2 ) Logarithmic Slope or Exponent R 2 p-Value
Figure 8. RDA ordination on forest structural-taxonomic attributes from landscape metrics and road characteristics in different road density (RD) regions.(a) low RD region and (b) medium and heavy RD regions.

Dependent Variables Final Model R 2 p-Value
Road Development Associated with Increased Forest Areas Characterized by More Patches, Complex Shapes, Smaller Tree Sizes, Lower Density, and More Diversified Species Compositions