Beijing is an international metropolis located in the North China Plain, with a history of over 1,000 years. In the past four decades, it has spread out rapidly with concentric ring roads, during which time its population has reached over 21 million. It has a temperate, semi-humid, continental monsoon climate with a mean annual temperature of 11–12°C and mean annual precipitation of 500 mm.

The study sites were located along a south–north transect that runs through the urban center (Figure 1) and covers both the fastest-growing urban regions and rural areas of Beijing (Peng et al., 2016). Along this transect, 64 Chinese pine stands located on 24 roadsides, in 25 parks, and in 15 neighborhoods were selected (Figure 1 and Table 2). A roadside is a strip of plants, mainly located between a roadway (carriageway) and a sidewalk (pavement). A park is an area of open space provided for recreational use, usually owned and maintained by a local government. A neighborhood is a district, especially one forming a community within a town or city and sharing the same services and management. The distance of the chosen sites from the city center ranged from 1.00 to 24.89 km. Study sites were sampled at 1-km intervals to make the distribution of sites spatially balanced.

Chinese pine (Pinus tabuliformis Carr.) is a monoecious, wind-pollinated, and predominantly outcrossing species (Wang et al., 2010). It is a light-loving tree species with low water consumption and slow growth but strong drought tolerance (Zheng and Fu, 1978). Because of its high adaptivity to poor environments, Chinese pine has been widely planted in urban areas of Beijing for a long time.

Between July and August 2016, three well-developed pines at each site (a total of 64 sites and 192 trees) were selected to examine their traits. The age of each tree was estimated by counting the number of whorls of branches along the main trunk because pines typically grow in annual spurts, putting on one whorl of branches each year (Hartman et al., 2009; Rollinson, 2012). The average age of the trees chosen was 22 years, their average height was 5.85 m, and average diameter at breast height was 16.62 cm. To determine leaf functional traits, 80–100 current- and previous-year needles were collected randomly from each selected tree; these were fully expanded, sun-exposed, in good health, and at a height of approximately 1.8 m (Cornelissen et al., 2003).

Stomatal density was determined using scanning electron microscopy in accordance with the methods of Yang et al. (2017). Three fresh needles wiped with a wet towel were placed immediately in a plastic bottle with a fixative solution (70%, v/v, ethanol:glacial acetic acid:formalin; 90:5:5) to maintain the integrity of cellular structure. Before examination with a scanning electron microscope (SEM; FE-SEM-EDS, SU-8020, Hitachi, Tokyo, Japan), the middle segment of each needle, approximately 1 cm long, was dehydrated with ethanol, incubated overnight in isoamyl acetate, and then dried to the critical point. The number of stomata within the images was counted manually and then converted to stomatal density (number of stomata per square-millimeter of adaxial leaf surface). Images of needles captured by the scanning electron microscope were also used to measure leaf width using ImageJ software (US National Institutes of Health, Bethesda, MD, United States). Leaf length of three fresh needles was measured using a digital caliper. Twelve fresh needles were scanned using a digital scanner (HP Scanjet G3110, Hewlett-Packard Development Company, Beijing, China) to measure leaf area using ImageJ software. After leaf area measurement, needles were dried at 80°C to a constant weight. SLA was calculated as the leaf area divided by its oven-dried mass (Cornelissen et al., 2003). The remaining needles were dried and ground to measure leaf nutrient concentrations. Leaf N concentration was measured using an automatic elemental analyzer (Vario EL III, Elementar, Germany), and leaf P and K concentrations were measured using Inductively Coupled Plasma Optical Emission Spectroscopy (ICP-OES; Prodigy, Leeman, Hudson, NH, United States) after digestion with freshly distilled concentrated HNO₃/H₂O₂ in a microwave oven (Oliva et al., 2003).

