Comparing Negative Impacts of Prunus serotina , Quercus rubra and Robinia pseudoacacia on Native Forest Ecosystems

: The introduction of invasive alien plant species (IAPS) can modify plant-soil feedback, resulting in an alteration of the abiotic and biotic characteristics of ecosystems. Prunus serotina , Quercus rubra and Robinia pseudoacacia are IAPS of European temperate forests, where they can become dominant and suppress the native biodiversity. Assuming that the establishment of these invasive species may alter native forest ecosystems, this study comparatively assessed their impact on ecosystems. This study further investigated plant communities in 12 forest stands, dominated by the three IAPS and native trees, Quercus robur and Carpinus betulus (three plots per forest type), in Northern Italy, and collected soil samples. The relationships between the invasion of the three IAPS and modiﬁcations of humus forms, soil chemical properties, soil biological quality, bacterial activity and plant community structure and diversity ( α -, β -, and γ -diversity) were assessed using one-way ANOVA and redundancy analyses (RDA). Our comparative study demonstrated that invaded forests often had unique plant and / or soil properties, relative to native forests, and the degree of dissimilarity depended on the invasive species. Particularly, Q. rubra is related to major negative impacts on soil organic horizons and low / modiﬁed levels of microarthropod and plant biodiversity. R. pseudoacacia is associated with an altered base content of soil and, in turn, with positive feedback to the soil biological quality (QBS-ar) and plant diversity, but with a high cover compared with other alien plant species. P. serotina is associated with intermediate impacts and exhibits a plant species assemblage that is more similar to those of native forest stands. Our work suggests impact-based management decisions for the three investigated IAPS, since their e ﬀ ects on the diversity and composition of resident ecosystems are very di ﬀ erent. were conducted to test the homogeneity of the and normality assumption, The correlation plots were performed the “corrplot” R


Introduction
The establishment of invasive alien plant species (IAPS) can significantly change the structure and function of ecosystems [1]. It is well known that many IAPS affect resident plant communities in terms of dynamics, dominance, diversity, and structure [2,3]. IAPS can also modify plant-soil feedback, which alters the abiotic and biotic characteristics of soil, such as pH, moisture, nutrient availability and faunal and microbial communities [4,5]. Particularly, soils colonized by alien species often show increased proportions of organic carbon, total nitrogen, phosphorus and biomass through litterfall, in comparison with those developed under the dominance of native species [6]. These factors may (a). To evaluate the possible differences in humus forms and subsequent variations in other soil physical and chemical properties among forest stands; (b). To determine if variations in biotic components occurred among forest stands, particularly ascertaining possible reductions in biodiversity levels and biological soil quality; and (c). To recognize the relationships among biotic and abiotic determinants across the different forest stands.

