*Article* **Does the Spatial Pattern of Plants and Green Space Affect Air Pollutant Concentrations? Evidence from 37 Garden Cities in China**

**Chengkang Wang 1,\*, Mengyue Guo <sup>1</sup> , Jun Jin <sup>2</sup> , Yifan Yang <sup>1</sup> , Yujie Ren <sup>3</sup> , Yang Wang <sup>4</sup> and Jiajie Cao 1,\***


**Abstract:** Relevant studies have demonstrated that urban green spaces composed of various types of plants are able to alleviate the morbidity and mortality of respiratory diseases, by reducing air pollution levels. In order to explore the relationship between the spatial pattern of urban green spaces and air pollutant concentrations, this study takes 37 garden cities with subtropical monsoon climate in China as the research object and selects the urban air quality monitoring data and land use type data in 2019 to analyze the relationship between the spatial pattern and the air pollutant concentration through the landscape metrics model and spatial regression model. Moreover, the threshold effect of the impact of green space on air pollutant concentrations is estimated, as well. The results showed that the spatial pattern of urban green space was significantly correlated with the concentrations of PM2.5 (PM with aerodynamic diameters of 2.5 mmor less), NO<sup>2</sup> (Nitrogen Dioxide), and SO<sup>2</sup> (Sulfur dioxide) pollutants in the air, while the concentrations of PM<sup>10</sup> (PM with aerodynamic diameters of 10 mmor less) pollutants were not significantly affected by the green space pattern. Among them, the patch shape index (LSI), patch density (PD) and patch proportion in landscape area (PLAND) of forest land can affect the concentration of PM2.5, NO<sup>2</sup> , and SO<sup>2</sup> , respectively. The PLAND, PD, and LSI of grassland and farmland can also have an additional impact on the concentration of SO<sup>2</sup> pollutants. The study also found that there was a significant threshold effect within the impact mechanism of urban green space landscape pattern indicators (LSI, PD, PLAND) on the concentrations of PM2.5, NO<sup>2</sup> , and SO<sup>2</sup> air pollutants. The results of this study not only clarified the impact mechanism of the spatial pattern of urban green space on air pollutant concentrations but also provided quantitative reference and scientific basis for the optimization and updating of urban green space to promote public health.

**Keywords:** public health; urban green spaces; landscape pattern; air pollution; quantitative analysis; threshold effect

### **1. Introduction**

The emergence and development of modern urban green spaces are closely related to the improvement of human living environments and public health [1]. As an important discipline for the development of human living environments, landscape architecture provides important support in the establishment of public safety systems in urban and rural spaces and promotes the health and wellness of residents [2]. With the accelerated global urbanization and rapid industrial development, including the frequent occurrence of public health incidents such as the new coronavirus pneumonia, the public health of residents is under unprecedented threat [3]. The severe challenges posed by the pandemic have prompted the entire society to deeply recognize the importance of public health and

**Citation:** Wang, C.; Guo, M.; Jin, J.; Yang, Y.; Ren, Y.; Wang, Y.; Cao, J. Does the Spatial Pattern of Plants and Green Space Affect Air Pollutant Concentrations? Evidence from 37 Garden Cities in China. *Plants* **2022**, *11*, 2847. https://doi.org/10.3390/ plants11212847

Academic Editors: Jie Gao, Weiwei Huang, Johan Gielis, Peijian Shi and Aneta Baczewska-D ˛abrowska

Received: 15 September 2022 Accepted: 24 October 2022 Published: 26 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

have triggered unprecedented attention to the urban built environment and the construction of healthy urban green spaces.