All indicators were tested for normality using the Shapiro–Wilk test and were log-transformed before analysis if they were distributed non-normally. Levene’s test and Q-Q plots of residuals were used to test the homogeneity of variance. The Mantel test was performed to test the spatial autocorrelation using the “mantel.rtest” function in the “ade4” package. Among the traits studied, only stomatal density showed significant spatial autocorrelation (Supplementary Table 1). Analysis of variance was used to evaluate leaf age, land-use type, and the effects of their interactions on traits. Analysis of variance was performed using the “Anova” function with type II tests in the “car” package to control unequal sample size (Lewsey et al., 1997; Langsrud, 2003). The “emmeans” function in the “emmeans” package, which provides the t-value and corrected p-values automatically for unequal sample size, was executed for multiple comparisons of trait indicators (Russell, 2018).

For traits without significant spatial autocorrelation, to test changes in traits along the urban–rural gradient, a generalized linear mixed model (GLMM) was fitted with traits as the response variable, distance of the site from the city center as the fixed effect, and land-use type as the random intercept (the gamma distribution with LOG LINK function) in R software. Likelihood ratio tests of a full model against a null model were used to measure the significance of the fixed effect (Fajardo and Siefert, 2018). The explanatory power of the model was evaluated by calculating the conditional R² (Nakagawa and Schielzeth, 2013). For stomatal density with significant spatial autocorrelation, we performed generalized linear mixed models via Penalized Quasi-Likelihood estimation (GLMMPQL) using “glmmPQL” function in “MASS” package to test changes in stomatal density along the urban–rural gradient. The GLMMPQL enables the building of spatial models with dependent data not normally distributed (Dormann et al., 2007). The response variable, fixed effect, random intercept, the error distribution and link function in the GLMMPQL were set to be the same as in the GLMM. Taking into account the unequal sample size, these generalized linear mixed models were also performed in SAS software (PROC GLIMMIX procedure, METHOD = RMPL, DDFM = KENWARDROGER) followed by a type II sum of squares ANOVA (Langsrud, 2003; Spilke et al., 2005). To test whether changes in traits along the urban–rural gradient varied with leaf age and land-use type, general linear models with traits as the response variable and distance as the predictor variable were fitted for each leaf age class and land-use type.

Variation partitioning analysis was performed using the ‘‘varpart’’ function in the ‘‘vegan’’ package to compare the percentage of trait variations explained by leaf age, land-use type, urban–rural gradient, and their intersections. Statistical analyses were performed using the statistical software R v. 3.6.2¹ and SAS v. 9.4 (SAS Institute, Cary, NC, United States).

The mean trait values of current- and previous-year needles were 65.37 and 59.61 n/mm² for stomatal density, 13.73 and 14.25 cm for leaf length, 0.25 and 0.27 cm for leaf width, 2.09 and 2.15 cm² for leaf area, 57.45 and 48.34 cm²/g for SLA, 15.69 and 14.48 mg/g for leaf N concentration, 1.24 and 0.81 mg/g for leaf P concentration, and 7.38 and 4.48 mg/g for leaf K concentration, respectively.

Stomatal density, leaf width, SLA, and leaf N, P, and K concentrations differed significantly between leaves of different ages (Figure 2). Compared with previous-year needles, current-year needles had 10% greater stomatal density, 8% smaller leaf width, 20% greater SLA, 8% higher leaf N concentration, 53% higher leaf P concentration, and 51% higher leaf K concentration (p < 0.05, Figure 2).

The mean trait values of plants on roadsides and in parks and neighborhoods were 62.65, 67.90, and 65.50 n/mm² for stomatal density, 12.55, 14.31, and 14.65 cm for leaf length, 0.22, 0.25, and 0.28 cm for leaf width, 1.87, 2.22, and 2.21 cm² for leaf area, 59.81, 57.09, and 54.28 cm²/g for SLA, 15.36, 15.80, and 16.02 mg/g for leaf N concentration, 1.33, 1.19, and 1.18 mg/g for leaf P concentration, and 7.66, 6.83, and 7.83 mg/g for leaf K concentration, respectively.

Stomatal density, leaf length, width, and area, and leaf P and K concentrations varied significantly with land-use types (Figure 2). Compared with plants on roadsides, those in parks and neighborhoods had 9 and 9% shorter leaf length, 18 and 21% smaller leaf width, and 13 and 12% smaller leaf area, but 7 and 9% higher leaf P concentration, respectively (p < 0.05, Figure 2). In addition, compared with plants in parks, those in neighborhoods had 11% smaller leaf width and 13% lower leaf K concentration, while the values of other traits showed no significant differences between plants in parks and those in neighborhoods (Figure 2).