Study Area and Sites
The study was carried out in temperate deciduous woodland patches of the Po plain, in the Lombardy region (Italy; 45 •  The original land cover of the area is Padanian-Illyrian hardwood forest in transition to mesophytic Padanian mixed oak forest [27], dominated by deciduous oaks, Quercus spp., and hornbeam Carpinus betulus. Nowadays, the area is highly urbanized, with a residual ancient forest cover that is, in part, threatened by the invasion of neophyte alien species, such as P. serotina, Q. rubra and R. pseudoacacia. Q. rubra, native to North America, was introduced in Europe between the XVII and XVIII centuries, as ornamental and timber species, while the first evidence of its cultivation in Italy dates back to 1812. P. serotina, native to central-eastern North America, was introduced in Italy and Lombardy in 1922 [28]. Both species have been studied for forestry purposes since 1922, with different results. Red oak has been used in tree plantations, artificially introduced into the native forest, after tree cutting and by planting or sowing. From there, it quickly spread as consequence of its ability to compete with shade-tolerant species and its rapid growth (faster than the common oak) and resistance to water stress and disease. The experimental fields of black cherry were instead rapidly abandoned due to the poor quality of the produced wood. From these plantations, the surviving individuals started an intensive invasion, spreading in the same way that the red oak did, but also randomly spreading with the help of birds.
R. pseudacacia, native to eastern North America, was introduced into Europe in 1601 and into Italy in 1662. The forestry interest in this species favored its spreading due to its suitability for the woodland coppice. As it provided good quality timber, with a high growth rate, the native species were progressively cut and displaced.
Over an area of approximately 250 km 2 , four forest types were investigated: (a) Q. robur and Carpinus betulus native forest (QF, Quercus forest), (b) P. serotina forest (PS), (c) Q. rubra forest (QR) and (d) R. pseudoacacia forest (RP). For each forest, 3 randomly selected plots (=sites) of 20 × 20 m were surveyed, with a total of 12 independent sites. Each forest type was defined when the dominant species had a canopy cover higher than 80%. The distance between the plots ranged from approximately 0.15 to 20.5 km. With regard to the forest stands invaded by alien species, all the forests are between 60 and 65 years old, except for one stand (i.e., plot) of the RP forest, which is 40 years old, as estimated by dendrocronologiacal analyses of collected increment cores using a Pressler borer. All the study sites belonged to the same pedological landscape unit [29], allowing the assumption that they had the same abiotic and biotic characteristics before the invasion. They are located within the fluvioglacial würmian plain of the Ticino river basin, with acidic pebbly-gravelly substrates [30].
According to the World Reference Base (WRB) classification [31], the soils mainly belong to the Cambisols reference group. They are desaturated in bases, often rich in fragments and have a subacid to acid reaction. The elevation of the study sites ranged from 192 to 209 m a.s.l.
The average annual precipitation and temperature at the Busto Arsizio meteorological station (close to the study area) were 1081 mm and 13.8 • C, respectively (ARPA Lombardia, mean for the period 2006-2016).

Soil Abiotic Characteristics
In each plot, a soil minipit was opened down to C horizon and described for the taxonomic classification [31]. The soil sampling was performed during the spring of 2016 at three georeferenced points located along the diagonal (28.28 m) of the plot. The soil samples were collected from 0-10 cm, 10-30 cm and 30-60 cm layers using a cylindrical core sampler (5.4 cm diameter), for the first two layers, and a gouge auger (2.5 cm diameter) for the deeper one.
The soil bulk density (BD) was determined for the first two layers using the samples, collected in the manner mentioned above and considering the volume of stones. For the third layer, the BD was determined by collecting the core samples from the soil mini-pit.
The soil samples were combined to form a composite sample of each layer, which represented the average soil characteristics of the plot. A total of 36 mineral soil samples were obtained. The list and abbreviations of all investigated parameters are shown in Table 1. The composite samples were air dried, sieved (2 mm mesh) and analyzed to determine the soil organic carbon (SOC) and total nitrogen (N tot ) content (Flash EA 1112 NCSoil, Thermo Fisher Scientific elemental analyzer, Pittsburgh, PA, USA), pH (soil to water ratio of 1:2.5), and particle-size distribution (by sieving and sedimentation; [32] and available P (P av ) [33,34]. CEC (extracted by barium chloride, pH 8.1), exchangeable cations (Ca, Mg, Na, and K) and BS were determined for the first layer only. Through the BD and stone volume, the SOC content was computed on an area basis (kg m −2 ) for mineral layers. The humus form was described in the field and then classified, according to [35]. The organic horizons (OL, OF and OH) were sampled using a 30 × 30 cm frame, collecting 27 composite samples of biomass for each horizon, which where oven-dried at 70 • C for 48 h, weighed and ground using a planetary ball mill. The SOC and N tot contents of organic layers were determined with a CN elemental analyzer.
To evaluate the level of accumulation of organic residues on the soil surface, the littering index (Litt ind ), obtained by the ratio between SOC Org and SOC of the first investigated layer, was calculated.
A quantitative assessment of forest humus forms was used to calculate the humus index (Hum ind ) [36,37]. The humus forms were classified and scaled in order of increasing accumulation of organic matter in the O horizons and decreasing burrowing activity in the A horizon, and the order was as follows: 1: Eumull, 2: Mesomull, 3: Oligomull, 4: Dysmull, 5: Hemimoder, 6: Eumoder, and 7: Dysmoder. Eumull, Oligomull and Eumoder were not present in our study sites. In the text and figures, the average value and standard error are used for the comparison of the different forest types.

Soil Biotic Characteristics: Microbial Activity and Microarthropod Communities
At each sampling plot, three core soil samples (100 cm 3 ) were collected in the 0-5 cm layer of the A horizon and were combined to form a composite sample and stored at 5-6 • C for microbial activity analyses. The potential soil respiration by titration [38,39] and bacteria count were determined (see Supplementary Materials Text S1).
Microarthropod samples were collected in May 2016, during a period of peak activity [40,41]. In each plot, after litter removal, three random soil samples were collected with a cylindrical corer (10.5 cm diameter), inserted at a depth of 10 cm. The samples were placed in black fabric bags, labeled, and processed in the laboratory on the day of collection to maximize organism survival and extraction. Soil fauna was extracted using Berlese-Tullgren funnels, following the standard QBS-ar methodology of [40], with an extraction duration of 10 days. Microarthropods were mainly determined on the order level. The community was characterized using: (a) Individual abundance, calculated as the mean of the three samples collected in each plot, for the more representative taxa (i.e., with at least 40 individuals); (b) richness within the plot, in terms of the number of euedaphic taxa, i.e., those well-adapted to the soil environment; and (c) the QBS-ar index, according to [40]. QBS-ar is an index based on the biological form approach and assesses the adaptation of organisms to an edaphic habitat through an index (eco-morphological index-EMI) that ranges from 1 to 20, with higher values assigned to microarthropod taxa that are better adapted to the soil environment. QBS-ar is then the sum of the EMI scores of microarthropod groups recorded in the sample, following the concept that the higher the number of groups that are well-adapted to the edaphic habitat, the higher the soil quality.

Plant Communities
A floristic survey was carried out in each 20 × 20 m plot. For each species (trees, grasses, forbs and moss layer), the percentage cover was visually estimated (expert-based, subdividing the whole plot into 4 sub-plots of 10 × 10 m to reduce the error). In order to determine if the three IAPS affect the native forest herb species (i.e., typical of the woodland habitat), the number of individuals of the forest herb species, Polygonatum multiflorum (L.) All., was recorded. This last species was chosen as an indicator species, since it was present in all forest types, and it is typical of the Central Europe oak or oak-hornbeam forests of the Carpinion betuli alliance, association Polygonato multiflori-Quercetum roboris Sartori 1980. The α-diversity was measured as the mean number of species registered in the plots and using the Shannon H' diversity index for each woodland type. The β-diversity was described using two different metrics: (a) The Jaccard's coefficient, a similarity measure based on co-occurrence of species between plots; and (b) the Whittaker index based on the ratio between αand γ-diversity. The γ-diversity, according to Whittaker, was calculated as the total number of species recorded in each woodland type (Table 1). In each plot, the number and cover of alien species were recorded, excluding the studied dominant alien species in each kind of woodland.

Data Analysis
An ANOVA model was fitted to assess differences in the means of abiotic and biotic factors among four woodland types. When necessary, the data were log-transformed to normalization. In particular, the Duncan post-hoc test was performed to assess the pairwise differences at a significance level of p < 0.05.
In order to analyze the strength of the ANOVA design and reduce the risk of Type II error (retaining a null hypothesis, when it was false), a post-hoc power analysis was carried out with the Gpower software [42]. The analyses showed a power higher than the critical value of 0.8 (or a very similar value of 0.79 for both SOCorg and SOCol), which is sufficient to detect significant effect sizes.
The ecological determinants of the structure of the different woodland types were investigated by a redundancy analysis (RDA). The RDA allowed the comparison the different woodlands, in terms of plant and microarthropod communities (composition and abundance), and the identification of the gradient trajectories for the different biotic and abiotic factors. The Monte Carlo permutation test (based on 999 iterations) was performed in order to assess both the significance of the environmental variables and the ordination axes. The Spearman's rho correlation analysis was performed to reduce the variables and noise in the variance components used for ordination. One of each pair of highly correlated/redundant variables (i.e., correlation coefficient ≥ 0.75) were removed from the analysis [43]. The inflation factor of the remnant variables was checked in order to detect for possible multicollinearity in the predictors. The outcomes are shown in a correlation matrix, the coefficient and significance of which are provided in Supplementary Materials Table S1.
The software employed for the RDA was CANOCO version 4.5. The univariate (ANOVA) and correlation analysis was performed in the R environment, version 3.3.1 [44]. The Levene test and Shapiro-Wilk test (car package in R; [45]) were conducted to test the homogeneity of the variance and normality assumption, respectively. The correlation plots were performed using the "corrplot" package in R [46].

Humus Forms
The main humus forms found in the study plots were, in order of decreasing biological activity, Mull and Moder ( Figure 1). The QF and RP plots were characterized by Mull forms: Dysmull for QF, Mesomull and Dysmull for RP (Humind equal to 4.0 and 2.7 ± 0.7 for QF and RP, respectively; Figure 2a). Dysmoder was the main form found in the QR plots, except for one plot, where Dysmull was present (average Humind of 6.0 ± 1.0). PS plots were characterized by the Hemimoder form (Humind of 5).
The higher SOCOrg (Figure 2b) was found in the Moder forms of QR and PS. The average SOCOrg at QR was significantly higher than the average SOCOrg at QF and RP (F3,8 = 5.12; p = 0.025 and p = 0.009, respectively), while non-significant differences were found with PS. The lower SOCOrg was found in the RP plots, and it was comparable with values at the QF plots and significantly different from QR and PS (F3,8 = 5.12; p = 0.077).
The OH horizons were found only at QR and PS, characterized by an average SOCorg of 0.19 ± 0.15 kg m −2 and 0.41 ± 0.31 kg m −2 in PS and QR, respectively.
The differences between the woodlands also concerned the type of organic matter, as shown by the CN values (Figure 2c-d). In particular, at QR, the average CNOL was statistically higher (F3,8 = 4.48; p = 0.0399) than that of all the other woodland types. CN was lower in the OF horizon. The values at QF and QR were statistically higher than CNOF at RP (F3,8 = 5.017; p = 0.026 and 0.011, respectively), and the CNOF of QR was also marginally different from that at PS (F3,8 = 5.017; p = 0.055). The QF and RP plots were characterized by Mull forms: Dysmull for QF, Mesomull and Dysmull for RP (Hum ind equal to 4.0 and 2.7 ± 0.7 for QF and RP, respectively; Figure 2a). Dysmoder was the main form found in the QR plots, except for one plot, where Dysmull was present (average Hum ind of 6.0 ± 1.0). PS plots were characterized by the Hemimoder form (Hum ind of 5).
The higher SOC Org (Figure 2b) was found in the Moder forms of QR and PS. The average SOC Org at QR was significantly higher than the average SOC Org at QF and RP (F 3,8 = 5.12; p = 0.025 and p = 0.009, respectively), while non-significant differences were found with PS. The lower SOC Org was found in the RP plots, and it was comparable with values at the QF plots and significantly different from QR and PS (F 3,8 = 5.12; p = 0.077).
The OH horizons were found only at QR and PS, characterized by an average SOC org of 0.19 ± 0.15 kg m −2 and 0.41 ± 0.31 kg m −2 in PS and QR, respectively.
The differences between the woodlands also concerned the type of organic matter, as shown by the CN values (Figure 2c-d). In particular, at QR, the average CN OL was statistically higher (F 3,8 = 4.48; p = 0.0399) than that of all the other woodland types. CN was lower in the OF horizon. The values at QF and QR were statistically higher than CN OF at RP (F 3,8 = 5.017; p = 0.026 and 0.011, respectively), and the CN OF of QR was also marginally different from that at PS (F 3,8 = 5.017; p = 0.055).  Figure S2. For abbreviations, see Table 1.

Mineral Soil
The soils almost always showed an A horizon of approximately 5 cm thick, followed by a AB or BA transition horizon up to 20-25 cm and often by a Bw up to 30-60 cm; below, a C horizon was present. From the taxonomic point of view (IUSS Working Group WRB, 2015), the first two horizons together reached the properties of an umbric diagnostic horizon, while the B horizon was often a cambic horizon.
The soils of investigated plots were Cambisols and Umbrisols. They are all moderately developed soils, with low BS except for RP stands, sandy loam soil texture, high rock fragment content and OC-rich topsoil.
Considering the whole investigated thickness, the soils showed a SOCmin with no statistical difference among woodland types. Approximately one third of SOC was in the first investigated layer (Supplementary Material Figure S2,b-d).
The main statistically significant differences (p < 0.05) in soil properties among woodland types concerned pHmin, CNmin, BS, CEC and Mg and Ca cation content (Figure 2). The pHmin was significantly lower at QF, QR and PS than at RP.
QF and QR showed similar CNmin values, higher than those at PS and RP (Figure 2e) F3,8 = 7.721, p < 0.05). All the forest types were characterized by very low BS (<35%) except for RP (35-50%) that showed higher values compared to the other forest stands (Figure 2h). The CEC showed medium values in all the forest types and higher values at PS (Figure 2m). In this woodland type, exchangeable Ca and Mg were higher than at QF and QR, but not significantly at RP for Ca (Figure 2f-g).
The exchangeable K content was similar for all the woodland types averaging to 0.2 ± 0.03 cmol(+) kg −1 (Supplementary Material Figure S2,g).  Figure S2. For abbreviations, see Table 1.

Mineral Soil
The soils almost always showed an A horizon of approximately 5 cm thick, followed by a AB or BA transition horizon up to 20-25 cm and often by a Bw up to 30-60 cm; below, a C horizon was present. From the taxonomic point of view (IUSS Working Group WRB, 2015), the first two horizons together reached the properties of an umbric diagnostic horizon, while the B horizon was often a cambic horizon.
The soils of investigated plots were Cambisols and Umbrisols. They are all moderately developed soils, with low BS except for RP stands, sandy loam soil texture, high rock fragment content and OC-rich topsoil.
Considering the whole investigated thickness, the soils showed a SOCmin with no statistical difference among woodland types. Approximately one third of SOC was in the first investigated layer (Supplementary Materials Figure S2b-d).
The main statistically significant differences (p < 0.05) in soil properties among woodland types concerned pHmin, CNmin, BS, CEC and Mg and Ca cation content (Figure 2). The pHmin was significantly lower at QF, QR and PS than at RP.
QF and QR showed similar CNmin values, higher than those at PS and RP (Figure 2e) F 3,8 = 7.721, p < 0.05). All the forest types were characterized by very low BS (<35%) except for RP (35-50%) that showed higher values compared to the other forest stands (Figure 2h). The CEC showed medium values in all the forest types and higher values at PS (Figure 2i). In this woodland type, exchangeable Ca and Mg were higher than at QF and QR, but not significantly at RP for Ca (Figure 2f-g).
The exchangeable K content was similar for all the woodland types averaging to 0.2 ± 0.03 cmol(+) kg −1 (Supplementary Materials Figure S2g).