As a global environmental and public health problem, air pollution has severe adverse effects on human health. The pollutants can enter the body through the respiratory system and affect the lungs and heart, causing cardiovascular and respiratory diseases [4]. The World Health Organization (WHO) report demonstrated that respiratory diseases caused by air pollution in 2019 ranked fourth in the top 10 causes of death worldwide (https://www.who.int/zh/ (accessed on 15 July 2021)). Numerous studies have verified that short-term or long-term exposure to air pollutants, including PM2.5, PM10, NO2, and SO2, increased the risks of mortality and morbidity, thereby posing a severe threat to the public health of the residents [5]. As an important part of the urban ecosystem, in addition to providing residents with green spaces for recreation and entertainment, urban green spaces also improve air quality. The ability of different plant species, plant communities, and green spaces along roads to retain dust and mitigate air pollution has been widely demonstrated [6]. In addition, by optimizing the structure of urban green spaces, increasing the areas of green spaces and green coverage could effectively reduce the concentration of particulate matter and gas pollutants in the air, which effectively improves the urban environment, and ultimately play an important role in promoting public health [7].

Domestic and international research on the mitigation of air pollutants in green spaces at the micro-level has primarily focused on the ability of urban garden green spaces to reduce air pollutant concentrations, using individual plant microstructure processes such as sedimentation, retardation, adsorption, and absorption. For example, Latha and Highwood demonstrated that changes in structures, such as the roughness of plant leaf surfaces, affect the sedimentation pattern of dust particles [8]. Beckett et al. deduced that plants maintained a higher humidity in a certain range during transpiration, and dust sedimentation is more likely to occur when it increases in weight after absorbing moisture, while the ability of the leaves to adsorb dust increases with an increase in their humidity [9]. In addition, numerous studies have demonstrated the variability in the absorption of different gaseous pollutants by different landscape plants [10–12].

The meso-level was adopted to investigate the reduction effect of the size, shape, plant configuration, vertical structure, etc., of small- and medium-scale urban green spaces on different pollutants. In his study on the ecological mechanism of urban open space planning, Wang Shaozeng pointed out that the more practical the mix of green space levels, the better the filtration effect on the atmospheric particulate matter [13]. By analyzing the relationship between the three-dimensional green volume of green spaces along roads and PM2.5 concentration, Sheng found that high 3D green volume did not indicate low PM2.5 concentration, and such green spaces with the uniform vertical distribution of biomass and diverse vegetation were more effective in reducing PM2.5 concentration [14]. Fan analyzed the correlation between the daily PM<sup>10</sup> and PM2.5 concentrations and particulate matter concentrations of seven typical land cover types and different scales of land cover patterns and determined that pavement-type and low-to-medium canopy density vegetation exerted a more significant effect on PM<sup>10</sup> levels, while PM2.5 concentrations were more sensitive to the response of building-type and low-to-medium canopy density vegetation [15]. In several ways, the aforementioned studies verified that green spaces, as living media, are one of the most important vehicles for mitigating air pollution and play a significant role in the public health of residents and the urban environment.

Recently, studies on how to optimize the landscape pattern of urban green spaces to reduce air pollution concentrations at the macro-level have increasingly garnered considerable attention, and the current research results primarily focused on exploring the correlation between the concentrations of PM2.5, PM10, and other pollutants, including land use or the landscape patterns of land cover. Ye et al. [16] explored the relationship between PM2.5 growth and land use changes in China from 1998–2015 and inferred that PM2.5 concentrations were higher in the eastern plains and Taklamakan Desert in China, and higher PM2.5 concentrations existed on artificial land surfaces, croplands, and deserts, while forests, grasslands, and unused land usually contained lower PM2.5 concentrations. Simultaneously, the average annual increase in PM2.5 concentrations in a desert land and artificial land surfaces was higher than that of other land types. Yue et al. [17] studied the quantitative relationship between vegetation coverage and atmospheric particulate matter based on remote sensing inversion and deduced that vegetation coverages of ≤10% and >45% exhibited a significant effect on mitigating atmospheric particulate matter pollution. Lei et al. [18] explored the effect of green spaces on particulate matter pollution by studying the landscape patterns of urban green spaces at multiple spatial scales and determined that increasing the biodiversity of green spaces and increasing the number of large green spaces significantly reduced PM<sup>10</sup> concentrations at almost all scales. Zhao et al. [19] conducted a land use regression (LUR) analysis on the green spaces of lakes and wetlands, including the surrounding 500-m built environment in Wuhan City, and their results demonstrated that the lakes, wetlands, and nearby greenery exerted a positive and significant effect on PM2.5 concentrations within a buffer zone of 300 m or closer. Via a multivariate linear regression modeling of the daily average PM<sup>10</sup> concentration data with land use pattern information for the cities of Vienna and Dublin, McNabola et al. [20] determined that adding transboundary air pollution and traffic activity representations to the predictor variables significantly improved the accuracy of LUR-based methods. Zhang et al. [21] applied the LUR model to analyze the correlation between air pollution levels and childhood asthma hospitalization rates, by establishing a spatial distribution LUR model of the daily pollutant concentration data of PM<sup>10</sup> and SO<sup>2</sup> with associated influencing factors in Shenyang, China and deduced that the number of childhood asthma hospitalizations was highly correlated with PM<sup>10</sup> and SO<sup>2</sup> pollution levels. Lee and Koutrakis adopted satellite ozone monitoring instruments together with land use parameters to develop a mixed-effects model through which they estimated the daily NO<sup>2</sup> concentrations in New England, and explored the source areas of emissions (e.g., high population or traffic areas) in the study area and elucidated the seasonal characteristics of NO2, based on NO<sup>2</sup> spatial distribution patterns [22]. Lu et al. [23] analyzed the landscape pattern indices and PM2.5 data in China and demonstrated that differences exist in the significant impact indicators of different land use types on PM2.5 concentration and that the landscape pattern indices exhibit a significant effect on PM2.5 concentrations. De Jalón et al. [24] conducted a comparative analysis on the effect of the dry deposition of trees on air pollution reduction for different land use types in the Basque Country, and finally deduced that coniferous forests are the most effective in eliminating air pollution. These research results mainly elucidated the effect of urban land use changes on univariate air pollutants, and the pollutants were mainly PM2.5, PM10, and other fine particulate matter. A fewer number of studies have been conducted on NO2, SO2, and other gaseous pollutants, and these studies rarely involved any comprehensive air pollution research on the quantification and regulation strategies of the concentration of mixed pollutants of gases and particulate matter, including the spatial distribution of urban landscape patterns.

On this basis, this study takes 37 garden cities with subtropical monsoon climate in China as the research object, selects the urban air quality monitoring data and land use remote sensing data in 2019 to analyze the relationship between the landscape pattern index and the air pollutant concentration through the spatial regression model method. Also, a new threshold effect estimation method based on a polynomial model is designed to explore the threshold effect of green space landscape patterns on air quality [25,26].

### **2. Results**

In this paper, regression modeling tools in GeoDa software are used to conduct SEM regression analysis, and the results obtained are shown in Table 1.


**Table 1.** Conclusion of thresholds for landscape pattern indices.

\*\*\* *p* < 0.01, \*\* *p* < 0.05, \* *p* < 0.1. Correlations significant at the level of 0.1 are marked in bold.

### *2.1. Regression Analysis of Landscape Pattern Indices and Air Pollutants*

PM2.5 concentrations were significantly and negatively correlated with the LSI of forestlands; each unit increase in the landscape shape index of the forestlands was followed by a 0.68 unit decrease in the PM2.5 pollutant concentration. NO<sup>2</sup> concentration was significantly and negatively correlated with the PD of forestlands. With each unit of increase in the PD of forestlands, the concentration of the NO<sup>2</sup> pollutant was subsequently reduced by 258.409 units. However, Grasslands and farmlands had no significant correlation with PM2.5 and NO<sup>2</sup> pollutant concentrations. SO<sup>2</sup> concentrations were significantly and positively correlated with the PD of forestlands, PLAND, and LSI of grasslands, where the relationship between SO<sup>2</sup> pollutant concentrations and the PD of forestlands and LSI of grasslands was highly significant (*p* < 0.01). SO<sup>2</sup> concentration was significantly and negatively correlated with the PLAND of forestlands, PD of grasslands, and PLAND of farmlands, where the PD of grasslands and PLAND of farmlands were highly significantly correlated with the SO<sup>2</sup> pollutant concentration (*p* < 0.01). When the PD of forestlands, PLAND, and LSI of grasslands increased by one unit, the SO<sup>2</sup> concentration increased by 149.939, 0.752, and 0.429 units, respectively. When each unit of PLAND of forestlands, PD of grasslands, and PLAND of farmlands increased, the SO<sup>2</sup> concentration was reduced by 0.073, 214.564, and 0.172 units, respectively. In contrast, SO<sup>2</sup> concentrations were not significantly affected by the LSI of forestlands, and the PD and LSI of farmlands.

From the results of the spatial correlation analysis between the landscape pattern indices of green spaces and air pollutants in the 37 cities, it was determined that the landscape pattern of urban green spaces was significantly associated with the concentrations of PM2.5, NO2, and SO<sup>2</sup> pollutants in the air, while the concentrations of PM<sup>10</sup> pollutants were not significantly affected by the pattern of green spaces. Meanwhile, LSI, PD, and PLAND of forestlands influenced the concentration of the three types of air pollutants, respectively. In addition, the PLAND, PD, and LSI of grasslands and PLAND of farmlands exert an additional influence on the concentration of SO<sup>2</sup> pollutants.

### *2.2. Threshold Effect of the Impact of Landscape Pattern Indices on Air Quality*

To further validate the relationship between the landscape pattern indices and the concentrations of the four pollutants, this study conducted regression analysis and linear fitting on the four pollutant concentrations and their landscape pattern indices, with EXCEL revealing significant effects on them. After comparing the R<sup>2</sup> of the function's trend line, it was deduced that the overall effect of the polynomial fit was better than several other types of functions in the experiment; hence, the polynomial was chosen for linear fitting in this study.

### 2.2.1. Threshold Effect of the Impact of Landscape Pattern Indices on PM2.5

Table 1 indicated that the LSI of forest land has a significant negative effect on the concentration of PM2.5 (*p* < 0.1). Therefore, this study analyzes the influence trend of PM2.5

concentration change to determine the threshold range of forest land LSI. It is found that the concentration of PM2.5 will first decrease and then increase with the increase of forest land LSI. According to the fitting calculation of the polynomial function, the coordinate value of the polynomial vertex of the forest landscape shape index is [18.02, 37.94]. Therefore, when the landscape shape index of the forest is 18.02, the minimum concentration of PM2.5 reaches 37.94 µg/m<sup>3</sup> . Therefore, this study found that the LSI value of 18.02 is the threshold of forest land (Figure 1). the concentration of PM2.5 will first decrease and then increase with the increase of forest land LSI. According to the fitting calculation of the polynomial function, the coordinate value of the polynomial vertex of the forest landscape shape index is [18.02, 37.94]. Therefore, when the landscape shape index of the forest is 18.02, the minimum concentration of PM2.5 reaches 37.94 ug/m³. Therefore, this study found that the LSI value of 18.02 is the threshold of forest land (Figure 1).

2.2.1. Threshold Effect of the Impact of Landscape Pattern Indices on PM2.5

To further validate the relationship between the landscape pattern indices and the concentrations of the four pollutants, this study conducted regression analysis and linear fitting on the four pollutant concentrations and their landscape pattern indices, with EX-CEL revealing significant effects on them. After comparing the R2 of the function's trend line, it was deduced that the overall effect of the polynomial fit was better than several other types of functions in the experiment; hence, the polynomial was chosen for linear

Table 1 indicated that the LSI of forest land has a significant negative effect on the concentration of PM2.5 (*p* < 0.1). Therefore, this study analyzes the influence trend of PM2.5 concentration change to determine the threshold range of forest land LSI. It is found that

*Plants* **2022**, *11*, x FOR PEER REVIEW 5 of 21

fitting in this study.

**Figure 1.** Threshold of significant landscape pattern indices for PM2.5 pollutant. **Figure 1.** Threshold of significant landscape pattern indices for PM2.5 pollutant.

2.2.2. Threshold Effect of the Impact of Landscape Pattern Indices on NO<sup>2</sup>

2.2.2. Threshold Effect of the Impact of Landscape Pattern Indices on NO2 From Table 1, it was clear that the PD of forestland exerted a significant effect on NO2 concentration (*p* < 0.1); therefore, the study determined the range of both thresholds by analyzing the trend of the effect of NO2 concentration changes. From the equation curve in Figure 2, it can be observed that the NO2 concentration exhibited a decreasing trend, and then started increasing with the increase in the PD of the forestland. From the function calculation, the value of the polynomial vertex coordinates of the PD of forestland was [0.072, 29.161]; hence, when the forestland patch density was 0.072, the NO2 concentration From Table 1, it was clear that the PD of forestland exerted a significant effect on NO<sup>2</sup> concentration (*p* < 0.1); therefore, the study determined the range of both thresholds by analyzing the trend of the effect of NO<sup>2</sup> concentration changes. From the equation curve in Figure 2, it can be observed that the NO<sup>2</sup> concentration exhibited a decreasing trend, and then started increasing with the increase in the PD of the forestland. From the function calculation, the value of the polynomial vertex coordinates of the PD of forestland was [0.072, 29.161]; hence, when the forestland patch density was 0.072, the NO<sup>2</sup> concentration reached the minimum of 29.161 µg/m<sup>3</sup> . Therefore, 0.072 was deduced to be the threshold value for the PD of forestland in this study (Figure 2).

reached the minimum of 29.161 ug/m³. Therefore, 0.072 was deduced to be the threshold

value for the PD of forestland in this study (Figure 2).

**Figure 2.** Threshold of significant landscape pattern indices for NO2 pollutant. **Figure 2.** Threshold of significant landscape pattern indices for NO<sup>2</sup> pollutant.

2.2.3. Threshold Effect of the Impact of Landscape Pattern Indices on SO<sup>2</sup>

2.2.3. Threshold Effect of the Impact of Landscape Pattern Indices on SO2 PD and PLAND of the forest, PD, and LSI of grassland, and PLAND of farmland will affect the concentration of SO2 pollutants, significantly. Through linear analysis of SO2 concentration change influence trend, the threshold results of 6 indicators are determined PD and PLAND of the forest, PD, and LSI of grassland, and PLAND of farmland will affect the concentration of SO<sup>2</sup> pollutants, significantly. Through linear analysis of SO<sup>2</sup> concentration change influence trend, the threshold results of 6 indicators are determined as follows.

as follows. For forest land, the concentration of SO2 will increase first and then decrease in the forest. Through polynomial fitting, it is found that when the PLAND of forest value is 50%, the SO2 concentration reaches the maximum value of 8.7 ug/m³. When the PD value of forest land is 0.038, the SO2 concentration reaches the maximum value of 8.72 ug/m³. In For forest land, the concentration of SO<sup>2</sup> will increase first and then decrease in the forest. Through polynomial fitting, it is found that when the PLAND of forest value is 50%, the SO<sup>2</sup> concentration reaches the maximum value of 8.7 µg/m<sup>3</sup> . When the PD value of forest land is 0.038, the SO<sup>2</sup> concentration reaches the maximum value of 8.72 µg/m<sup>3</sup> . In consideration of air quality and public health of residents, the optimal range of PD and PLAND of forest land is greater than or less than 50% and 0.038.

consideration of air quality and public health of residents, the optimal range of PD and PLAND of forest land is greater than or less than 50% and 0.038. For grassland, the concentration of SO2 will increase first and then decrease with the increase of PD and PLAND of grassland. When the PLAND value of grassland is 3.33%, the SO2 concentration reaches the maximum value of 8.72 ug/m³. When the PD value of grassland is 0.121, the SO2 concentration reaches the maximum value of 9.545 ug/m³. Also, the SO2 concentration will decrease first and then increase with the increase of LSI of grass-For grassland, the concentration of SO<sup>2</sup> will increase first and then decrease with the increase of PD and PLAND of grassland. When the PLAND value of grassland is 3.33%, the SO<sup>2</sup> concentration reaches the maximum value of 8.72 µg/m<sup>3</sup> . When the PD value of grassland is 0.121, the SO<sup>2</sup> concentration reaches the maximum value of 9.545 µg/m<sup>3</sup> . Also, the SO<sup>2</sup> concentration will decrease first and then increase with the increase of LSI of grassland, when the grassland LSI value is 14.13, the SO<sup>2</sup> concentration reaches the minimum value of 7.84 µg/m<sup>3</sup> .

land, when the grassland LSI value is 14.13, the SO2 concentration reaches the minimum value of 7.84 ug/m³. For farmland, the concentration of SO2 will increase first and then decrease with the increase of PLAND of farmland. When the PLAND value of farmland is 32.56, the PLAND For farmland, the concentration of SO<sup>2</sup> will increase first and then decrease with the increase of PLAND of farmland. When the PLAND value of farmland is 32.56, the PLAND value of 8.78 µg/m<sup>3</sup> . Therefore, this study found that 32.56 is the threshold of the PLAND value of farmland (Figure 3).

value of 8.78 ug/m³. Therefore, this study found that 32.56 is the threshold of the PLAND

value of farmland(Figure 3).

**Figure 3.** Threshold of significant landscape pattern indices for SO2 pollutant.(**a**−**f** represents the threshold effect of SO2 and forest land PLAND, forest land PD, grassland PLAND, grassland PD, **Figure 3.** Threshold of significant landscape pattern indices for SO<sup>2</sup> pollutant ((**a**–**f**) represents the threshold effect of SO<sup>2</sup> and forest land PLAND, forest land PD, grassland PLAND, grassland PD, grassland LSI, and farmland PLAND landscape pattern index).

### grassland LSI, and farmland PLAND landscape pattern index). **3. Discussions**

#### **3. Discussions**  *3.1. Impact Mechanism of Landscape Pattern Index on Air Pollutant Concentration*

### *3.1. Impact Mechanism of Landscape Pattern Index on Air Pollutant Concentration*  3.1.1. Impact Mechanism of Landscape Pattern Index on PM2.5

3.1.1. Impact Mechanism of Landscape Pattern Index on PM2.5 The results of spatial regression analysis showed that the LSI of forest land had a significant negative effect on the concentration of PM2.5, which was consistent with the results obtained by previous studies and indicated that increasing the contact area between the edges of the green space patches and surrounding urban areas at large spatial scales significantly reduces PM2.5 concentrations [18]. PM2.5 primarily originates from industrial emissions, traffic emissions, and the burning of biomass, and is mainly present on roads, factories, and surrounding farmland in built-up areas [27]. When the LSI of forestland was elevated, the complexity of its patch edge also increased, and the contact area between the forestland patch edge and urban construction land was subsequently elevated. As the contact area between vegetation and PM2.5 particulate matter increased, the effect of dust reduction via vegetation leaf surface villi retardation, stem adsorption, and The results of spatial regression analysis showed that the LSI of forest land had a significant negative effect on the concentration of PM2.5, which was consistent with the results obtained by previous studies and indicated that increasing the contact area between the edges of the green space patches and surrounding urban areas at large spatial scales significantly reduces PM2.5 concentrations [18]. PM2.5 primarily originates from industrial emissions, traffic emissions, and the burning of biomass, and is mainly present on roads, factories, and surrounding farmland in built-up areas [27]. When the LSI of forestland was elevated, the complexity of its patch edge also increased, and the contact area between the forestland patch edge and urban construction land was subsequently elevated. As the contact area between vegetation and PM2.5 particulate matter increased, the effect of dust reduction via vegetation leaf surface villi retardation, stem adsorption, and stomatal absorption by plants was enhanced, while the concentration of PM2.5 pollutants reduced [28,29].

#### stomatal absorption by plants was enhanced, while the concentration of PM2.5 pollutants 3.1.2. Impact Mechanism of Landscape Pattern Index on NO<sup>2</sup>

reduced [28,29]. 3.1.2. Impact Mechanism of Landscape Pattern Index on NO2 The results of spatial regression analysis showed that the PD of forest land had a The results of spatial regression analysis showed that the PD of forest land had a significant negative effect on the concentration of NO2, which indicated that increasing the contact area of forestland with NO<sup>2</sup> pollutants was beneficial to the reduction of NO<sup>2</sup> pollutant concentration.

significant negative effect on the concentration of NO2, which indicated that increasing

The main sources of NO<sup>2</sup> in the atmosphere are industrial and motor vehicle exhaust emissions [30]. Therefore, NO<sup>2</sup> air pollution mainly exists in urban industrial areas with large traffic emissions. Studies have demonstrated that plants can absorb NO<sup>2</sup> through their leaves [31]. The forestland patches outside the urban constructed land have a large area; however, NO<sup>2</sup> was mainly generated in the built-up area and blocked by high-density buildings; hence, the contact area between urban forestland and NO<sup>2</sup> pollutant gas was small, and increasing the forestland area exhibited no significant effect on the weakening effect of the overall urban NO<sup>2</sup> pollutant concentration. Increasing the density of forestland patches in the built-up area enhanced the number of public green spaces in the city, and measures were taken to insert greenery in the edges, especially the greening of industrial areas, and street greenery was appropriately sufficient in moving plant leaves to come into contact with more polluting gases while playing a more significant role in the improvement of urban NO<sup>2</sup> pollutant gases.

### 3.1.3. Impact Mechanism of Landscape Pattern Index on SO<sup>2</sup>

The results of spatial regression analysis showed that the PLAND of forest land, PD of grassland, and PLAND of farmland had a significant negative effect on the concentration of SO2. And the PD of forest land, PLAND of grassland, and LSI of grassland had a significant positive effect on the concentration of SO2.

The production of SO<sup>2</sup> in urban air pollution mainly originates from coal combustion and industrial emissions, and studies have demonstrated that SO<sup>2</sup> concentrations were lower in areas with rich vegetation under different spaces in cities [32]; hence, the increase in the area share of forestland patches effectively enhanced the efficiency of SO<sup>2</sup> absorption by urban plants and reduced the SO<sup>2</sup> air content in cities. Simultaneously, as a gas, SO<sup>2</sup> was more susceptible to wind speed, temperature, and humidity; hence, the single vertical structure of grass greenery was more conducive to the circulation of urban winds and the transportation of SO<sup>2</sup> pollutants from the inner city to the outer city, to achieve a lower average SO<sup>2</sup> concentration [33]. The diffusion of SO<sup>2</sup> is mainly influenced by wind speed and direction, and in areas with high building density due to the blockage of tall buildings resulting in low internal wind speed, SO<sup>2</sup> concentration was not easily diffused. Therefore, with the increase in the proportion of farmland area, the degree of SO<sup>2</sup> diffusion was accelerated and the concentration was reduced [34]. Meanwhile, cities with relatively large areas of farmland generally have industries positioned as agricultural cities, such as Henan Province, which is a largely agricultural province in China, including Xinyang and Nanyang, and these areas have low average SO<sup>2</sup> concentrations.

As the density of forestland patches increased, SO<sup>2</sup> concentration increased, and related studies have demonstrated that the larger the average area of green space patches in urban landscapes, the lower the fragmentation index, and the greater the role of green spaces in air pollution purification. In addition, when the degree of fragmentation of forestland was more severe, it could not function as an effective urban green heart, and the fragility of the landscape structure reduced the absorption capacity of SO<sup>2</sup> [35]. As the ratio of the area occupied by grass increased, SO<sup>2</sup> concentration increased. It was speculated that the reason for this result might be that the decrease in the number of plant leaves (the main organ of SO<sup>2</sup> absorption by plants) in grass patches resulted in the subsequent decrease in the efficiency of SO<sup>2</sup> absorption by green spaces [36]. Therefore, it was inferred that in addition to improving the air purification capacity, increasing the density of grass patches with small areas and low edge complexity also reduces economic costs.

### *3.2. Threshold Mechanism of Landscape Pattern Index on Air Pollutant Concentration* 3.2.1. Threshold Effect of PM2.5

The LSI of forest land is related to the edge complexity and patch density of green space patches. The increase in LSI value can enhance the energy flow and exchange between green patches and surrounding patches and create more interaction opportunities for source and sink landscapes, so as to absorb and settle more pollution particles and reduce the concentration of PM2.5 [37]. However, as the LSI of forest land continues to increase, the patch density then increases the degree of green space fragmentation increases, and the abatement between forest land edges and air pollutants cannot offset the increasing amount of air pollution due to urban green space fragmentation, and PM2.5 concentrations then continue to rise to show a trend of first decreasing and then increasing [38].

### 3.2.2. Threshold Effect of NO<sup>2</sup>

As the PD of forest land keeps increasing, NO<sup>2</sup> concentration shows a trend of first decreasing and then increasing. When the density of forested land patches increases continuously from 0, the intra-urban green space gradually transitions from single-core large green space to multi-core green space. With the connection of streets, rivers, and other green channels, the multi-core urban green space plays a greater ecosystem service function, more contact area of air pollutants, and higher efficiency of material exchange and energy circulation, which is conducive to the reduction of NO<sup>2</sup> pollutant concentration [39]. However, with the increasing density of forest land patches, the urban green space system transitions from multi-core green space to green space patch fragmentation, and the concentration of NO<sup>2</sup> pollutants then increases [40].

### 3.2.3. Threshold Effect of SO<sup>2</sup>

Except for the LSI of grassland, with the increasing values of PLAND and PD of forest land, PLAND and PD of grassland, and PLAND of agricultural land, the trend between the landscape pattern index and SO<sup>2</sup> concentrations of different green space types was first increasing and then decreasing.

Green space is a sink landscape for mitigating urban SO<sup>2</sup> pollution, which can adsorb and absorb and deter SO<sup>2</sup> through plant leaves and branches, while urban construction land is a source landscape for SO<sup>2</sup> pollutants [41]. The decrease in the area share of forest land patches is often associated with cities with high levels of development. According to relevant studies showing more, developed cities with relatively well-developed environmental measures and more rational urban master plans, as well as industrial structures with the upgrading and transformation of heavy industries to light industries and high-tech industries, tend to have lower pollution levels. The opposite is often true in medium-sized developing cities. Such cities are still at the stage of economic development, with more types of industries and insufficient attention to the environment, resulting in the reduction of SO<sup>2</sup> air pollution by vegetation in such cities being much lower than the SO<sup>2</sup> emissions, and therefore the pollution concentration is higher. With the increasing area of forest land emissions and reductions gradually reaching a balance or even vegetation reduction exceeding the local pollution emissions air pollution levels are reduced again [42,43].

### *3.3. Implications for Urban Planning and Management Policies*

In the context of optimizing and renewing the landscape patterns of urban green spaces, several studies have demonstrated that regulating land use patterns by carefully planning the morphology and layout of green space networks effectively improves air quality and enhances the public health of residents. Based on the aforementioned findings, we made the following recommendations for the improvement of air pollution in cities with subtropical monsoon climates.

Forestland, grassland, and farmland can inhibit the concentration of PM2.5, NO2, and SO2. The grassland area in the urban built-up area should be properly controlled, the forest coverage should be gradually increased, the restoration and reconstruction of damaged forest land should be accelerated, the integrity and stability of the ecosystem should be improved, the urban green space landscape pattern should be reasonably controlled through scientific basis, and the production and living ecological space layout should be coordinated.

For cities with NO<sup>2</sup> and SO<sup>2</sup> as the main pollutants, the pollutants can be reduced systematically by reasonably arranging the land use pattern of forest land, grassland, farmland, and other green areas, and cooperating with each other. When the PD value of forest land is about 0.072, the overall value of urban NO<sup>2</sup> pollutant concentration reaches the optimum. When the land and PD values of forest land are away from 50% and 0.038, the land and PD values of grassland are away from 3.33% and 0.121, and the LSI of grassland reaches 14.13, the urban SO<sup>2</sup> pollutant concentration reduction effect is the best. Therefore, reasonable urban green space planning should be carried out to optimize the density of green space patches and evenly distribute the types of green space patches, to achieve a close relationship between various land uses and better alleviate air pollution. Through the way of inserting green in the gap, the coverage of urban green space can be improved economically and efficiently, and the contact area between green space and air pollutants can be increased to better improve the air quality level.

PM2.5 mainly comes from industrial waste gas emission, traffic emission, and biomass combustion, and mainly exists in roads and factories in built-up areas. China's subtropical monsoon climate zone is in the stage of rapid urban development, and it is difficult to reduce traffic emissions and industrial emissions. However, PM2.5 pollution can be improved by optimizing the landscape pattern of urban green space. When the LSI of forest land reaches 18.02, the overall value of PM2.5 pollution concentration is the best. It is suggested to reasonably arrange the urban green space area. For the streets and factories with large traffic and industrial emissions, we will focus on improving the greening level of the streets, appropriately increasing the contact area between the street green space and PM2.5, improving the vertical structure of the urban green belt, block, and absorb PM2.5 pollutants to the greatest extent, and reduce the transmission route. In addition, the treatment efficiency of polluted gas in the factory shall be strictly controlled to reduce air pollution.

The comprehensive management of unused land and inefficient land should be promoted to form a reasonable and efficient urban green landscape pattern. At the same time, it is suggested to promote the development of residents' lifestyle and consumption concepts towards green, healthy, and low-carbon. The scientific achievements in the prevention and control of PM2.5, NO2, and SO<sup>2</sup> pollution should be actively disseminated to the residents so that the relevant departments can implement the relevant pollution control measures more smoothly, and the residents can support the air pollution control.

In addition, the results of the spatial regression model in this study also show that air pollution has obvious spatial effects, so it is not feasible to carry out internal air pollution control for a single city. From the perspective of urban agglomerations, we should coordinate and plan the air pollution control policies among cities, jointly improve the regional air quality level and reduce the threat of air pollution to residents' public health. At the same time, our results show that the careful planning of urban green space landscape patterns has brought some positive benefits to air pollution, but compared with traffic emissions, industrial emissions, and biomass combustion, the role of pollutant gas emissions is still limited. Therefore, improving the rationality of urban green space landscape patterns is on the one hand. On the other hand, we should also attach great importance to energy efficiency and traffic management, further promote the reform of industrial structure, and eliminate backward industries with high pollution and high consumption. Finally, it aims to increase the proportion of high-tech industries, promote China's transformation from incremental expansion to stock revitalization, improve development quality and resource utilization efficiency, and reduce air pollutant emissions from the source.

### *3.4. Research Innovations and Limitations*

As early as the end of the last century Wickham et al. [44] proposed an integrated assessment of the environmental condition of a large region, by combining data on land cover, population, roads, rivers, air pollution, and topography. In 2009, Rafiee et al. [45] applied a combination of remote sensing image classification, landscape indicator assessment, and vegetation indices, to explore changes in urban landscape patterns and provide an assessment of changes and trends in urban living environments. Recently, there has been

an increasing interest in applying remote sensing images to investigate the role of changes in urban land use patterns on air pollution concentrations. In 2015, Wu et al. [46] adopted landscape indicators such as PLAND, PD, ED, SHEI, and CONTAG to explore the effect of urban landscape patterns on PM2.5 pollution and determined that vegetation and water bodies were significant landscape components that reduced PM2.5 concentrations. In 2016, Xu et al. [47] explored the quantitative relationship between land use and air quality (SO2, NO2, and PM10) through binary correlation analysis, and the results indicated that land use had a significant effect on air quality. For each standard deviation increase in construction land, NO<sup>2</sup> concentrations increased by 2%. For one standard deviation increase in water bodies, SO<sup>2</sup> or PM<sup>10</sup> concentrations decreased by 3–6%. Although this study quantified the role of land use type on air pollution levels, it did not explore the quantitative relationship between land use effects on PM2.5. In 2020, Li et al. [48] studied the non-linear effects of land use distribution on PM2.5, using the boosted regression tree method to capture the effects of land use scale on PM2.5 in different seasons. It was inferred that when the grassland and forestland areas were below 8% and 20%, respectively, the air quality improved significantly with the increase in grassland and forestland areas. When the distribution of construction land was greater than approximately 10%, PM2.5 pollution increased significantly with the increase in the construction land area. This study was the first to identify the threshold for the effect of land use type on PM2.5 concentrations. However, the role of the effect on the three air pollutants such as PM10, NO2, and SO<sup>2</sup> is yet to be addressed. Most of the aforementioned studies adopted linear regression models to explore the role of land use types on air pollutants, which could not comprehensively incorporate the spatial effects between cities into the study, and also failed to demonstrate how the effects of landscape patterns on multiple air pollutants varied with changes in landscape pattern indices. On this premise, this study employed a spatial regression model approach to regress air pollutant concentrations and landscape pattern indices in 37 subtropical monsoon climate garden cities of China in 2019, to explore the core landscape pattern indices that exhibited significant effects on air pollution concentrations and the optimal core indicator thresholds for mitigating urban air pollution, to quantify the spatial effects of the landscape pattern indices on air pollutants, using urban spatial effects as a starting point, and then a general framework that differed from existing studies was proposed. In addition to focusing on quantifying the negative or positive effects of landscape pattern indices on air pollutant concentrations, to a certain extent, this study also reflected the appropriate threshold values of landscape pattern indices for reducing air pollution concentrations, which provides quantitative reference and technical support for urban planning in a more targeted manner.

However, in general, this study had some limitations. Firstly, the spatial resolution of satellite remote sensing images in this study is low, and there may be a certain misclassification of patch types, which may cause deviation to the impact indicators of air pollution. At the same time, because the research object is a large city with a subtropical monsoon climate, it is not universal for other climate belt cities. In addition, the image data used in our study provide a reference for the overall land use pattern of the whole city. The conclusion is that it provides a reference for the overall planning of urban green space landscape patterns from a macro perspective. The role of small-scale urban green space in reducing pollutants needs to be further discussed. At the same time, in order to ensure the accuracy of the urban green space landscape pattern indicators, the image data selected in this study are all summer image data, but the air pollution data is the annual average value of pollutants. Due to the seasonal changes, the green space coverage of a few cities in winter is reduced, so the impact analysis of this study on air pollutants is still insufficient. Therefore, based on the above shortcomings, in future research, we will expand the sample number of research objects, improve the accuracy of remote sensing data, and analyze the urban landscape pattern and air pollutant concentration in different seasons and climate zones respectively, so as to reduce variables, improve universality and research accuracy, and put forward more targeted urban green space landscape pattern planning strategies.

### **4. Materials and Methods** In this study, firstly, we take 37 garden cities with subtropical monsoon climate as the research unit and take the annual average concentration of PM2.5, PM10, NO2, and SO2

**4. Materials and Methods** 

planning strategies.

In this study, firstly, we take 37 garden cities with subtropical monsoon climate as the research unit and take the annual average concentration of PM2.5, PM10, NO2, and SO<sup>2</sup> and the landscape pattern index of these 37 cities in 2019 as the dependent variable and independent variable, respectively, to carry out the spatial regression model. Based on the output of the model, the impact mechanism between landscape pattern index and PM2.5, PM10, NO2, and SO<sup>2</sup> pollution was explored. Finally, a new threshold effect estimation method based on a polynomial model is designed to explore the threshold effect of green space landscape patterns on air quality (Figure 4). and the landscape pattern index of these 37 cities in 2019 as the dependent variable and independent variable, respectively, to carry out the spatial regression model. Based on the output of the model, the impact mechanism between landscape pattern index and PM2.5, PM10, NO2, and SO2 pollution was explored. Finally, a new threshold effect estimation method based on a polynomial model is designed to explore the threshold effect of green space landscape patterns on air quality (Figure 4).

The conclusion is that it provides a reference for the overall planning of urban green space landscape patterns from a macro perspective. The role of small-scale urban green space in reducing pollutants needs to be further discussed. At the same time, in order to ensure the accuracy of the urban green space landscape pattern indicators, the image data selected in this study are all summer image data, but the air pollution data is the annual average value of pollutants. Due to the seasonal changes, the green space coverage of a few cities in winter is reduced, so the impact analysis of this study on air pollutants is still insufficient. Therefore, based on the above shortcomings, in future research, we will expand the sample number of research objects, improve the accuracy of remote sensing data, and analyze the urban landscape pattern and air pollutant concentration in different seasons and climate zones respectively, so as to reduce variables, improve universality and research accuracy, and put forward more targeted urban green space landscape pattern

*Plants* **2022**, *11*, x FOR PEER REVIEW 12 of 21

**Figure 4.** Flow chart Schemes follow the same formatting.

#### **Figure 4.** Flow chart Schemes follow the same formatting. *4.1. Study Region*

*4.1. Study Region*  The subtropical monsoon climate zone is one of the four major climate zones in China, covering approximately a quarter of China's land area. It supports almost 600 million people in China, is the most densely populated region in China, and mainly includes The subtropical monsoon climate zone is one of the four major climate zones in China, covering approximately a quarter of China's land area. It supports almost 600 million people in China, is the most densely populated region in China, and mainly includes the Yangtze River Delta region, Pearl River Delta region, and middle reaches of the Yangtze River, including other areas of China's economic center, where the level of economic development and urban landscaping are at a relatively high level in the country. With the advancement of urbanization, the original urban ecological network structure has been destroyed, owing to the drastic changes in urban land use patterns caused by the population concentration and rapid industrial development; in addition, the subtropical monsoon climate region has also become the region most threatened by air pollution in China, which exposes the public health of residents to tremendous pressure [49,50]. Considering the substantial influence of natural climate on urban air quality, to improve the accuracy of this study, 37 "national garden cities" with a population size of more than 4 million people and complete air quality data for 2019 were selected from the subtropical monsoon climate zone of China as the research objects of this study. According to the "National Garden City Standards" adopted by the Ministry of Housing and Urban-Rural Development of the People's Republic of China, garden cities are cities with balanced distribution, practical

structure, perfect functions, beautiful landscapes, fresh and comfortable human living, and ecological environments, safe, and pleasant. In addition, they were adopted as examples to investigate the effect of landscape patterns of urban green spaces on air pollution levels (see Figure 5, Appendix A). living, and ecological environments, safe, and pleasant. In addition, they were adopted as examples to investigate the effect of landscape patterns of urban green spaces on air pollution levels (see Figure 5, Appendix A).

the Yangtze River Delta region, Pearl River Delta region, and middle reaches of the Yangtze River, including other areas of China's economic center, where the level of economic development and urban landscaping are at a relatively high level in the country. With the advancement of urbanization, the original urban ecological network structure has been destroyed, owing to the drastic changes in urban land use patterns caused by the population concentration and rapid industrial development; in addition, the subtropical monsoon climate region has also become the region most threatened by air pollution in China, which exposes the public health of residents to tremendous pressure [49,50]. Considering the substantial influence of natural climate on urban air quality, to improve the accuracy of this study, 37 "national garden cities" with a population size of more than 4 million people and complete air quality data for 2019 were selected from the subtropical monsoon climate zone of China as the research objects of this study. According to the "National Garden City Standards" adopted by the Ministry of Housing and Urban-Rural Development of the People's Republic of China, garden cities are cities with balanced distribution, practical structure, perfect functions, beautiful landscapes, fresh and comfortable human

*Plants* **2022**, *11*, x FOR PEER REVIEW 13 of 21

**Figure 5.** Distribution map of the selected 37 cities. **Figure 5.** Distribution map of the selected 37 cities.

### *4.2. Data Sources*

*4.2. Data Sources*  Air quality index data were obtained from the "2019 National Air Quality Monthly Report" published by the Ministry of Ecology and Environment of the People's Republic of China (http://www.mee.gov.cn/ (accessed on 13 September 2020)). In addition, the specific distributions of the air pollution index data (PM2.5, PM10, SO2, and NO2) for the 37 cities are presented in Figure 6, while the population data were obtained from the "2020 City Yearbook" for each city provided by the National Bureau of Statistics of China (http://www.stats.gov.cn/ (accessed on 1 September 2020)), and the number of permanent residents in 2019 was selected as the population index (Figure 6). The 37 urban land use classification maps in this study were obtained from the MCD12Q1.006 land use type product of the 2019 MODIS/Terra, provided by the official USGS online platform with a resolution of 500 m (https://www.usgs.gov/ (accessed on 20 September 2020)). The urban Air quality index data were obtained from the "2019 National Air Quality Monthly Report" published by the Ministry of Ecology and Environment of the People's Republic of China (http://www.mee.gov.cn/ (accessed on 13 September 2020)). In addition, the specific distributions of the air pollution index data (PM2.5, PM10, SO2, and NO2) for the 37 cities are presented in Figure 6, while the population data were obtained from the "2020 City Yearbook" for each city provided by the National Bureau of Statistics of China (http://www.stats.gov.cn/ (accessed on 1 September 2020)), and the number of permanent residents in 2019 was selected as the population index (Figure 6). The 37 urban land use classification maps in this study were obtained from the MCD12Q1.006 land use type product of the 2019 MODIS/Terra, provided by the official USGS online platform with a resolution of 500 m (https://www.usgs.gov/ (accessed on 20 September 2020)). The urban land use map was converted to a tiff format by ArcGIS10.3 and then imported into FRAGSTATS for calculations, to obtain the data for patch-level landscape pattern indices (Appendix B).

pendix B).

**Figure 6.** Distribution map of PM2.5, PM10, NO2, and SO2 concentration. (**a**) Distribution map of PM2.5 concentration in the study city. (**b**) Distribution map of PM10 concentration in the study city. (**c**) Distribution map of NO2 concentration in the study city. (**d**) Distribution map of SO2 concentration in the study city). **Figure 6.** Distribution map of PM2.5, PM10, NO2, and SO<sup>2</sup> concentration. (**a**) Distribution map of PM2.5 concentration in the study city. (**b**) Distribution map of PM<sup>10</sup> concentration in the study city. (**c**) Distribution map of NO<sup>2</sup> concentration in the study city. (**d**) Distribution map of SO<sup>2</sup> concentration in the study city.

land use map was converted to a tiff format by ArcGIS10.3 and then imported into FRAG-STATS for calculations, to obtain the data for patch-level landscape pattern indices (Ap-

### *4.3. Landscape Pattern Indices*

*4.3. Landscape Pattern Indices*  To some extent, the landscape pattern indices capture the intrinsic spatial structure of the environment and enhance the interpretation of spatial patterns and characteristics of the landscape and are now widely used to measure landscape patterns [51]. Based on the results obtained from previous related studies and validation experiments, we selected three landscape indicators to measure the urban landscape pattern of the study region, these indexes are patch proportion in landscape area (PLAND), patch density (PD), and landscape shape index (LSI). The reason why these three indicators are selected is that they are not only widely used to describe the fragmentation, irregularity, and complexity of urban landscape patterns, but also have a more direct macro control effect on the adjustment and optimization of urban green space landscape pattern in the later stage, which is conducive to clarifying and standardizing the follow-up policymaking. Patch proportion of landscape area (PLAND) is the proportion of different landscape patch types in the overall land area percentage measurement, which can help us judge the proportion of this patch type in the overall spatial pattern The patch density (PD) represents the number of patches of a certain type within 100 hectares and can reflect the spatial pattern of the landscape. Its value has a positive correlation with the fragmentation of patch types. The greater the density index, the higher the fragmentation of the patch type. The landscape To some extent, the landscape pattern indices capture the intrinsic spatial structure of the environment and enhance the interpretation of spatial patterns and characteristics of the landscape and are now widely used to measure landscape patterns [51]. Based on the results obtained from previous related studies and validation experiments, we selected three landscape indicators to measure the urban landscape pattern of the study region, these indexes are patch proportion in landscape area (PLAND), patch density (PD), and landscape shape index (LSI). The reason why these three indicators are selected is that they are not only widely used to describe the fragmentation, irregularity, and complexity of urban landscape patterns, but also have a more direct macro control effect on the adjustment and optimization of urban green space landscape pattern in the later stage, which is conducive to clarifying and standardizing the follow-up policymaking. Patch proportion of landscape area (PLAND) is the proportion of different landscape patch types in the overall land area percentage measurement, which can help us judge the proportion of this patch type in the overall spatial pattern The patch density (PD) represents the number of patches of a certain type within 100 hectares and can reflect the spatial pattern of the landscape. Its value has a positive correlation with the fragmentation of patch types. The greater the density index, the higher the fragmentation of the patch type. The landscape shape index (LSI) is a robust index used to describe the complexity of urban morphology through the ratio of urban patch perimeter. The larger the LSI value, the more irregular the patch shape, the higher the landscape complexity, and the lower the stability. The above indicators are calculated using FRAGSTATS 4.2 software (Table 2).

shape index (LSI) is a robust index used to describe the complexity of urban morphology through the ratio of urban patch perimeter. The larger the LSI value, the more irregular


**Table 2.** Landscape pattern indices.

### *4.4. Spatial Regression Modeling Methods*

In this study, to investigate the quantitative influence of the landscape patterns of urban green spaces on air pollutants, the PLAND, PD, and LSI of three land use types (forestland, grassland, and farmland) were selected as the independent variables of air pollutant concentration that influence factors for spatial regression analysis.

First, in the selection of the model, after the regression analysis. The calculation of Moran's I index for air pollutants using the GeoDa software revealed that the Moran's I index was 0.57, and the model exhibited a significant spatial correlation. This indicated that air pollutants exhibit spatial interactions, and their effects could spread through adjacent regions, which confirmed that local and regional landscape pattern changes directly or indirectly affected their surrounding air pollutant concentrations. Therefore, the spatial regression model can better reveal the relationship between air pollutant concentration and landscape green space pattern.

Secondly, the regression results of the spatial autoregression model (SAR) and spatial error model (SEM) were compared in Table 3. The evaluation indicators for the goodness of fit of the spatial regression model included the coefficient of determination R2, the natural logarithm of the likelihood function (log-likelihood, logL), Akaike information criterion (AIC), and Schwarz Criterion (SC). Among them, the value range of R2 is (0,1), and the closer R2 is to 1, the better the regression fit of the model. In addition, the higher the logL value, and the smaller the AIC and SC values, the better the regression effect of the spatial regression model. For the selection of the optimal spatial regression model, the R2 and logL values of the SEM model were higher than those of the SAR model in the regression model information of the four pollutants, while the AIC and SC values of the SEM model were smaller for the remaining three pollutants, except for PM10 of the SEM model, which was 0.04 greater than the SAR model, thus implying that the regression results of the SEM model simulating air pollutants were better. thus, implying that the regression results of the SEM model simulating air pollutants were better. After comprehensive analysis and comparison, the SEM was selected for the spatial regression of the four atmospheric pollutants in this study. Therefore, this research established an SEM using the air pollutant concentrations of 37 Chinese cities as the dependent variables and the landscape pattern indices of the green spaces of each city as the independent variables [52,53].

$$\mathbf{Y} = \beta\_1 \mathbf{X}\_1 + \beta\_2 \mathbf{X}\_2 + \beta\_3 \mathbf{X}\_3 + \dots \dots \dots + \beta\_{\theta} \mathbf{X}\_{\theta} + \varepsilon \tag{1}$$

where Y denotes the dependent variable, i.e., pollutant concentration; Wy is the spatial weight matrix; ρ denotes the spatial regression coefficient of the spatial weight matrix WY, which is adapted to indicate the spatial interaction of air pollution; X1–X<sup>9</sup> represents the impact factors of landscape pattern indexes (PLAND, PD, LSI) for three types of urban green space: forest land, grassland, and farmland, respectively (The independent variables); β1–β<sup>9</sup> denote the regression coefficients of independent factors such as landscape pattern

indices; ε denotes the random error. This was ultimately utilized to establish the SEM spatial regression model between pollutant concentrations and landscape patterns in 2019.


**Table 3.** Information comparison of SEM and SAR models.

### **5. Conclusions**

Urban green space is an effective tool to improve air quality. Quantifying the relationship between urban green space landscape patterns and air pollution concentration is of great significance to promote the public health level and sustainable development of residents in high-density population areas. This study takes 37 garden cities with subtropical monsoon climate in China as the research object, selects the urban air quality monitoring data and land use remote sensing data in 2019, carries out regression analysis on urban air pollutant concentration and landscape pattern index through spatial regression model method, and explores the relationship between landscape pattern index and air pollutant concentration, According to the above regression analysis results, the landscape pattern index threshold with significant correlation with air pollutant concentration was explored. The specific conclusions are as follows:


The results of this study not only clarify the impact mechanism of the landscape pattern of urban green space on air quality but also propose a polynomial-based threshold effect estimation method, which provides quantitative reference and scientific basis for the optimization and updating of urban green space landscape patterns to promote public health.

**Author Contributions:** C.W.: Conceptualization, Methodology, Software; M.G.: Writing—Reviewing and Editing; J.J.: Data curation, Writing—Original draft preparation; Y.Y.: Visualization, Investigation; Y.R.: Software, Methodology; Y.W.: Validation; J.C.: Supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D Program of China (Funder: Jun Jin; Funding No.: 2020YFD1100405); the National Natural Science Foundation of China (Funder: Jiajie Cao; Funding No.: 32071832); the Humanity and Social Sciences Youth Foundation of the Ministry of

Education of China (Funder: Chengkang Wang; Funding No.: 18YJCZH167) and the Natural Science Foundation of Jiangsu Province (Funder: Chengkang Wang; Funding No.: BK20190751).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A**

**Table A1.** The effect of landscape patterns of urban green spaces on air pollution levels.


**Appendix B**



### **References**


**Jie Zhao 1,2, Kunlun Xiang <sup>3</sup> , Zhitao Wu <sup>1</sup> and Ziqiang Du 1,\***


**\*** Correspondence: duzq@sxu.edu.cn; Tel.: +86-25-8542-7231

**Abstract:** The distribution of global warming has been varying both diurnally and seasonally. Little is known about the spatiotemporal variations in the relationships between vegetation greenness and day- and night-time warming during the last decades. We investigated the global inter- and intra-annual responses of vegetation greenness to the diurnal asymmetric warming during the period of 1982–2015, using the normalized different vegetation index (NDVI, a robust proxy for vegetation greenness) obtained from the NOAA/AVHRR NDVI GIMMS3g dataset and the monthly average daily maximum (Tmax) and minimum temperature (Tmin) obtained from the gridded Climate Research Unit, University of East Anglia. Several findings were obtained: (1) The strength of the relationship between vegetation greenness and the diurnal temperature varied on inter-annual and seasonal timescales, indicating generally weakening warming effects on the vegetation activity across the global. (2) The decline in vegetation response to Tmax occurred mainly in the mid-latitudes of the world and in the high latitudes of the northern hemisphere, whereas the decline in the vegetation response to Tmin primarily concentrated in low latitudes. The percentage of areas with a significantly negative trend in the partial correlation coefficient between vegetation greenness and diurnal temperature was greater than that of the areas showing the significant positive trend. (3) The trends in the correlation between vegetation greenness and diurnal warming showed a complex spatial pattern: the majority of the study areas had undergone a significant declining strength in the vegetation greenness response to Tmax in all seasons and to Tmin in seasons except autumn. These findings are expected to have important implications for studying the diurnal asymmetry warming and its effect on the terrestrial ecosystem.

**Keywords:** varying response; diurnal warming; vegetation activity; NDVI

### **1. Introduction**

Surface vegetation cover is a core component of terrestrial ecosystems, as it plays an important role in connecting material circulation and energy flow in the atmosphere, hydrosphere, and soil circles [1–4]. Current satellite observations have revealed long-term greening and browning changes in the vegetation greenness (i.e., the normalized difference vegetation index, NDVI) across the global from the 1980s until the present time [2,5–8]. Such variations in the vegetation activity are closely related to climate change, particularly the increasing global temperatures over the last several decades [1,8–11], which can be characterized by the increase in both daytime and night-time temperature [12,13]. The link between global warming and vegetation activity is overwhelming. Recent studies have investigated the responses of satellite-derived vegetation productivity to temperatures during daytime and nighttime between regions and ecosystem [12,14–20]. However, most of the measurements were static and insufficient to clarify the varying responses of vegetation dynamics to climate warming [21].

The relationship between vegetation activity and temperature variations is not fixed; it constantly changes over time because of the interventions of other environmental and

**Citation:** Zhao, J.; Xiang, K.; Wu, Z.; Du, Z. Varying Responses of Vegetation Greenness to the Diurnal Warming across the Global. *Plants* **2022**, *11*, 2648. https://doi.org/ 10.3390/plants11192648

Academic Editors: Jie Gao, Weiwei Huang, Johan Gielis and Peijian Shi

Received: 13 September 2022 Accepted: 5 October 2022 Published: 8 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

anthropogenic factors, in the face of persistent global warming [22,23]. Therefore, the current focus is on the dynamic responses of vegetation greenness to temperature change over time. For instance, Andreu-Hayles et al. [24] explored the varying boreal forest responses to arctic environmental change near the Firth River, Alaska. Piao et al. [23] found a weakening relationship between inter-annual temperature variability and northern vegetation activity. Fu et al. [25] reported declining global warming effects on the phenology of spring leaf unfolding. Cong et al. [26] examined the varied responses of vegetation activity to climate change on the grasslands of the Tibetan Plateau. He et al. [21] identified obvious shifts in the relationship between vegetation growth and temperature in China's temperate desert and rainforest areas. Although these findings confirm the dynamic relationship between vegetation greenness and temperature change, the link between diurnal warming and vegetation activity, particularly the variations in their correlation over time, remains unclear.

The main objective of this study was to explore the varying relationships between vegetation greenness and day- and night-time warming on temporal and spatial scale, by simultaneously employing time-series satellite-derived vegetation NDVI data and gridded meteorological data. The understandings drawn from the findings are expected to have important implications for studying the diurnal asymmetry warming and its effect on the terrestrial ecosystem.

### **2. Results**

### *2.1. Trends of Correlations between Vegetation Greenness and Diurnal Warming*

### 2.1.1. Inter-Annual Changes in RNDVI-Tmax and RNDVI-Tmin

Recent changes in satellite-based vegetation greenness across different areas were closely related to global warming. Figure 1 shows the inter-annual variation in the partial correlation coefficients of the NDVI and diurnal temperature for the period of 1982–2015. In the mid-latitudes of the southern hemisphere, RNDVI-Tmax is approximately 0.323 for the first window of 1982–1998 and decreased to −0.436 for the last window of 1999–2015. In the high latitudes of the northern hemisphere, RNDVI-Tmin shows a significant increasing trend (β = 0.098, *p* < 0.01) in the first six windows and then a significant decreasing trend (β = −0.056, *p* < 0.01) in the remaining windows. Similarly, in the low latitudes of the southern hemisphere, RNDVI-Tmin shows a significant increasing trend (β = 0.039, *p* < 0.01) in the first seven windows, followed by a significant decreasing trend (β = −0.055, *p* < 0.01) in the remaining windows. In contrast, in the mid-latitudes of the southern hemisphere, RNDVI-Tmin is generally on an upward trend, increasing from 0.022 in the window of 1982–1998 to 0.737 in the window of 1998–2014. These findings indicate a notable shift in the responses of vegetation activity to diurnal temperatures over the last three decades.

Overall, RNDVI-Tmax exhibits only a significant downward trend in the mid-latitudes of the southern hemisphere (β = −0.051, *p* < 0.05). In contrast, RNDVI-Tmin shows a significant upward trend in the mid-latitudes of the southern hemisphere (β = 0.032, *p* < 0.05), but a significant downward trend in the low latitudes of the southern hemisphere (β = −0.021, *p* < 0.05). RNDVI-Tmax does not exhibit a significant trend in the northern hemisphere (*p* > 0.05), and RNDVI-Tmin exhibits a significant downward trend only in the high latitudes of the northern hemisphere (β = −0.021, *p* < 0.05).

**Figure 1.** Temporal variations in the partial correlation coefficients between mean annual NDVI and the diurnal temperature (Tmax and Tmin) for each 17-year moving window across latitudes intervals. (a, b, c, d, and e represents latitudes intervals at 60º~ 90ºN, 30º~ 60ºN, 0º~ 30ºN, −30º~ 0ºS, and −60º~ 0ºS, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998,…, 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients). **Figure 1.** Temporal variations in the partial correlation coefficients between mean annual NDVI and the diurnal temperature (Tmax and Tmin) for each 17-year moving window across latitudes intervals. (**a**–**e** represents latitudes intervals at 60◦~90◦ N, 30◦~60◦ N, 0◦~30◦ N, −30◦~0◦ S, and −60◦~0◦ S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998, . . . , 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients).

#### Overall, RNDVI-Tmax exhibits only a significant downward trend in the mid-latitudes of 2.1.2. Inter-Annual Changes in RNDVI-Tmax and RNDVI-Tmin

the southern hemisphere (β = −0.051, *p* < 0.05). In contrast, RNDVI-Tmin shows a significant upward trend in the mid-latitudes of the southern hemisphere (β = 0.032, *p* < 0.05), but a significant downward trend in the low latitudes of the southern hemisphere (β = −0.021, *p* < 0.05). RNDVI-Tmax does not exhibit a significant trend in the northern hemisphere (*p* > 0.05), and RNDVI-Tmin exhibits a significant downward trend only in the high latitudes of the north-Variations in the temperature on a seasonal scale may not correspond to similar changes on inter-annual scale. The photosynthesis of vegetation in different growth stages depends on the seasonal cycle of the temperature. Therefore, the sensitivity of vegetation greenness to day- and nighttime warming and variations in seasonality is likely to be highly complex.

ern hemisphere (β = −0.021, *p* < 0.05). 2.1.2. Inter-annual Changes in RNDVI-Tmax and RNDVI-Tmin Variations in the temperature on a seasonal scale may not correspond to similar changes on inter-annual scale. The photosynthesis of vegetation in different growth stages During spring, RNDVI-Tmax shows a significant downward trend at the mid-latitudes in both the northern hemisphere (β = −0.031, *p* < 0.05) and the southern hemisphere (β = −0.055, *p* < 0.05) (Figure 2). RNDVI-Tmin shows a significant downward trend (β = −0.040, *p* < 0.05) at the low latitudes in the northern hemisphere while a significant upward trend (β = 0.035, *p* < 0.05) at the mid-latitudes in the southern hemisphere (Figure 2).

depends on the seasonal cycle of the temperature. Therefore, the sensitivity of vegetation greenness to day- and nighttime warming and variations in seasonality is likely to be highly complex. During spring, RNDVI-Tmax shows a significant downward trend at the mid-latitudes in both the northern hemisphere (β = −0.031, *p* < 0.05) and the southern hemisphere (β= −0.055, *p* < 0.05) (Figure 2). RNDVI-Tmin shows a significant downward trend (β = −0.040, *p* < During summer, RNDVI-Tmax shows a significant downward trend not only at the low latitudes in the northern hemisphere (β = −0.020, *p* < 0.05), but also at the low and midlatitudes in the southern hemisphere (β = −0.029, *p* < 0.05) (Figure 3). RNDVI-Tmin shows a significant downward trend (β = −0.022, *p* < 0.05) at the high latitudes in the northern hemisphere, but a significant upward trend (β = 0.023, *p* < 0.05) at the mid-latitudes in the southern hemisphere (Figure 3).

0.05) at the low latitudes in the northern hemisphere while a significant upward trend (β=

0.035, *p* < 0.05) at the mid-latitudes in the southern hemisphere (Figure 2).

**Figure 2.** Temporal variations in the partial correlation coefficients between mean NDVI and the diurnal temperature (Tmax and Tmin) in spring for each 17-year moving window across latitudes intervals. (**a**–**e** represents latitudes intervals at 60°~90° N, 30°~60° N, 0°~30° N, −30°~0° S, and −60°~0° S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998,…, 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients.). During summer, RNDVI-Tmax shows a significant downward trend not only at the low **Figure 2.** Temporal variations in the partial correlation coefficients between mean NDVI and the diurnal temperature (Tmax and Tmin) in spring for each 17-year moving window across latitudes intervals. (**a**–**e** represents latitudes intervals at 60◦~90◦ N, 30◦~60◦ N, 0◦~30◦ N, −30◦~0◦ S, and −60◦~0◦ S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998, . . . , 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients). *Plants* **2022**, *11*, x FOR PEER REVIEW 5 of 14

**Figure 3.** Temporal variations in the partial correlation coefficients between mean NDVI and the **Figure 3.** Temporal variations in the partial correlation coefficients between mean NDVI and the diurnal temperature (Tmax and Tmin) in summer for each 17-year moving window across latitudes in-

diurnal temperature (Tmax and Tmin) in summer for each 17-year moving window across latitudes intervals. (**a**–**e** represents latitudes intervals at 60°~90° N, 30°~60° N, 0°~30° N, −30°~0° S, and −60°~0° S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from

During autumn, RNDVI-Tmax shows a significant upward trend at the low-latitudes (β= 0.057, *p* < 0.01) and the mid-latitudes (β = 0.021, *p* < 0.01) in the northern hemisphere. (Fig.4) RNDVI-Tmin shows a significant downward trend at the mid-latitudes (β = –0.017, *p* < 0.05) and the low latitudes (β = –0.048, *p* < 0.01) in the northern hemisphere (Figure 4). However, in the southern hemisphere, we did not find any statistically significant variations in RNDVI-

Tmax and RNDVI-Tmin.

tervals. (**a**–**e** represents latitudes intervals at 60◦~90◦ N, 30◦~60◦ N, 0◦~30◦ N, −30◦~0◦ S, and −60◦~0◦ S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998, . . . , 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients).

During autumn, RNDVI-Tmax shows a significant upward trend at the low-latitudes (β = 0.057, *p* < 0.01) and the mid-latitudes (β = 0.021, *p* < 0.01) in the northern hemisphere. (Figure 4) RNDVI-Tmin shows a significant downward trend at the mid-latitudes (β = −0.017, *p* < 0.05) and the low latitudes (β = −0.048, *p* < 0.01) in the northern hemisphere (Figure 4). However, in the southern hemisphere, we did not find any statistically significant variations in RNDVI-Tmax and RNDVI-Tmin. *Plants* **2022**, *11*, x FOR PEER REVIEW 6 of 14

**Figure 4.** Temporal variations in the partial correlation coefficients between mean NDVI and the diurnal temperature (Tmax and Tmin) in autumn for each 17-year moving window across latitudes intervals. (**a**–**e** represents latitudes intervals at 60°~90° N, 30°~60° N, 0°~30° N, −30°~0° S, and −60°~0° S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998,…, 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients.). **Figure 4.** Temporal variations in the partial correlation coefficients between mean NDVI and the diurnal temperature (Tmax and Tmin) in autumn for each 17-year moving window across latitudes intervals. (**a**–**e** represents latitudes intervals at 60◦~90◦ N, 30◦~60◦ N, 0◦~30◦ N, −30◦~0◦ S, and −60◦~0◦ S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998, . . . , 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients).

During winter, RNDVI-Tmax shows a significant upward trend at the mid-latitudes in the northern hemisphere (β = 0.041, *p* < 0.01) and the low latitudes in the southern hemisphere (β = 0.016, *p* < 0.05), but a significant downward trend (β = –0.025, *p* < 0.01) at the midlatitudes in the southern hemisphere (Figure 5). RNDVI-Tmin shows a significant downward trend at the low (β = –0.010, *p* < 0.05) and the mid-latitudes (β = –0.036, *p* < 0.01) in the northern hemisphere and at the low latitudes (β= –0.025, *p* < 0.01) in the southern hemisphere, but a significant upward trend at the mid-latitudes (β = 0.027, *p* < 0.01) in the During winter, RNDVI-Tmax shows a significant upward trend at the mid-latitudes in the northern hemisphere (β = 0.041, *p* < 0.01) and the low latitudes in the southern hemisphere (β = 0.016, *p* < 0.05), but a significant downward trend (β = −0.025, *p* < 0.01) at the mid-latitudes in the southern hemisphere (Figure 5). RNDVI-Tmin shows a significant downward trend at the low (β = −0.010, *p* < 0.05) and the mid-latitudes (β = −0.036, *p* < 0.01) in the northern hemisphere and at the low latitudes (β = −0.025, *p* < 0.01) in the southern hemisphere, but a significant upward trend at the mid-latitudes (β = 0.027, *p* < 0.01) in the southern hemisphere (Figure 5). Although we tried to characterize the sensitivity of vegeta-

> southern hemisphere (Figure 5). Although we tried to characterize the sensitivity of vegetation greenness for each season, the possible mechanism leading to the seasonal divergence in the response of vegetation greenness to diurnal warming needs to be studied

further.

tion greenness for each season, the possible mechanism leading to the seasonal divergence in the response of vegetation greenness to diurnal warming needs to be studied further. *Plants* **2022**, *11*, x FOR PEER REVIEW 7 of 14

> **Figure 5.** Temporal variations in the partial correlation coefficients between mean NDVI and the diurnal temperature (Tmax and Tmin) in winter for each 17-year moving window across latitudes intervals. (**a**–**e** represents latitudes intervals at 60°~90° N, 30°~60° N, 0°~30° N, −30°~0° S, and −60°~0° S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998,…, 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients.). **Figure 5.** Temporal variations in the partial correlation coefficients between mean NDVI and the diurnal temperature (Tmax and Tmin) in winter for each 17-year moving window across latitudes intervals. (**a**–**e** represents latitudes intervals at 60◦~90◦ N, 30◦~60◦ N, 0◦~30◦ N, −30◦~0◦ S, and −60◦~0◦ S, respectively (N and S indicates the Northern and southern hemisphere, respectively). The x axis is the last year of the 17-year moving-window (for example, 1998 stands for a moving-window from 1982 to 1998, . . . , 2015 stands for a moving-window from 1999 to 2015). The Y axis is the partial correlation coefficients).

#### *2.2. Spatial Patterns of the Trends in the Correlations between Vegetation Greenness and 2.2. Spatial Patterns of the Trends in the Correlations between Vegetation Greenness and Diurnal Temperatures*

#### *Diurnal Temperatures*  2.2.1. Inter-Annual Patterns of RNDVI-Tmax and RNDVI-Tmin

.

2.2.1. Inter-annual Patterns of RNDVI-Tmax and RNDVI-Tmin Figure 6 shows a summary of the spatial variations in the correlations between vegetation activity and diurnal temperatures over 30 years. The areas with a significant increase in RNDVI-Tmax account for 26.01% of the total vegetated areas, mainly located in the southeastern part of Oceania, the Ganges Plain, the western plains of Eastern Europe, and the Amazonian Plain. The areas with a significant decrease in RNDVI-Tmax account for 35.65% of the total vegetated areas, mainly distributed in the northwestern part of Oceania, southern South America, eastern Africa, and the cold temperate zones of the northern hemisphere. About 28.62% of the total vegetated areas show a significant upward trend in RNDVI-Tmin, these areas are mainly located in southern South America, southern Africa, Oceania, and the northwestern boreal regions of North America. Over 32.49% of the total vegetated areas shows a significant downward trend in RNDVI-Tmin, these areas are mainly distributed in eastern Oceania, central Africa, and the tropical regions in the northern hemisphere. In particular, the number of pixels with a significant decrease in RNDVI-Tmax were Figure 6 shows a summary of the spatial variations in the correlations between vegetation activity and diurnal temperatures over 30 years. The areas with a significant increase in RNDVI-Tmax account for 26.01% of the total vegetated areas, mainly located in the southeastern part of Oceania, the Ganges Plain, the western plains of Eastern Europe, and the Amazonian Plain. The areas with a significant decrease in RNDVI-Tmax account for 35.65% of the total vegetated areas, mainly distributed in the northwestern part of Oceania, southern South America, eastern Africa, and the cold temperate zones of the northern hemisphere. About 28.62% of the total vegetated areas show a significant upward trend in RNDVI-Tmin, these areas are mainly located in southern South America, southern Africa, Oceania, and the northwestern boreal regions of North America. Over 32.49% of the total vegetated areas shows a significant downward trend in RNDVI-Tmin, these areas are mainly distributed in eastern Oceania, central Africa, and the tropical regions in the northern hemisphere. In particular, the number of pixels with a significant decrease in RNDVI-Tmax were greater than those with a significant increase in RNDVI-Tmax in each latitude interval.

greater than those with a significant increase in RNDVI-Tmax in each latitude interval.

This difference in the number of pixels between RNDVI-Tmax and RNDVI-Tmin was most pronounced at the mid-latitudes in the southern hemisphere and the mid- and high

RNDVI-Tmin shows the opposite trend.

latitudes in the northern hemisphere. At the mid- and high latitudes in the northern hemisphere and at the mid-latitudes in the southern hemisphere, the proportion of RNDVI-Tmin showing a significant upward trend is slightly higher than that of RNDVI-Tmin showing a significant downward trend. At low latitudes, however, the percentage of areas where RNDVI-Tmin shows a significant downward trend is much higher than that of areas where

**Figure 6.** The response of vegetation greenness to the diurnal temperature. (**a**. spatial distribution of the temporal trend of the partial coefficients between mean annual NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between mean annual NDVI and Tmin. Supplementary Table S1). **Figure 6.** The response of vegetation greenness to the diurnal temperature. (**a**. spatial distribution of the temporal trend of the partial coefficients between mean annual NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between mean annual NDVI and Tmin. Supplementary Table S1).

Overall, regions with a declining correlation between vegetation greenness and diurnal temperatures were more than those with an increasing correlation. This implies that worldwide vegetation activity on a larger scale became less sensitive to variations in the day- and night temperature with global warming during the last several decades. 2.2.2. Intra-annual Patterns of RNDVI-Tmax and RNDVI-Tmin The seasonal variation in the vegetation greenness response to global warming may increase the complexity of the intra-annual patterns of RNDVI-Tmax and RNDVI-Tmin. For spring, This difference in the number of pixels between RNDVI-Tmax and RNDVI-Tmin was most pronounced at the mid-latitudes in the southern hemisphere and the mid- and high latitudes in the northern hemisphere. At the mid- and high latitudes in the northern hemisphere and at the mid-latitudes in the southern hemisphere, the proportion of RNDVI-Tmin showing a significant upward trend is slightly higher than that of RNDVI-Tmin showing a significant downward trend. At low latitudes, however, the percentage of areas where RNDVI-Tmin shows a significant downward trend is much higher than that of areas where RNDVI-Tmin shows the opposite trend.

the areas with positive trends in RNDVI-Tmax accounts for 25.66% of the total vegetated areas, mostly in the northeastern part of South America, the eastern plains of Eastern Europe, the Iberian Peninsula, the southeastern parts of Australia, and the Ganges Plain (Figure 7). The areas with negative trends in RNDVI-Tmax accounts for 32.43% of the total vegetated Overall, regions with a declining correlation between vegetation greenness and diurnal temperatures were more than those with an increasing correlation. This implies that worldwide vegetation activity on a larger scale became less sensitive to variations in the day- and night temperature with global warming during the last several decades.

#### areas, mainly distributed in Western Australia, the Katanga Plateau, and the cold temper-2.2.2. Intra-Annual Patterns of RNDVI-Tmax and RNDVI-Tmin

ate zones in the northern hemisphere. Similarly, about 24.09% of the vegetated areas show positive trends in RNDVI-Tmin, these areas are mainly located in northern North America, southern South America, southern Europe, the Indian Peninsula, and southwestern Australia. Over 32.11% of the vegetated areas show negative trends in RNDVI-Tmin, these areas are mainly distributed in northern South America, central Africa, southeastern Australia, and the temperate regions of the northern hemisphere. The seasonal variation in the vegetation greenness response to global warming may increase the complexity of the intra-annual patterns of RNDVI-Tmax and RNDVI-Tmin. For spring, the areas with positive trends in RNDVI-Tmax accounts for 25.66% of the total vegetated areas, mostly in the northeastern part of South America, the eastern plains of Eastern Europe, the Iberian Peninsula, the southeastern parts of Australia, and the Ganges Plain (Figure 7). The areas with negative trends in RNDVI-Tmax accounts for 32.43% of the total vegetated areas, mainly distributed in Western Australia, the Katanga Plateau, and the cold temperate zones in the northern hemisphere. Similarly, about 24.09% of the vegetated areas show positive trends in RNDVI-Tmin, these areas are mainly located in northern North America, southern South America, southern Europe, the Indian Peninsula, and southwestern Australia. Over 32.11% of the vegetated areas show negative trends in RNDVI-Tmin, these areas are mainly distributed in northern South America, central Africa, southeastern Australia, and the temperate regions of the northern hemisphere.

For summer, RNDVI-Tmax shows significant positive trends in 29.80% of the vegetated areas, located mainly in the northern and eastern North America, southern South America, Armenian plateau, Tulan lowland, Kazakh hills, and Ganges plain (Figure 8). RNDVI-Tmax shows significant negative trends in 31.78% of the vegetated areas, located mainly in southern North America, Central European Plains, Scandinavia, Northeast Asia, India Peninsula, and Indochina. Similarly, RNDVI-Tmin shows significant positive trends in 26.68% of the vegetated areas, mainly in the southern Patagonia Plateau of South America, the Yukon Plateau in North America, the south-central United States, the northwestern part of Asia, the Indo-China Peninsula, and Southern Australia. RNDVI-Tmin shows significant negative trends in 34.71% of the vegetated areas, located mainly in northern South America,

southern Africa, northwestern Australia, and the cold temperate zones of the northern hemisphere (Figure 8). *Plants* **2022**, *11*, x FOR PEER REVIEW 9 of 14 **Figure 7.** The response of vegetation greenness to the diurnal temperature in spring. (**a**. spatial dis-

tribution of the temporal trend of the partial coefficients between spring NDVI and Tmax. **b**. spatial

*Plants* **2022**, *11*, x FOR PEER REVIEW 9 of 14

**Figure 7.** The response of vegetation greenness to the diurnal temperature in spring. (**a**. spatial distribution of the temporal trend of the partial coefficients between spring NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between spring NDVI and Tmin. Supplementary Table S2). **Figure 7.** The response of vegetation greenness to the diurnal temperature in spring. (**a**. spatial distribution of the temporal trend of the partial coefficients between spring NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between spring NDVI and Tmin. Supplementary Table S2). of Asia, the Indo-China Peninsula, and Southern Australia. RNDVI-Tmin shows significant negative trends in 34.71% of the vegetated areas, located mainly in northern South America, southern Africa, northwestern Australia, and the cold temperate zones of the northern hemisphere (Figure 8).

hemisphere (Figure 8). **Figure 8.** The response of vegetation greenness to the diurnal temperature in summer. (**a**. spatial distribution of the temporal trend of the partial coefficients between summer NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between summer NDVI and Tmin. Supplementary Table S3). **Figure 8.** The response of vegetation greenness to the diurnal temperature in summer. (**a**. spatial distribution of the temporal trend of the partial coefficients between summer NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between summer NDVI and Tmin. Supplementary Table S3).

**Figure 8.** The response of vegetation greenness to the diurnal temperature in summer. (**a**. spatial distribution of the temporal trend of the partial coefficients between summer NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between summer NDVI and Tmin. For autumn, about 27.43% of the total vegetated areas show significant positive trends in RNDVI-Tmax, these areas are mainly distributed in central Eurasia, central and southern North America, northern Africa, and central Australia (Figure 9). Over 34.07% of the total vegetated areas show significant negative trends in RNDVI-Tmax, these areas are mainly distributed in northern Eurasia, central North America, southern Africa, western Africa, southern South America, and northeastern Australia. In contrast, over 32.67% of the total vegetated areas show significant positive trends in RNDVI-Tmin, these areas are widely distributed in southern Africa, northern Australia, and the cold temperate zones of the northern hemisphere. About 28.76% of the total vegetated areas show significant For autumn, about 27.43% of the total vegetated areas show significant positive trends in RNDVI-Tmax, these areas are mainly distributed in central Eurasia, central and southern North America, northern Africa, and central Australia (Figure 9). Over 34.07% of the total vegetated areas show significant negative trends in RNDVI-Tmax, these areas are mainly distributed in northern Eurasia, central North America, southern Africa, western Africa, southern South America, and northeastern Australia. In contrast, over 32.67% of the total vegetated areas show significant positive trends in RNDVI-Tmin, these areas are widely distributed in southern Africa, northern Australia, and the cold temperate zones of the northern hemisphere. About 28.76% of the total vegetated areas show significant negative trends in RNDVI-Tmin, these areas are mainly distributed in southern North America, northwestern South America, central Asia, southern Europe, central Africa, and southern Europe (Figure 9).

Supplementary Table S3). For autumn, about 27.43% of the total vegetated areas show significant positive trends in RNDVI-Tmax, these areas are mainly distributed in central Eurasia, central and southern North America, northern Africa, and central Australia (Figure 9). Over 34.07% of the total vegetated areas show significant negative trends in RNDVI-Tmax, these areas are mainly distributed in northern Eurasia, central North America, southern Africa, western Africa, southern South America, and northeastern Australia. In contrast, over 32.67% of For winter, the positive trends in RNDVI-Tmax are found in 25.06% of the total vegetated areas, located mainly in central Asia, eastern Australia, Eastern Europe, and southern China (Figure 10). The negative trends in RNDVI-Tmax are found in 25.74% of the total vegetated areas, mainly distributed in Western Europe, central and southern North America, central South America, eastern Africa, Western Australia, Indian Peninsula, and Indochina. Similarly, the positive trends in RNDVI-Tmin are found in 23.32% of the total vegetated areas, located mainly in central North America, Western Europe, southwestern Australia, and

> the total vegetated areas show significant positive trends in RNDVI-Tmin, these areas are widely distributed in southern Africa, northern Australia, and the cold temperate zones

the north-central La Plata Plain in South America. The negative trends in RNDVI-Tmin are found in 26.71% of the total vegetated areas, mainly distributed in central Asia, central and eastern Australia, and Eastern Europe (Figure 10). southern Europe (Figure 9).

negative trends in RNDVI-Tmin, these areas are mainly distributed in southern North America, northwestern South America, central Asia, southern Europe, central Africa, and

negative trends in RNDVI-Tmin, these areas are mainly distributed in southern North America, northwestern South America, central Asia, southern Europe, central Africa, and

*Plants* **2022**, *11*, x FOR PEER REVIEW 10 of 14

*Plants* **2022**, *11*, x FOR PEER REVIEW 10 of 14

southern Europe (Figure 9).

**Figure 9.** The response of vegetation greenness to the diurnal temperature in autumn. (**a**. spatial distribution of the temporal trend of the partial coefficients between autumn NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between autumn NDVI and Tmin. Supplementary Table S4). **Figure 9.** The response of vegetation greenness to the diurnal temperature in autumn. (**a**. spatial distribution of the temporal trend of the partial coefficients between autumn NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between autumn NDVI and Tmin. Supplementary Table S4). located mainly in central North America, Western Europe, southwestern Australia, and the north-central La Plata Plain in South America. The negative trends in RNDVI-Tmin are found in 26.71% of the total vegetated areas, mainly distributed in central Asia, central and eastern Australia, and Eastern Europe (Figure 10).

**Figure 10.** The response of vegetation greenness to the diurnal temperature in winter. (**a**. spatial distribution of the temporal trend of the partial coefficients between winter NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between winter NDVI and Tmin. Supplementary Table S5). **Figure 10.** The response of vegetation greenness to the diurnal temperature in winter. (**a**. spatial distribution of the temporal trend of the partial coefficients between winter NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between winter NDVI and Tmin. Supplementary Table S5).

### **3. Discussion**

**3. Discussion** 

**Figure 10.** The response of vegetation greenness to the diurnal temperature in winter. (**a**. spatial distribution of the temporal trend of the partial coefficients between winter NDVI and Tmax. **b**. spatial distribution of the temporal trend of the partial coefficients between winter NDVI and Tmin. Supplementary Table S5). The global climate has undergone continuous warming since the 19th century. Such warming has been characterized by a faster change in the global land surface during the nighttime than during the daytime [12,27,28]. However, the current understanding of the sensitivity of vegetation activity to daytime and nighttime temperatures (not to mean The global climate has undergone continuous warming since the 19th century. Such warming has been characterized by a faster change in the global land surface during the nighttime than during the daytime [12,27,28]. However, the current understanding of the sensitivity of vegetation activity to daytime and nighttime temperatures (not to mean temperature) during the last decades is still poor. Therefore, we applied the movingwindow-based partial correlation analysis and ordinary least squares linear regression method within a moving 17-year window to the GIMMS NDVI3g dataset and gridded meteorological data covering the global vegetated areas to reveal the variations in the response of vegetation greenness to the day- and nighttime warming during the period of 1982–2015.

**3. Discussion**  The global climate has undergone continuous warming since the 19th century. Such warming has been characterized by a faster change in the global land surface during the nighttime than during the daytime [12,27,28]. However, the current understanding of the sensitivity of vegetation activity to daytime and nighttime temperatures (not to mean We reported a wide reduction in the vegetation NDVI and diurnal temperature sensitivity since the early 1980s across the global. We also showed that the strength of the relationship between vegetation greenness and diurnal temperature varied among seasons and not just at inter-annual intervals. We further confirmed declining warming effects on vegetation activity, in line with some earlier studies conducted from local to global scale [21–26].

We tried to reveal the dynamic responses of vegetation greenness to day and night warming and explore the difference in the responses of vegetation to such diurnal warming. However, some uncertainties in this study still constrained our understanding of these responses. For example, the changing response could be related to any one, or a combination, of various factors, such as increasing drought stress [29–31], rising CO2 fertilization [22,23], fire disturbances [32,33], global dimming [29], and human-associated activities [21]. These factors were not considered in our analyses because of the lack of sufficient records. Regardless of the cause, our findings on the change in the response of vegetation greenness to diurnal warming have important implications for studies on climate change related to the past as well as the future. Further studies are needed, combining process studies (e.g., ecosystem warming experiments) with multi-scale data, to clarify the mechanisms of the observed decline in the correlation between the diurnal warming and vegetation activity. Our study mainly examined the immediate responses of terrestrial vegetation activity to asymmetric warming. In fact, such asymmetric warming has a non-uniform time-lag impact on terrestrial vegetation growth [34]. Previous studies have also found different responses of vegetation activity to warming across terrestrial ecosystems [12,23]. We did not pay enough attention to this difference in the biome involved. Although the combined effects of hydrothermal regime, ecosystem types, and the time-lag effects increase the uncertainties associated with the weakening relationship between vegetation greenness and current asymmetric warming during day- and night-time, our findings demonstrate the need to consider asymmetric warming and its uneven effect rather than variations in the annual mean temperature when modeling the future behavior of vegetation activity. Hence, significant efforts are required to examine the varying responses of vegetation activity to climate change in such a dynamic climate system in light of the substantial spatial and temporal variations in the diurnal warming and its effects.

### **4. Materials and Methods**

### *4.1. Datasets*

The third-generation NDVI data for the period of 1982–2015 was produced by the Global Inventory Modelling and Mapping Studies group (GIMMS NDVI3g) from the NOAA/AVHRR series satellite imageries, with a spatial resolution of ~0.083º and fortnightly interval. To further minimize the atmospheric noise, the value composites method was used to reconstruct the monthly and yearly NDVI datasets [35]. These NDVI data were aggregated to 0.5º×0.5º to match the resolution of Meteorological data. Pixels with the mean NDVI (annual and monthly) below 0.1 served as non-vegetation areas and were employed as a mask [5,26]. The GIMMS NDVI3g can be used as a proxy for vegetation greenness and has been widely applied to explore the diverse spatiotemporal-scale vegetation activity [5,12,23,26,34].

Climatological data for the same period, including the monthly average daily maximum (Tmax) and minimum temperature (Tmin) and precipitation, with a spatial resolution of ~ 0.5◦ × 0.5◦ were obtained from the gridded Climate Research Unit, University of East Anglia (CRU TS 4.01). The CRU TS 4.01 data were produced using angular-distance weighting interpolation. These meteorological materials have been well used in studies on the relationship between global and regional vegetation cover and climate change [7,12,23,34,36].

### *4.2. Methods*

### 4.2.1. Moving-Window-Based Partial Correlation Analysis

The relationships between vegetation NDVI and temperature were examined using partial correlation analysis (the partial correlation coefficient of NDVI and day- and nighttime temperature were denoted by RNDVI-Tmax and RNDVI-Tmin, respectively), to eliminate the covariate effects of other factors [12,37–39]. In this case, for instance, the partial correlation coefficient between NDVI and Tmax was calculated when conditioning the precipitation and Tmin; and that between NDVI and Tmin was calculated when conditioning the precipitation and Tmax. The partial correlations of the NDVI with temperature were

calculated in the first 17-year window (1982–1998). The window was then moved forward by 1 year to represent the second 17-year window (1983–1999), and so on until the last 17-year window (1999–2015). Thus, we obtained the partial correlation coefficients within all 17-year windows, which were applied to describe the responses of vegetation greenness to day- and nighttime warming. The NDVI and relevant variables were detrended by a linear detrending method prior to the partial correlation analysis [12,17].

### 4.2.2. Ordinary Least Squares Linear Regression

To explore the dynamics in the correlation between global vegetation greenness and day- and night temperature, the Ordinary least-squares linear regression analysis method was used to calculate the slope of the partial correlation coefficients of NDVI and Tmax and Tmin within each 17-year window over the entire temporal domain. The partial correlation coefficient and P-value were used to determine the relationship and its significance. The slope of the linear regression represents the changing trend in the correlations between the vegetation greenness and diurnal warming during the period 1982–2015. A positive slope indicates a positive growth response to temperature, while a negative slope indicates a weakening relationship with diurnal warming.

$$y\_{(t)} = \alpha + \beta t + \varepsilon \tag{1}$$

Here, α and β were the fitted intercept and slope, respectively. A two-tailed *t*-test was used to examine the significance level of the trend. Thus, we described the change in vegetation greenness responses to the diurnal warming from 1982 to 2015.

### **5. Conclusions**

In the most recent period (i.e., 1982–2015), most parts of the world exhibited an overall significant change in the response of vegetation greenness to diurnal warming. On an interannual timescale, the percentage of areas showing a significantly negative trend in the partial correlation coefficient between vegetation greenness and day- and nighttime temperatures was greater than that of areas showing a significant positive trend. For daytime temperatures, this decline in the vegetation response occurred mainly at the midlatitudes of the world and at the high latitudes of the northern hemisphere. For nighttime temperatures, the decline was primarily concentrated at the low-latitudes across the global. Our study also showed that the strength of the associations between vegetation greenness and diurnal warming varied among seasons. Although the trends in the correlation between vegetation greenness and diurnal warming showed a complex spatial pattern, most of the study areas had undergone a significant declining strength in the vegetation greenness response to daytime temperatures in all seasons and to nighttime temperatures in seasons except autumn. This further implies that the recent reduction in the strength of the response of vegetation greenness to diurnal temperature fluctuations became widespread throughout the world as such as the inter-annual tendency. These findings shew new light on our understanding of the response of vegetation greenness to global warming. We demonstrate the need to carefully consider the asymmetric diurnal warming when unraveling the influence of global warming on terrestrial ecosystem.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11192648/s1, Table S1: Slopes of the partial correlation coefficients between mean annual NDVI and the duirnal temperature (Tmax and Tmin) for each 17-year moving window across latitudes intervals; Table S2: Slopes of the partial correlation coefficients between mean NDVI and the duirnal temperature (Tmax and Tmin) in spring for each 17-year moving window across latitudes intervals; Table S3: Slopes of the partial correlation coefficients between mean NDVI and the duirnal temperature (Tmax and Tmin) in summer for each 17-year moving window across latitudes intervals; Table S4: Slopes of the partial correlation coefficients between mean NDVI and the duirnal temperature (Tmax and Tmin) in autumn for each 17-year moving window across latitudes intervals; Table S5: Slopes of the partial correlation coefficients between mean NDVI

and the duirnal temperature (Tmax and Tmin) in winter for each 17-year moving window across latitudes intervals.

**Author Contributions:** Methodology, J.Z., and Z.D.; formal analysis, J.Z. and Z.W.; writing—original draft preparation, Z.D. and K.X.; writing—review and editing, Z.D. and Z.W.; funding acquisition, Z.D. and Z.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** Z.D. and Z.W. were funded by the National Natural Science Foundation of China (grant number: U1810101, 41161066, 41977412, and 418711934), and Shanxi Province Disciplines Group Construction Project of Service to Industrial Innovation: Ecological Remediation of Soil Pollution Disciplines Group (grant number: 20181401).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in the present work have been listed in the Supplementary Materials.

**Acknowledgments:** The authors thank the editor, the associate editor, and anonymous reviewers for their helpful and constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Xiao Zheng 1,2, Karl J. Niklas 3,\*, David A. Ratkowsky <sup>4</sup> , Yabing Jiao <sup>2</sup> , Hui Ding <sup>1</sup> and Peijian Shi 2,\***


**Abstract:** Leaf shape and size can vary between hybrids and their parents. However, this has seldom been quantitatively tested. *Photinia* × *fraseri* is an important landscaping plant in East Asia as a hybrid between evergreen shrubs *P. glabra* and *P. serratifolia*. Its leaf shape looks like that of *P. serratifolia*. To investigate leaf shape, we used a general equation for calculating the leaf area (*A*) of broad-leaved plants, which assumes a proportional relationship between *A* and product of lamina length (*L*) and width (*W*). The proportionality coefficient (which is referred to as the Montgomery parameter) serves as a quantitative indicator of leaf shape, because it reflects the proportion of leaf area *A* to the area of a rectangle with *L* and *W* as its side lengths. The ratio of *L* to *W*, and the ellipticalness index were also used to quantify the complexity of leaf shape for elliptical leaves. A total of >4000 leaves from *P.* × *fraseri* and *P. serratifolia* (with >2000 leaves for each taxon) collected on a monthly basis was used to examine: (i) whether there is a significant difference in leaf shape between the two taxa, and (ii) whether there is a monotonic or parabolic trend in leaf shape across leaf ages. There was a significant difference in leaf shape between the two taxa (*p* < 0.05). Although there were significant differences in leaf shape on a monthly basis, the variation in leaf shape over time was not large, i.e., leaf shape was relatively stable over time for both taxa. However, the leaf shape of the hybrid was significantly different from its parent *P. serratifolia*, which has wider and more elliptical leaves than the hybrid. This work demonstrates that variations in leaf shape resulting from hybridization can be rigorously quantified and compared among species and their hybrids. In addition, this work shows that leaf shape does not changes as a function of age either before or after the full expansion of the lamina.

**Keywords:** leaf ellipticalness index; leaf length; leaf width; Montgomery equation; Montgomery parameter

### **1. Introduction**

Leaves are the primary photosynthetic organs of the majority of land plants. Prior research has shown that leaf area and structure are correlated with photosynthetic rates [1–4]. It is important therefore to accurately calculate leaf area. Prior studies have shown that leaf area follows a general equation called the Montgomery equation (denoted henceforth as ME) for many broad-leaved species with different leaf shapes. This equation assumes a simple proportional relationship between leaf area (*A*) and the product of leaf length (*L*) and width (*W*) [5–13]. Leaf age and area have been demonstrated to have little influence on the validity of ME in calculating leaf area, although leaf age and size can influence to some degree the proportionality coefficient of ME among different leaf age or size groups [13,14].

Leaf shape and area are adaptively correlated and therefore can vary as a function of local climatic conditions even among conspecifics [15–18]. In previous studies, the

**Citation:** Zheng, X.; Niklas, K.J.; Ratkowsky, D.A.; Jiao, Y.; Ding, H.; Shi, P. Comparison of Leaf Shape between a *Photinia* Hybrid and One of Its Parents. *Plants* **2022**, *11*, 2370. https://doi.org/10.3390/ plants11182370

Academic Editor: Yasutomo Hoshika

Received: 12 August 2022 Accepted: 8 September 2022 Published: 11 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

leaf roundness index [19–22] has been widely used to measure leaf shape. However, this index is more suitable for leaves whose length and width are nearly equal. Consequently, Li et al. [23] proposed the 'leaf ellipticalness index' for elliptical and oblong leaves. For an elliptical leaf, the leaf ellipticalness index approximates 1. Additionally, in spite of the simplicity of using the leaf width/length ratio to quantify leaf shape, Shi et al. [24] have shown that the leaf width/length ratio is significantly positively correlated with the fractal dimension of the lamina perimeter based on a large sample of leaves differing in size and shape among nine Magnoliaceae species. al. [23] proposed the 'leaf ellipticalness index' for elliptical and oblong leaves. For an elliptical leaf, the leaf ellipticalness index approximates 1. Additionally, in spite of the simplicity of using the leaf width/length ratio to quantify leaf shape, Shi et al. [24] have shown that the leaf width/length ratio is significantly positively correlated with the fractal dimension of the lamina perimeter based on a large sample of leaves differing in size and shape among nine Magnoliaceae species. Hybridization across closely related plant species can also result in significant varia-

influence to some degree the proportionality coefficient of ME among different leaf age or

Leaf shape and area are adaptively correlated and therefore can vary as a function of local climatic conditions even among conspecifics [15–18]. In previous studies, the leaf roundness index [19–22] has been widely used to measure leaf shape. However, this index is more suitable for leaves whose length and width are nearly equal. Consequently, Li et

Hybridization across closely related plant species can also result in significant variation in leaf shape. However, the effect of hybridization on leaf size and shape has seldom been quantitatively examined, particularly as leaves mature during the growing season. In order to examine this phenomenology, *Photinia* × *fraseri*, which is a hybrid between the evergreen shrubs *P. glabra* and *P. serratifolia*, and one of its parents *P. serratifolia*, were examined to determine (i) whether there is a significant difference in leaf shape between *P.* × *fraseri* and *P. serratifolia*, and (ii) whether there is a monotonic or parabolic trend in leaf shape as a function of age. Because both the hybrid and one of its parents have elliptical or obovate leaves (Figure 1), we used the ratio of leaf width to length and the leaf ellipticalness index to quantify leaf shape in addition to the ME. tion in leaf shape. However, the effect of hybridization on leaf size and shape has seldom been quantitatively examined, particularly as leaves mature during the growing season. In order to examine this phenomenology, *Photinia × fraseri*, which is a hybrid between the evergreen shrubs *P. glabra* and *P. serratifolia*, and one of its parents *P. serratifolia*, were examined to determine (i) whether there is a significant difference in leaf shape between *P. × fraseri* and *P. serratifolia*, and (ii) whether there is a monotonic or parabolic trend in leaf shape as a function of age. Because both the hybrid and one of its parents have elliptical or obovate leaves (Figure 1), we used the ratio of leaf width to length and the leaf ellipticalness index to quantify leaf shape in addition to the ME.

*Plants* **2022**, *11*, x FOR PEER REVIEW 2 of 9

size groups [13,14].

**Figure 1.** Representative leaves for the two *Photinia* taxa (the hybrid *P. × fraseri* and one of its parents *P. serratifolia*) at different times in the growing season during 2021. **Figure 1.** Representative leaves for the two *Photinia* taxa (the hybrid *P.* × *fraseri* and one of its parents *P. serratifolia*) at different times in the growing season during 2021.

### **2. Materials and Methods 2. Materials and Methods**

### *2.1. Plant Materials*

*2.1. Plant Materials*  We marked newly emerging leaves on 36 trees of *P. × fraseri* (the hybrid) and 3 trees of *P. serratifolia* (one of the two parents) growing at the Nanjing Forestry University campus (118°48′35″ E, 32°4′67″ N) from late February to early March 2021. Representative specimens of the second parent (*P. glabra*) were not found for study on the University campus, which was under quarantine. Each month, 320–380 leaves were randomly sampled from three individual trees of *P. × fraseri* from mid-March to mid-August 2021 and 320–350 leaves were randomly sampled from three individual trees of *P. serratifolia* from mid-April to mid-November 2021. Figure 1 provides images of representative leaves across the different months for the two taxa. We had planned to sample the leaves of *P. serratifolia* from spring to autumn. However, the marked branches in early March were We marked newly emerging leaves on 36 trees of *P.* × *fraseri* (the hybrid) and 3 trees of *P. serratifolia* (one of the two parents) growing at the Nanjing Forestry University campus (118◦4803500 E, 32◦4 06700 N) from late February to early March 2021. Representative specimens of the second parent (*P. glabra*) were not found for study on the University campus, which was under quarantine. Each month, 320–380 leaves were randomly sampled from three individual trees of *P.* × *fraseri* from mid-March to mid-August 2021 and 320–350 leaves were randomly sampled from three individual trees of *P. serratifolia* from mid-April to mid-November 2021. Figure 1 provides images of representative leaves across the different months for the two taxa. We had planned to sample the leaves of *P. serratifolia* from spring to autumn. However, the marked branches in early March were pruned by gardeners at the beginning of September 2021. Fortunately, the leaves growing from late February to early March 2021 were fully mature in August 2021, and the leaves sampled from mid-March to mid-August manifest temporal changes in leaf shape

from young to mature leaves of this hybrid. When a leaf matures, its leaf shape seldom changes significantly.

### *2.2. Data Acquisition*

The fresh leaves were scanned by a photo scanner (V550, Epson, Batam, Indonesia) to obtain .tiff images at a 600-dpi resolution, which were transferred to black-white .bmp files using the Photoshop software (CS6, version: 13.0; Adobe, San Jose, CA, USA). The Matlab (version ≥ 2009a; MathWorks, Natick, MA, USA) procedures developed by refs. [25,26] were used to extract the planar coordinates of each leaf boundary. The 'biogeom' package (version 1.1.1) [27] based on the statistical software R (version 4.2.0) [28] was then used to calculate leaf area, length, and width using the planar coordinates.

### *2.3. Methods*

The Montgomery equation (ME) [5] was used to describe the relationship between leaf area (*A*) and the product of leaf length (*L*) and width (*W*):

$$A = m\_p \cdot L\mathcal{W}\_\prime \tag{1}$$

where *m<sup>p</sup>* is the proportionality coefficient, which is the Montgomery parameter. To normalize the data, we log-transformed the data of both sides of Equation (1):

$$y = a + \mathfrak{x},\tag{2}$$

where *y* = ln(*A*); *x* = ln(*LW*); *a* = ln(*mp*). We used ordinary least-squares protocols to estimate the parameter *a*, and its 95% confidence interval (95% CI). To test whether there is a significant difference in *m<sup>p</sup>* between the two taxa, we used the bootstrap percentile method [29,30] to calculate the 95% CI of the differences in the 4000 bootstrapping replicates of *m<sup>p</sup>* between *P.* × *fraseri* and *P. serratifolia*. If the 95% CI of the differences includes 0, there is no significant difference in *m<sup>p</sup>* between the two taxa; if it does not include 0, a significant difference is detected. We used the Pearson's product moment correlation coefficient (*r*) to measure the validity (i.e., statistical robustness) of the linear relationship between *x* and *y*:

$$\sigma = \frac{\sum\_{i=1}^{n} \left(\boldsymbol{x}\_{i} - \overline{\boldsymbol{x}}\right) \left(\boldsymbol{y}\_{i} - \overline{\boldsymbol{y}}\right)}{\sqrt{\sum\_{i=1}^{n} \left(\boldsymbol{x}\_{i} - \overline{\boldsymbol{x}}\right)^{2}} \sqrt{\sum\_{i=1}^{n} \left(\boldsymbol{y}\_{i} - \overline{\boldsymbol{y}}\right)^{2}}} \tag{3}$$

where *x* and *y* represent the means of *x*- and *y*-values, respectively; and *n* is the sample size. The test of the significance of the correlation coefficient is based upon a statistic that theoretically follows a *t* distribution with *n* − 2 degrees of freedom if the samples of *x* and *y* are normally distributed.

The leaf ellipticalness index (EI) was used to quantify leaf shape [23]:

$$\text{EI} = \frac{A}{(\pi/4)LW}.\tag{4}$$

The EI has a close relationship with *mp*, i.e., EI = *mp*/(π/4). A linear regression method was used to estimate the *m<sup>p</sup>* value based on the number of leaves. However, EI can also be calculated directly using *A*, *L* and *W*. The EI quantifies the extent of a leaf shape deviating from an ellipse regardless of the eccentricity of the ellipse. Thus, we used the *W/L* ratio as another leaf-shape index. ANOVA with Tukey's honestly significant difference (HSD) test [31] and a 0.05 significance level were used to test the significance of the differences in EI as well as the *W/L* ratio between the two taxa at the combined data level and at the intraspecific level across different dates of collection. For any two groups, we calculated the 95% confidence interval (CI) of the Tukey's test confidence interval, which equals the difference in the means between the two groups ± HSD, and we determined its significance by checking whether the 95% CI included 0. We also used the linear and parabolic equations

to fit the EI (and the *W/L* ratio) vs. collection month data, where month was set to a numeric variable to test whether leaf shape has a monotonic or parabolic trend across different investigation dates (i.e., leaf ages). to be fitted (referred to as regression coefficients). For the simple linear regression, <sup>2</sup> β was set to be zero. Each regression coefficient's statistic in Equation (5) follows a *t* distribution with *n* − *p* degrees of freedom, where *n* is the sample size and *p* is the number of

*Plants* **2022**, *11*, x FOR PEER REVIEW 4 of 9

different investigation dates (i.e., leaf ages).

$$\text{LS} = \beta\_0 + \beta\_1 \text{Month} + \beta\_2 \text{Month}^2,\tag{5}$$

where LS is the leaf-shape index of interest (EI or the *W/L* ratio); Month is the month of sampling (ranging from 3 to 8 [i.e., March to August 2021] for *P. × fraseri*, and from 4 to 11 [i.e., April to November 2021] for *P. serratifolia*); and <sup>0</sup> β , <sup>1</sup> β and <sup>2</sup> β are the parameters

as well as the *W/L* ratio between the two taxa at the combined data level and at the intraspecific level across different dates of collection. For any two groups, we calculated the 95% confidence interval (CI) of the Tukey's test confidence interval, which equals the difference in the means between the two groups ± HSD, and we determined its significance by checking whether the 95% CI included 0. We also used the linear and parabolic equations to fit the EI (and the *W/L* ratio) vs. collection month data, where month was set to a numeric variable to test whether leaf shape has a monotonic or parabolic trend across

> 2 01 2 LS =+ + β β Month β Month *,* (5)

where LS is the leaf-shape index of interest (EI or the *W/L* ratio); Month is the month of sampling (ranging from 3 to 8 [i.e., March to August 2021] for *P.* × *fraseri*, and from 4 to 11 [i.e., April to November 2021] for *P. serratifolia*); and β<sup>0</sup> , β<sup>1</sup> and β<sup>2</sup> are the parameters to be fitted (referred to as regression coefficients). For the simple linear regression, β<sup>2</sup> was set to be zero. Each regression coefficient's statistic in Equation (5) follows a *t* distribution with *n* − *p* degrees of freedom, where *n* is the sample size and *p* is the number of parameters [32]. The statistical software R (version 4.2.0) [28] was used to carry out calculations and to draw figures. The 'agricolae' package (version 1.3-5) was used to implement the Tukey's HSD test. **3. Results** 

The statistical software R (version 4.2.0) [28] was used to carry out calculations and to draw figures. The 'agricolae' package (version 1.3-5) was used to implement the Tukey's HSD test. The Montgomery equation (ME) was valid in calculating leaf area. For the pooled intraspecific data, the correlation coefficient was greater than 0.99 (*p* < 0.001). The esti-

### **3. Results** mated values of the Montgomery parameters of *P. × fraseri* and *P. serratifolia* were 0.6494

coefficient.

parameters [32].

The Montgomery equation (ME) was valid in calculating leaf area. For the pooled intraspecific data, the correlation coefficient was greater than 0.99 (*p* < 0.001). The estimated values of the Montgomery parameters of *P.* × *fraseri* and *P. serratifolia* were 0.6494 (95% CI = 0.6482, 0.6506) and 0.6981 (95% CI = 0.6971, 0.6990),respectively (Figure 2). There was a significant difference in leaf shape between the two taxa, i.e., the leaves of *P. serratifolia* were on average wider and more elliptical in shape (Figure 3). The means of the *W/L* ratios for the two taxa were 0.38 (*P.* × *fraseri*) and 0.40 (*P. serratifolia*); the means of the leaf elliptical indices (EIs) for the two taxa were 0.8276 (*P.* × *fraseri*) and 0.8894 (*P. serratifolia*). (95% CI = 0.6482, 0.6506) and 0.6981 (95% CI = 0.6971, 0.6990), respectively (Figure 2). There was a significant difference in leaf shape between the two taxa, i.e., the leaves of *P. serratifolia* were on average wider and more elliptical in shape (Figure 3). The means of the *W/L* ratios for the two taxa were 0.38 (*P. × fraseri*) and 0.40 (*P. serratifolia*); the means of the leaf elliptical indices (EIs) for the two taxa were 0.8276 (*P. × fraseri*) and 0.8894 (*P. serratifolia*).

**Figure 2.** Log-log bivariate plot of leaf area (*A*) vs. the product leaf length (*L*) and width (*W*). The symbol ෝ represents the estimated value of the Montgomery parameter, i.e., the estimated proportionality coefficient of the Montgomery equation. The subscripts Pf and Ps represent the hybrid **Figure 2.** Log-log bivariate plot of leaf area (*A*) vs. the product leaf length (*L*) and width (*W*). The symbol *m*ˆ *<sup>p</sup>* represents the estimated value of the Montgomery parameter, i.e., the estimated proportionality coefficient of the Montgomery equation. The subscripts Pf and Ps represent the hybrid *P.* × *fraseri*, and one of its parents *P. serratifolia*, respectively. \*\*\* denotes *p* < 0.001 for a correlation coefficient.

*P. × fraseri*, and one of its parents *P. serratifolia*, respectively. \*\*\* denotes *p* < 0.001 for a correlation

**Figure 3.** Comparison of leaf shape between the two *Photinia* taxa for the pooled data. (**A**) Ratio of leaf width to length; (**B**) Leaf ellipticalness index. In each panel, the letters on the whiskers show differences between the two taxa (i.e., taxa with different letters are significantly different at the α = 0.05 significance level using Tukey's HSD); the values at the top of whiskers represent the coefficients of variation (%) for each taxon; the horizontal solid line represents the median; and the red asterisk represents the mean. The whiskers extend to the most extreme data point, which is no more than 1.5 times the interquartile range from the box. The difference in mean leaf *W/L* ratios is equal to 0.0193 with HSD = 0.0017, and the corresponding 95% Tukey's test confidence interval is 0.0176, 0.0210, which does not include 0, indicating a significant difference; the difference in mean leaf ellipticalness indices is equal to 0.0617 with HSD = 0.0019, and the corresponding 95% Tukey's test **Figure 3.** Comparison of leaf shape between the two *Photinia* taxa for the pooled data. (**A**) Ratio of leaf width to length; (**B**) Leaf ellipticalness index. In each panel, the letters on the whiskers show differences between the two taxa (i.e., taxa with different letters are significantly different at the α = 0.05 significance level using Tukey's HSD); the values at the top of whiskers represent the coefficients of variation (%) for each taxon; the horizontal solid line represents the median; and the red asterisk represents the mean. The whiskers extend to the most extreme data point, which is no more than 1.5 times the interquartile range from the box. The difference in mean leaf *W/L* ratios is equal to 0.0193 with HSD = 0.0017, and the corresponding 95% Tukey's test confidence interval is 0.0176, 0.0210, which does not include 0, indicating a significant difference; the difference in mean leaf ellipticalness indices is equal to 0.0617 with HSD = 0.0019, and the corresponding 95% Tukey's test confidence interval is 0.0598, 0.0636, which also does not include 0, indicating a significant difference.

confidence interval is 0.0598, 0.0636, which also does not include 0, indicating a significant difference. There were significant differences in leaf shape (measured as either the *W/L* ratio or the EI) among the different months of collecting (Figure 4). Three of the four slopes of the linear model of leaf shape vs. month were statistically significant (with three *p* values < 0.05) with the exception of the EIs of *P. serratifolia*. However, the four coefficients of deter-There were significant differences in leaf shape (measured as either the *W/L* ratio or the EI) among the different months of collecting (Figure 4). Three of the four slopes of the linear model of leaf shape vs. month were statistically significant (with three *p* values < 0.05) with the exception of the EIs of *P. serratifolia*. However, the four coefficients of determination (i.e., *r* <sup>2</sup> values) were all less than 0.06, i.e., a very weak linear relationship was observed between leaf shape and sampling month. Thus, there is no robust evidence for a monotonic increase or decrease in leaf shape across leaf ages. The linear and quadratic coefficients in the parabolic model were statistically significant (*p* values < 0.05) for the leaf *W/L*

mination (i.e., *r*2 values) were all less than 0.06, i.e., a very weak linear relationship was observed between leaf shape and sampling month. Thus, there is no robust evidence for a monotonic increase or decrease in leaf shape across leaf ages. The linear and quadratic

leaf *W/L* ratio vs. month relationship and the EI vs. month relationship for *P. × fraseri*. However, the two coefficients were not statistically significant (*p* values > 0.05) for any leaf-shape indices for *P. serratifolia*. Three of the four coefficients of determination of the parabolic regression were smaller than 0.03, with the exception of the EI vs. month data of *P. × fraseri* (*r*2 = 0.1452). Thus, there was also no strong evidence for a parabolic trend in

leaf shape across the different leaf ages (see Table 1 for details).

*Photinia × fraseri*  Leaf width/length ratio

*Photinia × fraseri*  Leaf ellipticalness index

ratio vs. month relationship and the EI vs. month relationship for *P.* × *fraseri*. However, the two coefficients were not statistically significant (*p* values > 0.05) for any leaf-shape indices for *P. serratifolia*. Three of the four coefficients of determination of the parabolic regression were smaller than 0.03, with the exception of the EI vs. month data of *P.* × *fraseri* (*r* <sup>2</sup> = 0.1452). Thus, there was also no strong evidence for a parabolic trend in leaf shape across the different leaf ages (see Table 1 for details). *Plants* **2022**, *11*, x FOR PEER REVIEW 6 of 9

**Figure 4.** Comparison of leaf shapes (reflected by the ratio of leaf width to length and leaf ellipticalness index) across different sampling months for *P. × fraseri* (**A**,**B**) and *P. serratifolia* (**C**,**D**). In each panel, the letters on the whiskers represent the significance of the differences in leaf shape between any two sampling months among which the letters a, b, c, and d are used to represent the significance of differences in leaf shape among different sampling months; groups sharing a common letter are not significantly different at the 0.05 significance level; the values at the top of whiskers represent the coefficients of variation (%) for each sampling month; the horizontal solid line represents the median; the red asterisk represents the mean. The whiskers extend to the most extreme data point, which is no more than 1.5 times the interquartile range from the box. For the label of the *x*-axis, the months 3 to 11 represent March to November 2021, respectively. The differences in the means of leaf-shape indices, HSD values, and the corresponding 95% Tukey's test confidence intervals between any two groups for *P. × fraseri* and *P. serratifolia* are provided in Tables S3 and S4, respectively. **Figure 4.** Comparison of leaf shapes (reflected by the ratio of leaf width to length and leaf ellipticalness index) across different sampling months for *P.* × *fraseri* (**A**,**B**) and *P. serratifolia* (**C**,**D**). In each panel, the letters on the whiskers represent the significance of the differences in leaf shape between any two sampling months among which the letters a, b, c, and d are used to represent the significance of differences in leaf shape among different sampling months; groups sharing a common letter are not significantly different at the 0.05 significance level; the values at the top of whiskers represent the coefficients of variation (%) for each sampling month; the horizontal solid line represents the median; the red asterisk represents the mean. The whiskers extend to the most extreme data point, which is no more than 1.5 times the interquartile range from the box. For the label of the *x*-axis, the months 3 to 11 represent March to November 2021, respectively. The differences in the means of leaf-shape indices, HSD values, and the corresponding 95% Tukey's test confidence intervals between any two groups for *P.* × *fraseri* and *P. serratifolia* are provided in Tables S3 and S4, respectively.

*a* 0.3733 <0.001

*c* −0.0006 0.0031

*y = a + bx <sup>a</sup>*0.8549 <0.001 0.0591 *b* −0.0050 <0.001

*y = a + bx + cx*<sup>2</sup> *<sup>a</sup>*0.9668 <0.001 0.1452 *b* 0.0504 <0.001

*b* 0.0054 0.0365 0.0233

**Taxon/Leaf-Shape Index Equation † Parameter Estimate Significance** *r***<sup>2</sup>**

*y = a + bx + cx*<sup>2</sup>

**Table 1.** Regression statistics of the leaf-shape index vs. sampling month.


**Table 1.** Regression statistics of the leaf-shape index vs. sampling month.

† Here, *y* represents a specific leaf-shape index, and *x* represents the sampling month.

### **4. Discussion and Conclusions**

Shi et al. [33] used the Montgomery parameter (*mp*) to fit leaf *A* vs. *LW* of 101 bamboo taxa using > 10,000 bamboo leaves. In their study, the estimated value of *m<sup>p</sup>* was 0.6959 (95% CI = 0.6952, 0.6966). Li et al. [23] estimated a *m<sup>p</sup>* value of 0.6840 (95% CI = 0.6827, 0.9855) using > 2200 leaves from nine Magnoliaceae species, whereas Ma et al. [13] estimated a *m<sup>p</sup>* value of 0.7710 (95% CI = 0.7696, 0.7714) using > 6500 leaves of an alpine oak species (*Quercus pannosa*) from different tree size groups. The present study shows that the estimated value of *m<sup>p</sup>* of the two *Photinia* taxa also approximates 0.7, i.e., the leaf area approximately accounts for 70% of the area of a hypothetical rectangle with the leaf length and width as its two side lengths. This result is significantly different from the numerical value of the alpine oak provided in ref. [13], because the leaf shape of *Q. pannosa* is a special superellipse with a mean leaf ellipticalness index greater than unity [34–36], which is larger than those of the two *Photinia* taxa examined in this study. Since leaf shape is found to be closely related to climatic factors [18,37], it would be worth studying the link between the *m<sup>p</sup>* and climate for some widely distributed plants in the future.

The leaf shapes of two *Photinia* taxa were quantified using different metrics, and found to significantly differ in spite of the fact that the two types of leaves appear to be superficially very similar in appearance. The leaves of *P. serratifolia* are wider and more elliptical than those of its hybrid. In addition, leaf shape significantly varied across the different leaf ages, but manifested no monotonic or parabolic tendency across the different leaf ages. The variations in leaf shape among the different months of sampling may reflect variations resulting from random samplings of individual plants. This variation, however, does not mask the fact that the leaf shape of the *Photinia* taxa is relatively stable, and that the Montgomery equation, which assumes a proportional relationship between leaf area and the product of leaf length and width, is valid at different leaf growth stages. These results highlight an approach that permits the non-destructive estimation of leaf area, i.e., direct field measurements of leaf length and width are used to estimate lamina area by multiplying their product by a proportionality coefficient (*mp*).

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants11182370/s1, Table S1: The raw leaf data of *Photinia* × *fraseri*; Table S2: The raw leaf data of *P. serratifolia*; Table S3: Tukey multiple comparisons of the means of leafshape indices between any sampling months for *P.* × *fraseri*; Table S4: Tukey multiple comparisons of the means of leaf-shape indices between any sampling months for *P. serratifolia*.

**Author Contributions:** Methodology, K.J.N., D.A.R. and P.S.; formal analysis, K.J.N., D.A.R. and P.S.; investigation, X.Z., K.J.N., D.A.R., Y.J., H.D. and P.S.; writing—original draft preparation, X.Z. and P.S.; writing—review and editing, K.J.N. and D.A.R.; funding acquisition, X.Z. and H.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** X.Z. and H.D. were funded by the Project of Biological Resources Survey in Wuyishan National Park (grant number: HXQT2020120701).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in the present work have been listed in the online Supplementary Materials.

**Acknowledgments:** We thank Kexin Yu for her valuable help during the preparation of this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