Stomatal density, leaf width, and leaf P concentration increased significantly along the urban–rural gradient, while leaf K concentration decreased significantly (p < 0.05, Table 3 and Figure 3).

Importantly, urban–rural gradients of leaf functional traits varied with leaf age and land-use type (Figure 3). Urban–rural gradients of stomatal density, leaf width, and leaf P concentration were identified only for current-year needles, and gradients in leaf K concentration were only observed for previous-year needles. Urban–rural gradients in stomatal density and leaf K concentration were found on roadsides, gradients in leaf width and leaf P and K concentrations were found in parks, while only gradients in leaf P concentration were found in neighborhoods.

The results of variation partitioning analysis revealed that among leaf age, land-use type, and urban-rural gradient, leaf age explained the greatest variations in leaf nutrient traits (leaf N, P, and K concentrations), while land-use type explained the greatest variations in leaf morphological traits (stomatal density, and leaf length, width and area) (Figure 4). This suggests that leaf age and land-use type were the most important factors influencing leaf nutrient and morphological traits in urban environments, respectively.

Leaf functional traits are associated with plant resource acquisition, damage resistance, and mechanical strength (Wang et al., 2016; Blonder et al., 2020). Our results implied that previous-year needles tended to have greater mechanical strength (smaller SLA), but less production capacity (lower stomatal density, and lower N, P, and K concentrations) than current-year needles (Zheng et al., 2017; Bucher et al., 2019; Liu Z. et al., 2020).

Leaf functional traits changed with age, which is consistent with results reported previously. For example, as leaf age increases, SLA decreases in another Pinus species (Pinus sylvestris var. mongolica; Zheng et al., 2017; Liu X. et al., 2020), non-pine woody species (Römermann et al., 2016), and herbaceous species (Bucher et al., 2019). Stomatal density decreases with age in Tsuga heterophylla (Kouwenberg et al., 2004) and Trichosanthes cucumerina (Adebooye et al., 2012). Leaf nutrient concentrations also decrease with age in Pinus species (Pinus koraiensis; Liu X. et al., 2020; Liu Z. et al., 2020) and other tree species (Fraxinus mandshurica and Acer mono; Liu et al., 2019).

More organic matter might accumulate in leaves as they age, explaining the smaller SLA for previous-year needles than current-year needles (Liu Z. et al., 2020). A potential explanation for the decrease in stomatal density as leaf age increases may be that leaf growth and expansion widen the spaces between stomata, thereby leading to fewer stomata per given leaf area (Adebooye et al., 2012). Decreases in leaf nutrient concentrations as leaves grow may result from the transport of nutrients from old leaves to young leaves (Zheng et al., 2017; Liu et al., 2019) to maintain favorable nutrient status in young leaves, thus maintaining optimal metabolic activities and maximizing plant growth (Schreeg et al., 2014).

Leaf functional traits varied with land-use types (Figure 2). Plants growing on roadsides had traits (lower levels of stomatal density, and smaller leaf length, width and area) representative of a reduced capacity to acquire resources (e.g., light and carbon dioxide) compared with plants growing in parks and neighborhoods (Figure 2). This decrease in resource acquisition capacity is detrimental to plant growth and production potential (Fanourakis et al., 2015; Wright et al., 2017). Our results therefore imply that roadside environments do not favor the performance of plant functions (de Bello et al., 2010; Bussotti and Pollastrini, 2015; Fanourakis et al., 2015). Similarly, compared with trees of the same species growing on roadsides, the deciduous, broad-leaved tree species Zelkova serrata (Thumb.) Makino growing in parks in Tokyo, Japan exhibits lower levels of stomatal density (Osone et al., 2014), and the evergreen, broad-leaved tree species Quercus ilex growing in parks in Naples, Italy exhibits smaller leaves (Maisto et al., 2013).