Soil Biotic Characteristics
From the microbiological point of view, no significant differences were found among the forest types for both the soil respiration and bacteria count (Supplementary Materials Figure S2).
With regard to microarthropod communities, more than 8100 organisms, belonging to 17 taxa, were extracted from the soil samples, collected in May 2016 (Supplementary Materials Table S2). The euedaphic arthropod diversity, expressed as the number of taxa that are well-adapted to the soil environment, was the highest in RP, with no difference with PS, and lower in QF and QR (F 3,8 = 21.67; p < 0.001; Figure 3a-c). The mean number of Coleoptera larvae was the highest in QR (F 3,8 = 4.921; p = 0.032). No differences were detected in Collembola, Pauropoda, Protura and Symphyla abundance. However, it should be noted that, for the euedaphic taxa of Protura and Symphyla, QR resulted in no individuals at all and in a mean that was considerably lower than the other woodland types (Supplementary Materials Table S2). The soil quality, evaluated by the QBS-ar index, was higher in RP, compared with the other forest types (F 3,8 = 5.322; p = 0.026).

Soil Biotic Characteristics
From the microbiological point of view, no significant differences were found among the forest types for both the soil respiration and bacteria count (Supplementary Material Figure S2).
With regard to microarthropod communities, more than 8100 organisms, belonging to 17 taxa, were extracted from the soil samples, collected in May 2016 (Supplementary Materials Table S2). The euedaphic arthropod diversity, expressed as the number of taxa that are well-adapted to the soil environment, was the highest in RP, with no difference with PS, and lower in QF and QR (F3,8 = 21.67; p < 0.001; Figure 3a-c). The mean number of Coleoptera larvae was the highest in QR (F3,8 = 4.921; p = 0.032). No differences were detected in Collembola, Pauropoda, Protura and Symphyla abundance. However, it should be noted that, for the euedaphic taxa of Protura and Symphyla, QR resulted in no individuals at all and in a mean that was considerably lower than the other woodland types (Supplementary Material Table S2). The soil quality, evaluated by the QBS-ar index, was higher in RP, compared with the other forest types (F3,8 = 5.322; p = 0.026).  Figure S2. For abbreviations, see Table 1.

Plant Community
This study identified 82 species which were recorded within the plots and surveyed during May-July 2016. The plant diversity was different among the studied woodlands (Figure 3d Figure S2).  Figure S2. For abbreviations, see Table 1.

Plant Community
This study identified 82 species which were recorded within the plots and surveyed during May-July 2016. The plant diversity was different among the studied woodlands (Figure 3d-g). The highest values of α-(richness and Shannon) and γdiversity were in RP, and the lowest were in QR (richness:  Figure S2).

Relationships between Biotic Communities and Ecological Determinants
The RDA analyses (Figure 4a,b) indicated that several ecological determinants were associated with the different woodland types in relation to the biotic communities (microarthropods and plants). The RDA ordinations resulted in high eigenvalues and cumulative percentage variances in the species-environment data, which are indicative of the distinctive species assemblages across the different woodland types (Table 2). Particularly, the biplots of the first three axes accounted for 57.7% and 82.6% of this variance for the microarthropod and plant communities, respectively.   Concerning the microarthropod assemblage (Figure 4a, Table 2a), QBS-ar and the Acari abundance accounted for the highest level of variance (19% each) and greatly varied across forest types (QBS-ar: F = 2.62; p = 0.001; Acari: F = 2.41; p = 0.005). In addition, the Symphyla abundance explained significant amounts of variance (0.11%; F = 1.87; p = 0.050). The main gradient exhibited positive scores of the above-mentioned variables for the RP stands and negative scores for the QR stands.
Concerning the plant assemblage (Figure 4b, Table 2b), the woodland types varied in CN OL (F = 4.89; p = 0.001; 33% of explained variance) and CN min (F = 2.38; p = 0.042). Both factors were positively associated with QR and negatively associated with RP stands. The concentration of K in the soil was significant (F = 2.95; p = 0.037) and was associated with PS stands. Among the biotic factors, the different woodlands greatly differed in the total number of plant species (γ-div: variance of 16%; F = 3.56; p = 0.001), and a positive association with RP stands was evident. Again, the main gradient exhibited positive and negative scores of the variables for the RP and QR stands, respectively.
A significant correlation between the selected variables was found in 46 out of 264 possible combinations, explaining the variation across the dataset ( Figure 5; Supplementary Materials Table S1). Among the other variables, some strong correlations (correlation coefficient > 0.7; and significance ≤ 0.01) were identified between ecosystem components ( Figure 5). For instance, CN min was negatively correlated with Alien richness; CN min and CN OL were negatively correlated with Eu_rich (Euedaphic arthropod taxa α-diversity); and QBS-ar was positively correlated with the pH value and negatively correlated with CN OL . Concerning the plant assemblage (Figure 4b, Table 2b), the woodland types varied in CNOL (F = 4.89; p = 0.001; 33% of explained variance) and CNmin (F = 2.38; p = 0.042). Both factors were positively associated with QR and negatively associated with RP stands. The concentration of K in the soil was significant (F = 2.95; p = 0.037) and was associated with PS stands. Among the biotic factors, the different woodlands greatly differed in the total number of plant species (γ-div: variance of 16%; F = 3.56; p = 0.001), and a positive association with RP stands was evident. Again, the main gradient exhibited positive and negative scores of the variables for the RP and QR stands, respectively.
A significant correlation between the selected variables was found in 46 out of 264 possible combinations, explaining the variation across the dataset ( Figure 5; Supplementary Material Table  S1). Among the other variables, some strong correlations (correlation coefficient > 0.7; and significance ≤ 0.01) were identified between ecosystem components ( Figure 5). For instance, CNmin was negatively correlated with Alien richness; CNmin and CNOL were negatively correlated with Eu_rich (Euedaphic arthropod taxa α-diversity); and QBS-ar was positively correlated with the pH value and negatively correlated with CNOL.  Table S1. White stars indicate correlation coefficients > 0.7 and a significance level p ≤ 0.01. For abbreviations of the variables, see Table 1.