Different levels of leaf functional traits among land-use types are related to each type’s unique environment. Stomata are important channels for plant water and gas exchange. The potential explanation can be proposed for the lower levels of stomatal density in plants on roadsides than in parks and neighborhoods. Drier and hotter environmental conditions on roadsides than in parks and neighborhoods may cause lower levels of stomatal density in plants (Fraser et al., 2009; Lau et al., 2018; Lu et al., 2019). Bao et al. (2001) reported that the relative humidity of air on roadsides was 12 and 7% lower than that in parks and neighborhoods, respectively, while the air temperature on roadsides was 1.49 and 0.95°C higher than that in parks and neighborhoods, respectively. Lau et al. (2018) and Lu et al. (2019) even revealed the molecular mechanism of reduced stomatal density caused by high temperature, during which expression of SPEECHLESS, a basic-helix-loop-helix transcription factor, is repressed. Deformation and curling of leaves in a dry and hot environment may cause the surface of epidermal cells to twist and fold, and the stoma to become trapped in mesophyll cells, leading to a further decrease in stomatal density. In addition, damage to leaves may inhibit stomata formation, reducing the number of stomata and finally resulting in a decrease in stomatal density.

Leaf size determines plants’ ability to capture light resources and carbon directly, and is affected by a variety of environmental factors, such as temperature, precipitation, and CO₂ concentration (Mcdonald et al., 2003; Wright et al., 2017; Wang et al., 2018). Plants tend to have smaller leaves when environmental conditions are stressful, and those growing on roadsides had smaller leaves than those growing in parks and neighborhoods (Figure 2), which may be attributable to the relatively drier and hotter environmental conditions on roadsides (Bao et al., 2001; Maisto et al., 2013; Wright et al., 2017). High levels of pollutants may be another reason for reduced leaf size on roadsides; because of high traffic activity, roads are often more polluted than parks and neighborhoods. Mukherjee and Agrawal (2018) reported that total suspended particulate matter, inhalable particulate matter, and SO₂ concentration levels on roadsides were 64, 96, and 202% higher than those in neighborhoods. Pollutants cannot only hinder gas exchange in photosynthesis, but also enter the plant directly and destroy photosynthetic tissue, leading to a decline in photosynthesis and subsequent reduction in leaf size (Bhatti and Iqbai, 1988; Leghari and Zaidi, 2013). Leghari and Zaidi (2013) found that the leaf length, width, and area of plants in polluted areas decreased by 20–73, 18–51, and 23–58%, respectively, compared with those of plants in unpolluted areas. Finally, lower soil nutrient levels may also contribute to reduced leaf size on roadsides (Mcdonald et al., 2003). The level of soil organic carbon, calcium, and magnesium contents in parks and neighborhoods was 36 and 62, 70 and 113, and 81 and 168% higher, respectively, than that on roadsides in Beijing (unpublished data).

Mean values of functional traits in plant communities have been used widely to determine the way the urban environment filters species into urban plant communities (Williams et al., 2015; Kalusova et al., 2017). In this study, we explored plant responses and adaptations to urban environments based on trait values measured in situ at the species level. We found that stomatal density, leaf width, and leaf P concentration increased gradually along the urban–rural gradient, while leaf K concentration decreased (p < 0.05); the other traits studied did not change significantly (p > 0.05, Figure 3 and Table 3).

Similarly, other studies have found that plants tend to have lower stomatal density and smaller leaves in urban areas than in rural areas (Pourkhabbaz et al., 2010; Leghari and Zaidi, 2013; Alotaibi et al., 2020). However, opposite results have been reported (Table 1). In addition, significant differences in SLA and leaf N concentration between urban and rural sites have been reported (Table 1).