Discussion
This work examined how Q. rubra, P. serotina and R. pseudoacacia affected soil properties and ecosystem components of native forests dominated by the trees, Q. robur and C. betulus. The three IAPS are associated with a series of abiotic and biotic ecological alterations, affecting, in different ways and degree, the resident forest ecosystem.
While it is acknowledged that the sample size in this study is low, and to some extent, the data analysis could be affected by the effect of the individual sites, the authors underline that this sampling  Table S1. White stars indicate correlation coefficients > 0.7 and a significance level p ≤ 0.01. For abbreviations of the variables, see Table 1.

Discussion
This work examined how Q. rubra, P. serotina and R. pseudoacacia affected soil properties and ecosystem components of native forests dominated by the trees, Q. robur and C. betulus. The three IAPS are associated with a series of abiotic and biotic ecological alterations, affecting, in different ways and degree, the resident forest ecosystem.
While it is acknowledged that the sample size in this study is low, and to some extent, the data analysis could be affected by the effect of the individual sites, the authors underline that this sampling strategy considered only those forest sites falling in the same pedological landscape unit (i.e., with a particular geological origin, landform, topography, hydrography, historical land use and soil type; see: [47]). Thus, the studied sites shared very similar chemical-physical soil characteristics until vegetation replacement, which allowed the containment of the site effect and the investment of sampling efforts in comparing different ecosystem components in terms of a greater number of invaded forest types.