The reasons for the inconsistences in urban-rural gradients of plant traits between different this study and other studies (Figure 3 and Table 1) may be as follows. First, the urban–rural gradients of traits depend upon the species studied. In this study, we found no significant urban–rural gradients for leaf length or area, SLA, or leaf N concentration in Chinese pine, an evergreen coniferous species; however, significant urban–rural gradients for these traits have been found in deciduous species, such as Populus tremula and Quercus robur (Nikula et al., 2010; Searle et al., 2011, 2012). This is consistent with a study by Jiahao et al. (2018), who found more moderate variations in stoichiometric and morphological traits of leaves of evergreen conifers than those of evergreen and deciduous broad-leaved species with changing environmental conditions. Differences in leaf N concentration may also be attributable to the fact that urban–rural gradients of leaf N concentration in broad-leaved forests are related closely to urban–rural gradients in aerial N deposition (Nikula et al., 2010; Searle et al., 2012); however, Chinese pine, an evergreen coniferous-leaved tree species, captures less deposition-derived N and acquires this N earlier than broad-leaved species (Rothe et al., 2002; Kristensen et al., 2004; Gundersen et al., 2009). Second, urban–rural gradients in traits depend upon the local climate. In New York, with a temperate climate, leaf area and leaf N concentration levels in Quercus robur are 730 and 20% higher, respectively, in urban areas compared with rural areas (Searle et al., 2012); meanwhile, in Potchefstroom, South Africa, with a savanna climate, leaf area and leaf N concentration are 35 and 33% lower, respectively, in urban areas (vanRensburg et al., 1997). Third, urban–rural gradients in traits vary inter-annually. Nikula et al. (2010) found that the mean leaf N concentration level of Populus tremula leaves was similar in 2007 and 2006, approximately 2.19 and 2.22%, respectively; however, the urban–rural difference in leaf N concentration in 2007 was approximately 33 times the urban–rural difference in 2006. Fourth, land-use type can enhance or weaken the degree of urban–rural gradient or even change its direction. The difference in leaf area of Ficus altissima between urban and rural sites when roadsides were used as urban sites was approximately twice that when neighborhoods were the urban sites (Alotaibi et al., 2020). Moreover, when industrial areas or roadsides were used as urban sites, the SLA of Ficus nitida in urban areas was approximately 5% lower or 74% higher, respectively, than that in rural areas (El-Khatib et al., 2020).

The occurrence of urban–rural gradients in leaf functional traits varies with leaf age (Figure 3). Functional traits of current-year needles were more sensitive to the urban–rural gradient than were those of previous-year needles (Figure 3). This may be attributable to the following: (1) the trait plasticity of evergreen woody species decreases with increasing leaf age (Niinemets, 2016); (2) old needles function primarily in mechanical protection, while young needles are the primary location of rapid photosynthesis, transpiration, and growth, which are influenced strongly by environmental changes (Tissue et al., 2001; Sariyildiz and Anderson, 2006); and (3) redistribution of nutrients to young needles may cause greater sensitivity of young needles to environmental changes (Turner and Olson, 1976; Yuan et al., 2018). Therefore, environmental changes along an urban–rural gradient, such as temperature fluctuations (Xu et al., 2001) and aerial ammonia pollution (Judžentienė et al., 2006), are more likely to alter traits in younger needles than in older needles.

The occurrence of urban–rural gradients in leaf functional traits also varies with land-use type (Figure 3), probably because the range of environmental factors along the urban–rural gradient varies with land-use type. There are more traits that show significant variation along an urban–rural gradient in parks and on roadsides than there are in neighborhoods, perhaps because neighborhoods have a more homogeneous environment attributable to similar landscape structures and management practices (Polsky et al., 2014; Wang et al., 2015). Urban–rural gradients in stomatal traits appear only on roadsides, which may be related to the urban–rural gradient in pollution levels (Bhatti and Iqbai, 1988; Leghari and Zaidi, 2013; Osone et al., 2014). There is more traffic pollution on roadsides than in parks and neighborhoods, and vehicle traffic is heavier on urban than on rural roads. Leaf width showed urban–rural gradient changes only in parks, which may be attributable to their greater temperature gradients. During the planning and development stages, a certain proportion of green spaces and impervious surfaces for roadsides and neighborhoods needs to be met. However, the proportion of land cover is more diverse in parks; for example, parks located in the city center often have more historic buildings (impervious surfaces) than those located in rural areas, which would widen urban–rural temperature gradients in parks.