Soil Abiotic Characteristics: Humus Forms
The major changes in soil characteristics mainly involved the organic layers. The substitution of native forests by P. serotina and Q. rubra resulted in a shift toward a less active humus form, from Dysmull to Hemimoder and Dysmoder, respectively, as a consequence of organic matter that did not decompose as easily [48]. Conversely, RP produced litter that rapidly decomposed, resulting in the more active Mesomull form, compared to native forests. Leaf litter C:N and polyphenolic compounds, such as tannins, were linked to the rate of leaf litter decomposition [49,50]. The action of tannins, whose content in the leaves of red oak and black cherry is higher than in common oak [51], depends on the species. Further, [52] observed that the litter tannins of Q. rubra are important for lowering the rate of soil-nutrient cycling to a greater extent than the tannins of other species.
Litter quantity and quality influence the litter decomposition rate and the dynamics of nutrient mineralization and immobilization [53], which, in turn, may affect the microbial decomposer and detritivore communities [54] that are directly involved in litter decomposition. The observed sequences and CN values across organic horizons supported the different turnover of organic matter among woodlands, showing the capacity of the organic layers of QR stands to act as a carbon sink.

Soil Abiotic Factors: Mineral Soil Properties
The forest substitutions also impacted the mineral soil, although to a lesser extent than the organic layers. Forest types affected the soil organic matter quality, as shown by the lower soil CN min in PS and RP (absolute lowest value), compared to QF and QR stands. Rice et al. (2004) demonstrated that due to an elevated organic matter mineralization rate (low C:N) in the topsoil, R. pseudoacacia stands had a higher nitrogen concentration than stands dominated by other species, such as pine and oak. The organic matter mineralization in R. pseudoacacia is also favored by lower levels of canopy cover and, consequently, a higher soil temperature compared to the conditions under a very dense canopy, for instance, those of QR. Previous studies demonstrated that P. serotina can modify N and C cycles, increasing carbon turnover via a fast decomposition of litter and humus formation, accelerating nutrient circulation at the ecosystem level in comparison to native forest stands [18,51,55]. In this study, such effects of PS stands were less evident, and they generally showed an intermediate behavior between QR and QF.
The impact on the chemical properties of the topsoil exhibited species-specific patterns so that increasing contents of exchangeable Ca, Mg and K, as well as higher CEC or BS values, were mostly associated with RP and PS. RP stands exhibited a higher level of Ca than QF stands, as already observed by [56]. Similarly, [57] reported that R. pseudoacacia soils had elevated Ca levels. In the PS stands, the Ca content in the soil was similar to that in the RP stands and was much higher than the content observed by [18]. The soils at the QF and QR sites had a similar exchangeable Ca cation content, likely indicating similar pedoclimatic conditions. This pattern is in accordance with previous studies [19,50].
The authors claim that the similarities and differences in soil nutrients across different studies have to be interpreted in consideration of the confounding factors that can be present in the soil exchange complex. For instance, since the availability of some plant nutrients is greatly affected by soil pH [58], the variation in soil pH may explain the differences between the nutrient/micronutrient content presented in this study and that detected by [18,50]. For such reasons, the variations in soil elements across sites, invaded and uninvaded by IAPS, should be assessed and interpreted with caution.