High temperature, CO₂ and pollution, and low moisture content in cities may decrease stomatal density, leaf width, and leaf P concentration, and increase leaf K concentration in urban plants (Woodward and Kelly, 1995; Luomala et al., 2005; Fraser et al., 2009; Pourkhabbaz et al., 2010; Wright et al., 2017). Alterations in these traits may represent adaptive adjustment of plants to environmental changes. Reduction in stomatal density cannot only restrict pollutants and the entry of CO₂ into leaves (Woodward and Kelly, 1995; Verma et al., 2006; Pourkhabbaz et al., 2010), but also improve water use efficiency to prevent excessive water loss (Woodward and Kelly, 1995; Luomala et al., 2005; Wright et al., 2017). Narrowed leaves can evacuate heat quickly (Smith, 1978; Shin et al., 2000). Lower P concentration in leaves of urban plants can constrain high temperature-induced increases in metabolic reaction rates and growth rates to avoid rapid aging (Reich and Oleksyn, 2004; He et al., 2008; Pretzsch et al., 2017). In future studies, leaf functional traits and environmental factors (e.g., soil and air conditions) should be monitored simultaneously to reveal the mechanism of urban environmental influence on leaf functional traits.

Leaf age and land-use type were the most important factors influencing leaf nutrient and morphological traits, respectively, in urban environments (Figure 4). Leaf nutrients are transferred from old leaves to new leaves to prevent nutrient loss as old leaves wither (Zheng et al., 2017; Liu et al., 2019). Therefore, leaf age has a greater impact on leaf nutrients than land-use type or urban–rural gradient. Leaf morphological traits were more sensitive to land-use type than leaf age or urban–rural gradient. This might be because (1) the contribution of ontogeny to variations in leaf size and stomatal density are relatively small once leaves mature, and (2) changes in environmental conditions across land-use types are greater than those along the urban–rural gradient. For example, soil organic C ranged from 8.52 to 15.04 g/kg across land-use types in Beijing, China, but only from 10.19 to 11.68 g/kg along the urban–rural gradient (Mao et al., 2014).

Spatial heterogeneity is one of the main characteristics of urban ecosystems (Pickett et al., 2017). This characteristic is not only manifested in non-biological factors such as social economic activities and land use, but also reflected in biological components (Pouyat et al., 2007; Pickett et al., 2017). Within the city, there are great spatial variations in plants species, as well as in plant traits (Gregg et al., 2003; Guo et al., 2018). The result of this study shows that stomatal density, leaf width and leaf P concentration increased along the urban–rural gradient, while leaf K concentration decreased (Figure 3). These changes in plant traits would affect the ability of plants to use resources and strategies to respond to environmental changes, and subsequently ecosystem services beneficial to urban residents (Wright et al., 2004; de Bello et al., 2010).

The mechanism underlying changes in plant traits in cities are complex. In this study, the spatial variations in plant traits are affected by both human activities, such as land-use changes and urban-rural gradients caused by urbanization, and biological factors, such as leaf age (Figure 4). Human activities alter temperature, water and nutrients that are essential biophysical factors controlling plant development and growth, and then incur the changes in plant traits. Environmental changes alter plant physiological and ecological processes, then influencing plant production and allocation, and at last, plant traits would be altered. Therefore, measurements of the essential biophysical factors and related physiological parameters are necessary for understanding how human activities affect plant traits in urban settings.

Our study pinpoints that plant traits have changed with urban-rural gradient and across land use types and plant developmental stages. However, the genetic basis and fitness consequences of these trait changes are less explored. Manipulative experiments, such as common garden experiments or reciprocal translocation, would be useful to assess evolutionary and ecological consequences of plant trait changes.

The alteration of leaf functional traits by land-use type and urban–rural gradient indicates that plants are stressed in urban environments. It is therefore necessary to select tolerant plant species and conduct intensive management in urban planting. Pines we studied, even as a tolerant species, had smaller leaf length, width and area on roadsides than in parks and neighborhoods, implying that the production capacity of plants on roadsides may be reduced and that roadsides are poor to plant growth. The stresses faced by plants on roadsides might be mitigated by irrigation. The influence of urbanization on leaf traits is primarily dependent on the land-use type, management measures specific to the type of land use are necessary to guarantee plant health and growth.