Biotic Factors: Microarthropod and Plant Communities
Concerning the biotic components of the ecosystem, our results highlighted that RP and QR exhibited (as a general rule) the highest and the smallest biodiversity values, respectively, while QF and PS showed intermediate values. On the whole, the dominance of Q. rubra in forest ecosystems was found to negatively impact the biotic components of native forests, such as Acari densities pH [59], micromycetes and ammonifying microorganisms [60], and plant communities pH [21,61]. In this study, Q. rubra stands were characterized by a low-quality litter layer (e.g., a high C:N ratio) and low cover of a few herbaceous species that have likely negatively influenced the microarthropod community, as shown by the low number of eudaphic taxa, such as Acari, Protura and Symphyla. Indeed, litter quality is considered one of the most important factors determining soil fauna assemblages [62]. Recalcitrant litter requires a peculiar decomposer system [63]. Moreover, the simplification of the vegetation structure and composition could affect groups of herbivores and predators [64].
In our results, the microarthropod soil biological quality index, QBS-ar, showed the highest values in RP forest soil (with no differences from that in PS). This result is in contrast with the findings of [56], who reported a higher QBS-ar value in native Quercus stands than in stands invaded by R. pseudoacacia. This discrepancy can be due to differences in: (1) The biogeographical and pedological characteristics of the two study areas (i.e., Mediteranean and Continental regions); (2) past forest management of the sampling sites (e.g., logging and coppicing); (3) dominant native trees that, in [56], consisted of different Quercus species (Q. cerris L. and Q. pubescens Willd.). It is known how the differences in environmental conditions across the sampling sites can influence the arthropods abundance [64]. Indeed, most studies reported the negative effects of invasive species on arthropod communities [65]. However, research on detritivores showed an increase, or the stability, of this functional group in forests invaded by alien trees [64], which is probably due to the higher quality and larger amount of the litter [62,66]. Similarly, our RP stands hosted a high number of euedaphic taxa that feed on plant material and detritus, fungi and bacteria, such as Acari, Collembola, Protura, Symphila and Pauropoda, which may not have been negatively influenced by the presence of the IAPS.
Our results concerning the soil fauna in the different woodland types are related to the litter quality. However, as responses of microarthropods groups to the soil factors complex are very variable and poorly understood [64,67], further investigations are certainly needed.
With regard to the total number of herb forest species (richness and Shannon), this study (RDA and ANOVA results) is, as a general rule, in accordance with previous findings. Specifically, Q. rubra has a negative effect on the number and cover of herbaceous species, which is related to the slow decomposition of the litter layer [61]. Similarly, P. serotina has an overall negative effect on the structure and composition of native plant communities [15]. However, in our results (of RDA) and in previous studies, it has also been associated with forest stands that are rich in native species [68] (i.e., C. majalis, P. multiflorum, V. minor, etc.). However, in this study, no difference in typical forest herb species were detected across different wood types, probably due to the ability of perennial species (with rhizomes and bulbs) to resist disturbance and soil modifications.
On the other hand, R. pseudacacia has been reported to have a positive effect on the cover of generalist (nitrophilous, ruderal and alien) and graminoid species, which seem to slow down forest succession and regeneration of native trees [68,69], and this is consistent with our results. Indeed, the number of native species noted under the R. pseudoacacia canopy may be even higher than in natural forests. However, the species composition may differ significantly from the species list of natural forest communities [70].
Across woodland types, no differences or trends were detected with regard to bacterial activity (abundance and soil respiration), which, generally, can vary among sites dominated by different IAPS [71]. Further analysis of the composition of bacteria communities may better explain these first observations.

Conclusions
The investigations regarding the impacts of IAPS on ecosystem components are often only based on a single invasive species. This comparative study of three alien species (P. serotina, Q. rubra, and R. pseudoacacia), listed among the worst IAPS, improves our understanding of the changes in ecological factors that accompany the establishment of alien species. These IAPS can impact forest ecosystems in terms of humus forms, carbon and nutrient cycles, and species composition and abundance (plant and microarthropod communities). Particularly, Q. rubra was associated with the major negative impacts on organic layers and exhibited low/modified levels of microarthropod and plant biodiversity. In the long term, the invasion of the studied species, especially P. serotina and Q. rubra, are expected to further weaken forest ecosystems, which are already under pressure due to global change and forest fragmentation. Particularly, many native species may be replaced or lost due to the invasion of these alien species. Indeed, across Europe, both species have been reported to have greatly replaced native forests.
Future research on the impact of IAPS on forest ecosystems in highly anthropized regions should include multiple stressors, such as forest management and the degree of anthropic pressure. This could be useful in implementing appropriate and targeted management strategies to combat IAPS.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/10/10/842/s1. Author Contributions: R.G., C.F., R.C., E.C., and S.C. conceived the idea and the sampling design. R.G., C.F., E.C. led the writing of the manuscript and of the statistical analyses. R.G., C.F., R.C., E.C., C.M., and G.B. participated in the field data collection for the different abiotic and biotic ecosystem components. All authors critically contributed to the paper and approved the final version for publication.
Funding: This research was partially funded by the LIFE financial instrument of the European Commission to Lombardy Region LIFE14 IPE IT 018GESTIRE2020-Nature Integrated Management in 2020.