• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Strong controls of daily minimum temperature on the autumn photosynthetic phenology of subtropical vegetation in China

    2021-10-12 08:11:20PeixinRenZelinLiuXiaoluZhouChanghuiPengJingfengXiaoSonghanWangXingLiandPengLi
    Forest Ecosystems 2021年3期

    Peixin Ren ,Zelin Liu ,Xiaolu Zhou ,Changhui Peng,2 ,Jingfeng Xiao ,Songhan Wang ,Xing Li and Peng Li*

    Abstract Background:Vegetation phenology research has largely focused on temperate deciduous forests,thus limiting our understanding of the response of evergreen vegetation to climate change in tropical and subtropical regions.Results:Using satellite solar-induced chlorophyll fluorescence (SIF) and MODIS enhanced vegetation index (EVI)data,we applied two methods to evaluate temporal and spatial patterns of the end of the growing season (EGS) in subtropical vegetation in China,and analyze the dependence of EGS on preseason maximum and minimum temperatures as well as cumulative precipitation.Our results indicated that the averaged EGS derived from the SIF and EVI based on the two methods (dynamic threshold method and derivative method) was later than that derived from gross primary productivity (GPP) based on the eddy covariance technique,and the time-lag for EGSsif and EGSevi was approximately 2 weeks and 4 weeks,respectively.We found that EGS was positively correlated with preseason minimum temperature and cumulative precipitation (accounting for more than 73% and 62% of the study areas,respectively),but negatively correlated with preseason maximum temperature (accounting for more than 59%of the study areas).In addition,EGS was more sensitive to the changes in the preseason minimum temperature than to other climatic factors,and an increase in the preseason minimum temperature significantly delayed the EGS in evergreen forests,shrub and grassland.Conclusions:Our results indicated that the SIF outperformed traditional vegetation indices in capturing the autumn photosynthetic phenology of evergreen forest in the subtropical region of China.We found that minimum temperature plays a significant role in determining autumn photosynthetic phenology in the study region.These findings contribute to improving our understanding of the response of the EGS to climate change in subtropical vegetation of China,and provide a new perspective for accurately evaluating the role played by evergreen vegetation in the regional carbon budget.

    Keywords:Carbon cycle,Evergreen vegetation,Plant phenology,Solar-induced Fluorescence,Climate change,MODIS,Eddy covariance

    Introduction

    Vegetation phenology refers to the rhythm of growth and development in the life cycle of plants,which is closely associated with seasonal changes in the environment (Suepa et al.2016).The life cycle comprises biological events that occur throughout the year such as budding,flowering,fruiting,defoliation and dormancy during the growth of plants in a year (Helmut 1974;Vrieling et al.2018).Vegetation phenology is a key indicator of climate change,and it can have a significant impact on the cycles of carbon,water,and energy in terrestrial ecosystems (Xiao et al.2009;Richardson et al.2013).Apart from gaining a better understanding of plant responses to climate change,studies on vegetation phenology can help us better understand the mechanisms involved in the exchange of matter and energy between vegetation and the atmosphere,and more accurately evaluate the contribution of vegetation to the global carbon budget (Piao et al.2008;Penuelas et al.2009;Richardson et al.2013).

    A large number of studies based on ground observations and satellite remote sensing monitoring have reported an advance in the start of vegetation growing season (SOS)during spring (Menzel et al.2006;Fu et al.2015) and a delay in the end of growing season (EGS) in autumn (Liu et al.2016;Piao et al.2019).Evidence from multiple studies indicated earlier SOS was significantly related to the increase of preseason temperature (Ge et al.2015;Xu et al.2018).Compared with SOS,the EGS also determines the length of the entire growing season and plays an important role in maintaining the global carbon balance (Piao et al.2008;Garonna et al.2014).However,the mechanisms underlying the response of EGS to climate change remain unclear(Gallinat et al.2015;Wu et al.2018).Though multiple studies tried to investigate the relationships between EGS and climate factors,there is still no consistent conclusion(Yang et al.2017;Wu et al.2018).For example,recent studies revealed that increase in daytime maximum temperature and nighttime minimum temperature would cause contrasting effects on drought stress,which result in inconsistent relation between autumn vegetation phenology with preseason maximum temperature and minimum temperature(Wu et al.2018).

    Previous studies on vegetation phenology have focused on the deciduous forest in the middle and high latitudes(Fu et al.2014;Flynn and Wolkovich 2018),while fewer studies have been conducted in subtropical regions with evergreen vegetation.Remote sensing based indicators have been frequently applied in vegetation phenology monitoring,such as Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI)(Wen et al.2017;Yuan et al.2018).However,these vegetation indices based on“greenness”observations can only reflect the greenness information of vegetation,and represent the“potential photosynthesis”of vegetation(Liu et al.2018).Subtropical regions mainly consist of evergreen forests and do not show significant seasonal changes with respect to canopy“greenness”.Therefore,using vegetation index to understand vegetation phenology and the response mechanism of vegetation phenology to climate change in these forests can cause significant bias (Karkauskaite et al.2017).

    The emergence of solar-induced chlorophyll fluorescence (SIF) in the recent decade provides a new opportunity to monitor vegetation phenology from regional to global scales (Guanter et al.2014;Yoshida et al.2015;Sun et al.2017).SIF is measured as light released during photosynthesis between the wavelengths of 650–800 nm,which can directly reflect the dynamic changes of plant photosynthesis (Frankenberg et al.2014).Compared to traditional vegetation indices that are significantly affected by atmospheric aerosols,soil,snow,and clouds(Zhang et al.2003;Balzarolo et al.2016;Liu et al.2018),SIF can more directly reflect the dynamic changes that occur during photosynthesis and is less sensitive to the influence of clouds or the atmosphere (Joiner et al.2014).Therefore,SIF is considered to be a reliable remote sensing-based indicator for monitoring the phenology of tropical or subtropical evergreen vegetation,especially in forests that are in the early stages of environmental stress (K?hler et al.2018;Zuromski et al.2018).In addition,SIF was also considered more effective in retrieving vegetation phenology in high productivity areas (Guanter et al.2014;Yang et al.2015).

    A quarter of the land area in China is covered by subtropical forests,characterized by unique vegetation types,high biodiversity,and remarkable ecological functions.These forests play an important role in maintaining the ecological balance in the region.However,there were still few studies on the phenology of subtropical vegetation,especially in EGS.Therefore,understanding the response of subtropical vegetation to climate change is critical,especially in terms of autumn phenology.In this study,we examined the autumn photosynthetic phenology for five vegetation types in subtropical China,including evergreen coniferous forest (ECF),evergreen broadleaved forest (EBF),deciduous broadleaved forest(DBF),shrub and grassland.The SIF,NIRv(near-infrared reflectance of vegetation)and MODIS EVI remote sensing data from 2000 to 2018 were used to extract the autumn phenology of subtropical vegetation in China.The main objectives of the study were as follows:(1)to examine the spatiotemporal patterns of the EGS in the study region;(2) to compare the results of vegetation photosynthetic phenology based on SIF and vegetation index;(3) to analyze the responses of vegetation autumn phenology in subtropical China to maximum temperature (Tmax),minimum temperature(Tmin)and precipitation and then to explore the underlying mechanisms.

    Materials and methods

    Datasets

    Remote sensing data

    The satellite SIF data used in this study is provided by Xiao (2019),which is a global high spatio-temporal resolution (0.05°,8-day) SIF dataset (namely,GOSIF) based on the OCO-2 (Orbiting Carbon Observatory 2) satellite.GOSIF was developed by a data-driven method which established a predictive SIF model from discrete OCO-2 soundings,MODIS remote sensing data and meteorological reanalysis data (Li and Xiao 2019a).It has similar and reasonable seasonal period as the original OCO-2 SIF,but it has higher spatial and temporal resolution,global continuous coverage and longer data record (Li and Xiao 2019a).This dataset plays an important role in understanding the long-term trends in global photosynthesis,and which has been widely used to evaluate the inter-annual variation in ecosystem productivity (Li and Xiao 2019b;Li and Xiao 2020;Qiu et al.2020).In calculation,we changed the temporal resolution of SIF data from 8 to 16 days by using the maximum synthesis method.The EVI data used in this study were extracted from NASA Earth Science Data (NASA 2015).For the analysis,we used EVI data (2000–2018) extracted from the MOD13C1 v006 dataset,with a spatial resolution of 0.05° and 16 days’ interval.The NIRv was calculated by the product of normalized vegetation index (NDVI) and near infrared reflectance,with a spatial resolution of 0.05 degree (Wang et al.2021).

    Meteorological data

    Meteorological data were obtained from the National Tibetan Plateau Data Center (2019).The dataset is based on the fusion of remote sensing product,reanalysis data set and field station data and has a temporal resolution of 3 h and a spatial resolution of 0.1° (Yang et al.2010).The dataset provides seven near-surface meteorological elements,including air temperature,surface pressure,specific humidity,wind speed,downward shortwave radiation,downward long-wave radiation and precipitation rate (He et al.2020).In this study,we used daily precipitation,Tmaxand Tmindata to examine the response of EGS to climate change.

    Flux and vegetation data

    In this study,carbon flux data of half hour scale in terrestrial ecosystem were used to evaluate the performance of different remote sensing data for monitoring vegetation phenology.The flux data of the Dinghushan and Qianyanzhou research stations from 2003 to 2010 were obtained from the national flux network of China(China-FLUX 2013).Coordinate axis rotation and WPL (Webb-Pearman-Leuning) correction were used to eliminate the effects of topography,air hydrothermal transmission,and observation height on the observed data.Then,the CO2flux data were partitioned into gross primary productivity(GPP) and total ecosystem respiration.The vegetation data (spatial resolution,1 km) were obtained from the Joint Research Center of European Commission under the project of Global Land Cover 2000.The final regional vegetation classification data for China were obtained by preprocessing the corresponding data(Xu et al.2005).We excluded cultivated areas affected by anthropogenic activity.For the final analysis,we selected five vegetation types in the subtropical region of China,including evergreen coniferous forest (ECF),evergreen broadleaved forest(EBF),deciduous broadleaved forest (DBF),shrub and grassland.The distribution of the subtropical region and vegetation type were shown in Fig.1.

    Fig.1 The spatial distribution of vegetation types and flux stations in study area

    Estimation of EGS

    In order to eliminate the background noise,the Savitzky-Golay filter was applied to smooth the SIF,EVI,NIRv and GPP flux time series data (Zhang et al.2016).We then used the dynamic threshold method and the derivative method to determine EGS,which indicates autumn vegetation phenology (Liu et al.2016;Filippa et al.2016).Compared with the fixed threshold method,the dynamic threshold method is advantageous since it eliminates the influence of background noise by allowing a threshold to be set based on the conditions in the study area (White et al.2009).The equation of the dynamic threshold method is as follows:

    where Xtis the value of X at a given time t,and Xmaxand Xminare the maximum and minimum values of X in the annual X time series,respectively.X indicates SIF,EVI and NIRv remote sensing data.In this study,the EGS date was defined as the first day of the year in the descending period when the Xratioreduction value was less than 0.5 as the EGS (Wu et al.2018).

    The derivative method assumes that the EGS is the time point when the X time series decreases rapidly:this corresponds to the points with the minimum slope of the fitting curve (Forkel et al.2015).Here the Xratiois calculated as the change in X at time t,based on the following equation:

    where Xtrepresents the value of X corresponding to the time t.Δt is value of the time variation.X indicates SIF,EVI and NIRv remote sensing data.In order to make the retrieval phenology more accurate,EGS retrieved from two methods was both used for the analysis.

    Analysis

    We first extracted the autumn photosynthetic phenology(i.e.,EGS) for each year at each pixel using both derivative and dynamic threshold methods from the GOSIF,NIRv and MODIS EVI datasets for the period from 2000 to 2018.We then calculated the annual average of EGS and analyzed its spatial distribution in the subtropical region.Subsequently,we used flux data to retrieve the annual EGS by the two methods in the study area,in order to evaluate the EGS derived from different remote sensing data.A simple linear regression was used to analyze the spatial distribution and temporal trends of EGS for each pixel from 2000 to 2018.In this analysis,we also compared the trends of the EGS across different vegetation types.

    Previous studies showed that precipitation,minimum temperature and maximum temperature all play an important role in regulating vegetation phenology (Piao et al.2019;Wang et al.2019b).A partial correlation analysis was used to evaluate the response of the EGS to climatic factors,including precipitation,Tmax,and Tmin,during preseason 5 months (with 1 month step).We used the absolute values of the maximum partial correlation coefficients for each pixel to identify the preseason period that was significantly associated with EGS,which was designated as the optimal preseason periods in the study area.Then,we analyzed the relationship between climatic factors and the EGS during the most related preseason periods and determined its significance of each pixel.Using the correlation coefficients,we also assessed the relationship between EGS and climatic factors across different vegetation types.

    Finally,to further investigate the response of EGS to climate factors,we conducted multiple regressions to evaluate the sensitivity of EGS to the preseason Tmax,Tminand cumulative precipitation.The coefficients of each factor in the regression model indicate the sensitivity of EGS to corresponding climatic factors.In this study,we aggregated all data to a 0.1°×0.1° grid to match the coarsest resolution among all datasets.In all the calculation and analysis,we excluded the area with low vegetation coverage (EVI<0.08) and retained those areas covering the five types of vegetation in the study.

    Results

    Spatial and temporal patterns of EGS in subtropical China

    There was a distinct latitudinal variation of EGS derived from SIF and EVI data (Fig.2):an advance of the EGS at higher latitudes and a delay of the EGS at lower latitudes.Furthermore,EGS extracted from SIF data (EGSsif) was earlier than that derived from the EVI data(EGSevi).For the two datasets,the spatial distributions of EGS were similar between two methods.In the study area,EGSsifextracted by the derivative and dynamic threshold methods were 280.1 and 276 days,respectively,with an average of 278 days (Fig.2b and d).In contrast,EGSeviof the two methods were 294 and 289.6 days,respectively,with an average of 291.8 days (Fig.2a and c).In addition,some interesting information was found by comparing EGS estimated by the two types of remote sensing data and flux data from 2003 to 2010 in Fig.3.At the Qianyanzhou station,the EGS of EVI,SIF and GPP flux data (average of 8 years) by the derivative method were on day 303.9,291.6 and 278.9,respectively(Fig.3a).EGSevi,EGSsifand EGSGPPestimated by the dynamic threshold method were on day 304.1,289.6 and 279.5,respectively (Fig.3b).The average values of EGS retrieved from two kinds of remote sensing data were later than EGSGPP,and the time-lags were 25 days(EGSevi) and 11 days (EGSsif),respectively.At the Dinghushan station,EGSevi,EGSsifand EGSGPPby the derivative method were on day 300.4,288 and 273.9,respectively (Fig.3c).For the dynamic threshold method,EGSevi,EGSsifand EGSGPPwere on day 299.9,287.1 and 272.7,respectively (Fig.3d).The average of EGSGPPwas 27 days (EGSevi) and 14 days (EGSsif) ahead of the two remote sensing data.Compared with SIF and EVI,EGS retrieved by the NIRv was later for the two methods(Supplementary Material:Figs.S1 and S2).

    Fig.2 The spatial patterns of the end of growing season(EGS)in subtropical vegetation in China from 2000 to 2018:MODIS EVI dataset and SIF dataset.a-b,derivative method and c-d,dynamic threshold method.Inset plots (the bottom-left of each figure)display the frequency distribution of EGS

    Fig.3 Comparison of the end of growing season(EGS)retrieved from remote sensing data and flux data.Two ecological monitoring stations were listed,including Qianyanzhou station(QYZ)and Dinghushan station(DHS).Both methods were listed,including derivative and dynamic threshold methods

    Across the study area,the trends in the EGS extracted from the SIF and EVI by the two methods were similar(Fig.4).In terms of spatial distribution,the delays in EGSeviand EGSsifwere restricted to the central part of the study area (Fig.4).For the derivative method,the delay of EGSeviand EGSsifwas observed in more than 60% of the total study area,of which~20% showed significant delays (P<0.05) (Fig.4a and b).A delay of EGSeviand EGSsifextracted by the dynamic threshold method was observed in more than 55% of the total study area,and the pixels with significant delay (P<0.05) for EGSeviand EGSsifaccounted for 13.23%,and 17.77% of all the pixels in the study region,respectively(Fig.4c,d).Both EGSeviand EGSsifextracted by the two methods exhibited similar trends across different vegetation types.A delayed EGS was observed across most vegetation types (e.g.,evergreen forest,shrub,and grassland),except for the deciduous broadleaved forest (Fig.S3).In contrast,the averaged EGS of DBF retrieved from the two methods showed an advancing trend by 0.09 day·yr?1(EGSevi) or 0.37 day·yr?1(EGSsif).In addition,in the study area and different vegetation types,the trend in the EGS extracted from the NIRv was different from that of the SIF and EVI (Figs.S4,S5).

    Fig.4 The spatial patterns of the linear trend of the end of growing season(EGS)of subtropical vegetation in China from 2000 to 2018:MODIS EVI dataset and SIF dataset.a-b,derivative method and c-d,dynamic threshold method.A negative value indicates an advance,and a positive value indicates a delay.Inset plots(the bottom-left of each figure)display the frequency distribution of change trend.The proportions of positive(P)and negative(N)(proportions of significant in parentheses) trends are provided

    Response of EGS to climate drivers

    Based on the above analysis,the SIF showed a better performance than EVI and NIRv in capturing the EGS,and thus we chose EGSsifto explore the relationship between autumn phenology and climate factors in the study area.At the regional scale,for both two methods,EGSsifin subtropical China was correlated with Tminduring the period of 2–4 months prior to EGSsif,the median and mean of the period related to Tminwere in the 3 months prior to EGSsif.The mean of the period related to Tmaxwas in the 2 months to EGSsif.For the cumulative precipitation,the EGSsifwas most correlated with the period about 3 months prior to EGSsif(Fig.S6).

    Based on the partial correlation analysis,we found that a large proportion of the pixels showed positive correlations between EGSsifextracted by the two methods and Tminas well as cumulative precipitation (Fig.5).For the derivative and dynamic threshold methods,about 77.1%and 73.58% of the pixels covering the study area showed a positive partial correlation with Tmin,of which 17.61%and 16.63% showed statistically significant relationships(P<0.05),respectively.Similarly,on average,62.91% of the pixels covering the study area of the two methods showed a positive partial correlation with preseason cumulative precipitation,and this correlation was significant over 9.5% of pixels (P<0.05).Compared with the preseason Tminand cumulative precipitation,there were more pixels with negative correlation between the preseason Tmaxand EGSsif,accounting for 59.96% and 59.36% of the study area for the two methods,respectively (Fig.5).

    Fig.5 Spatial pattern and frequency distribution of partial correlation coefficient between the end of growing season(EGSsif)and climatic factors:a,d precipitation(Pre),b,e maximum temperature(Tmax)and c,f minimum temperature(Tmin).The methods of the left panel and the right panel are derivative and dynamic threshold methods respectively.The percentages of positive(P)and negative(N)correlations(P<0.05,percentage of significant correlations in parentheses)are provided

    The responses of EGS of different vegetation types to climate factors were different (Fig.6).For two methods,there was a positive partial correlation between EGSsifand Tminacross different vegetation types (RP>0.51,P<0.05;Fig.6),with the exception of deciduous broadleaved forest (RP0.05).Preseason cumulative precipitation was positively correlated with EGSsiffor all vegetation types,and the strongest relationship was observed in the shrub (RP>0.56,P<0.05).In contrast,Tmaxwas negatively correlated with EGSsifacross all vegetation types,except for DBF.Overall,EGSsifwas more strongly correlated with preseason Tminand cumulative precipitation than with preseason Tmax(Fig.6).To further test the results of the partial correlation analysis,we also analyzed the sensitivity of EGSsifto climate factors.Both methods showed similar results.The sensitivity of EGS to the preseason Tminwas the strongest across all vegetation types (>1.93 day·sd?1),except for DBF (Fig.S7).Conversely,the sensitivity of EGS to the preseason Tmaxwas the weakest across different vegetation types (Fig.S7).

    Fig.6 Correlation between the end of growing season(EGSsif)of different vegetation and climatic factors.Five main vegetation types in this study area were listed,such as evergreen coniferous forest(ECF),evergreen broadleaved forest(EBF),deciduous broadleaved forest(DBF),shrub,and grassland.a derivative method and b dynamic threshold method.*indicates statistically significant trends at the 95%(P<0.05)level

    Discussion

    Comparison of satellite-retrieved EGS based on SIF and EVI data

    By comparing SIF with GPP estimated across two flux tower sites in the study area,the SIF dataset used in this study exhibited strong seasonal and interannual dynamics that were consistent with those of daily GPP (Fig.S8),and thus the reconstructed SIF product has a great potential for monitoring the photosynthetic phonology in the study area.We also found that the SIF clearly had an advanced EGS than did EVI,and it was closer to GPP-derived EGS at the two evergreen forest sites.The differences in the seasonal cycle of EVI and SIF could be explained by the differences in the information contained in the two kinds of data resource.On the one hand,the SIF is deemed to be directly linked to photosynthetic activity,which contains major information on photosynthetically active radiation (Walther et al.2016).While the EVI is more of an approximation for fraction of photosynthetically active radiation,which indicates the photosynthetic potential of the terrestrial vegetation cover (Jiang et al.2008).Because of the intrinsic limitations of the photosynthetic machinery and external stress factors,the absorbed photosynthetically active radiation by vegetation cannot be completely used for carbon fixation (Baker 2008).In contrast,SIF contains information on not only absorbed photosynthetically active radiation but also environmental stresses that determine photosynthetic light use efficiency (Yoshida et al.2015;Li and Xiao 2020).Therefore,SIF can be used to track changes in physiological changes induced by environmental stresses in the absence of changes in greenness or structure.This may explain why EGS based on SIF occurred earlier than that based on the“greenness”phenology reflected in EVI data.

    On the other hand,SIF is an energy flux emitted from plant chlorophyll molecules a few nanoseconds after light absorption by vegetation (Baker 2008),whereas the EVI was calculated based on vegetation canopy reflectance.Multiple studies have reported that physiological vegetation indices (e.g.,chlorophyll) performance better than structural vegetation indices for detecting the autumn photosynthetic phenology for evergreen forests (Wong et al.2019;Yin et al.2020).The structural recession of evergreen forest in autumn is gradual,and the photosynthetic rate is mainly controlled by physiology (Gallinat et al.2015).Our result indicated that autumn photosynthesis in subtropical China forests is mainly stressed by minimum temperature variability,which causes photosynthesis to end before structural recession (Jeong et al.2017).Therefore,SIF outperformed EVI in detecting the photosynthetic phenology.In addition,although EVI is an advanced vegetation index,it is still affected by clouds and other atmospheric noise (Miura et al.2001;Huete et al.2002).In contrast,the SIF is not sensitive to the influence of cloud and atmosphere (Joiner et al.2014).This may also contribute to the difference in phenology monitoring for SIF and EVI in our study.

    Response of vegetation phenology to climate drivers

    Different from the earlier spring phenology which was mainly caused by global warming,there is no consistent conclusion on the change trend of EGS and its influencing factors.Therefore,it is difficult to explain the response mechanism of EGS to climate change.Especially for the photosynthesis phenology in the subtropical China,the climate response mechanism of autumn photosynthesis is largely unknown due to the limitation of monitoring methods.Our results indicated that the EGS in subtropical China was slightly delayed,and the change trend of EGS was biome dependent.In temperate regions,multiple studies have also confirmed that the EGS has reported a small delay in recent years;however,there were no widespread delaying trends in autumn phenology (Wang et al.2019a).Our findings indicated that global climate change can extend the growing season in subtropical vegetation,which can in turn enhance the carbon sink capacity of subtropical vegetation.

    Previous studies have shown that temperature plays a key role in regulating vegetation autumn phenology(Cleland et al.2007;Chuine et al.2010).For example,the increase of Tmaxcan advance or delay autumn phenology,while the increase of Tminhad the opposite effect(Wu et al.2018).Our results indicated that the increase of preseason Tminsignificantly delayed autumn photosynthetic phenology,extending the growing season in subtropical vegetation in China.The increase of the preseason Tminmay delay the coloring of leaves in autumn and reduce the chilling injury caused by low nighttime temperature(Yang et al.2017).Additionally,the increase of Tmincan result in warmer autumn weather,fewer frost days and a delay in the first frost (Liu et al.2016).In our study,the sensitivity of each vegetation type to the preseason Tminwas stronger.This could be due to the fact that subtropical regions are warm and humid,and the nighttime temperature is a limiting factor for vegetation growth,thereby leading to a greater impact on vegetation growth and development than daytime temperature and precipitation.On the other hand,the increase of the preseason Tmaxled to an earlier EGS and shortened the growing season in our study.The higher subtropical daytime temperature can decrease photosynthetic enzyme activity (Rossi et al.2017),which can in turn inhibit photosynthesis and shorten the growing season.In addition,the increase of daytime temperature can also lead to higher evapotranspiration and lower soil water use efficiency,resulting in earlier senescence of vegetation (Estiarte and Penuelas 2015;Wu et al.2018).

    Our results also indicated that preseason cumulative precipitation had a positive impact on the EGS,resulting in a longer growing season.Water is one of the important components of protoplasm,and the amount of water in vegetation affects its metabolic intensity and photosynthetic rate (Quetin and Swann 2017).The sensitivity of vegetation to water varied with stage of growth (Quetin and Swann 2017).In the subtropical region,the increase of cumulative precipitation will strengthen the absorption of nutrients by vegetation and promote the effective photosynthesis (Bertani et al.2017),delaying the EGS.Additionally,the impact of precipitation on vegetation phenology in the subtropical region has regional and biological specificity (Zhou et al.2018).Although there is abundant precipitation in the subtropical region of China,the evergreen forests are more suitable in wet areas.Therefore,the increase of preseason precipitation can enhance the physiological activity of vegetation and photosynthesis,thus delaying the vegetation autumn phenology.

    It should be noted that the preseason Tmaxand Tminhad opposite effects on EGS in deciduous broadleaved forests in the study area.Wu et al.(2018)also reported that the Tmaxand Tminmay have the opposite effects on vegetation autumn phenology in different regions.The elevated Tmincan increase nighttime respiration and consumption of organic compounds,leading to a shorter growing season;in contrast,higher Tmaxpromotes photosynthesis and delays the EGS (Wu et al.2018).We speculated that the increase of the maximum temperature delayed the EGS,which may be related to the higher latitude position of deciduous broadleaved forest.In this region,the temperature is relatively low,and the increase of the maximum temperature is beneficial for photosynthesis.Our results also indicated that the effect of minimum temperature on EGS was stronger than that of maximum temperature,which may provide an explanation for the advancing trend of EGS in deciduous broadleaved forest.Another possible explanation is that there had been a slight decrease in rainfall across the deciduous broadleaved forest region over the past two decades (Table S1),which may inhibit the extension of the vegetation growing season.Therefore,the comprehensive effect of climate factors on EGS could partly explain the advancing trend in EGS for deciduous broadleaved forest.

    Implications

    Using remote sensing data to retrieve the phenological information of deciduous forest has been well reported(Yang et al.2015;Liu et al.2016),while the research on phenological monitoring of evergreen forest was still scarce.Our findings indicated that the remote sensing phenological monitoring based on SIF was closer to the photosynthetic phenology in subtropical vegetation.Our results confirmed a lag between the autumn decrease of photosynthesis and the change in greenness in evergreen forests (Walther et al.2016).This implies that estimates of the EGS purely based on greenness indices will be biased in evergreen forest in subtropical regions,which translates into errors in the autumn carbon budget.Therefore,compared with traditional vegetation indices,SIF can better capture the autumn decrease stage of photosynthesis of subtropical vegetation and effectively improve the dynamic monitoring of photosynthetic activity in evergreen ecosystems of subtropical regions.Our results provide a new reference for the study of subtropical vegetation phenology,and demonstrate the potential of SIF for simulating carbon budget in evergreen ecosystems of subtropical regions.

    In the context of global climate change,it needs to be explained,that how vegetation phenology responds to climate change.However,the response mechanisms of autumn phenology of subtropical vegetation to climate change remain unclear.Our results indicated that the effects of precipitation,maximum and minimum temperature on autumn phenology of subtropical vegetation were discrepant,and elucidating this inconsistency is beneficial to the establishment of subtropical vegetation phenology models.Furthermore,our results suggested the most significant effect of the minimum temperature on autumn phenology of subtropical vegetation,which can improve the understanding of the control factors of subtropical vegetation phenology.Based on our study,we proposed to further explore the effects of climate change on autumn phenology of evergreen forests in other regions of the world.This is an important implication for the improvement of phenological parameterization of terrestrial ecosystem models.

    Conclusions

    In this study,we used SIF and EVI data to examine the spatial and temporal variation of autumn vegetation phenology and to analyze its response to climatic factors in subtropical vegetation in China.We found that the delay of EGS occurred in evergreen forests,shrub and grassland,but not in the deciduous broadleaved forest.We found that the preseason Tminand cumulative precipitation were positively associated with the delay of EGS (with positive correlation for more than 73% and 62% of the study areas for the two methods,respectively).Conversely,the preseason Tmaxwas negatively associated with the EGS (with negative correlation for more than 59% of the study areas for both methods).In all vegetation types except deciduous broadleaved forest,the increase of Tminalso caused the delay of EGS (>1.93 day·sd?1),while the increase of preseason Tmaxadvanced it.For precipitation,the increase of cumulative precipitation could delay the EGS across all vegetation types.Our study indicated that the preseason Tminhad a significant effect on the photosynthetic phenology of subtropical evergreen vegetation,providing new insights into how climate change affects the EGS.These results also provide a scientific basis for the development of phenology models for evergreen vegetation.

    Abbreviations

    SIF:Solar-induced Chlorophyll Fluorescence;EVI:Enhanced Vegetation Index;NIRv:Near-infrared Reflectance of Vegetation;EGS:The End of the Growing Season;SOS:The Start of the Growing Season;ECF:Evergreen Coniferous Forest;EBF:Evergreen Broadleaved Forest;DBF:Deciduous Broadleaved Forest;Tmax:Maximum Temperature;Tmin:Minimum Temperature;GPP:Gross Primary Productivity;NDVI:Normalized Difference Vegetation Index;WPL:Webb-Pearman-Leuning

    Supplementary Information

    The online version contains supplementary material available at https://doi.org/10.1186/s40663-021-00309-9.

    Additional file 1:Figure S1.The spatial patterns of the end of growing season (EGS) in subtropical vegetation in China from 2000 to 2018:NIRv dataset.a) derivative method and b) dynamic threshold method.Inset plots(the bottom-left of each figure) display the frequency distribution of EGS.Figure S2.Comparison of the end of growing season (EGS) retrieved from NIRv and flux data.Two ecological monitoring stations were listed,including Qianyanzhou station (QYZ) and Dinghushan station(DHS).Both methods were listed,including derivative and dynamic threshold methods.Figure S3.Linear trends of the end of growing season(EGS) across China’s Subtropical biomes from 2000 to 2018.Five main vegetation types in study area were listed,including evergreen coniferous forest(ECF),evergreen broadleaved forest(EBF),deciduous broadleaved forest(DBF),shrub,and grassland.a) derivative and b) dynamic threshold method.A negative value indicates an advance,and a positive value indicates a delay.* indicates statistically significant trends at the 90% (P<0.1)level,**indicates statistically significant trends at the 95% (P<0.05)level.Figure S4.The spatial patterns of the linear trend of the end of growing season (EGS) of subtropical vegetation in China from 2000 to 2018:NIRv dataset.a) derivative method and b) dynamic threshold method.Inset plots(the bottom-left of each figure) display the frequency distribution of change trend.The proportions of positive (P)and negative (N) (proportions of significant in parentheses) trends are provided.Figure S5.Linear trends of the end of growing season (EGS) across China’s Subtropical biomes from 2000 to 2018:NIRv dataset.Five main vegetation types in study area were listed,including evergreen coniferous forest (ECF),evergreen broadleaved forest (EBF),deciduous broadleaved forest (DBF),shrub,and grassland.a) derivative and b) dynamic threshold method.A negative value indicates an advance,and a positive value indicates a delay.* indicates statistically significant trends at the 90% (P<0.1) level,**indicates statistically significant trends at the 95% (P<0.05) level.Figure S6.Optimal preseason periods depicting correlations between the end of growing season (EGS) derived from SIF data and climatic factors:Precipitation (Pre),maximum temperature (Tmax)and minimum temperature (Tmin).Two methods (derivative and dynamic threshold methods) were listed.Figure S7.Sensitivity of end of growing season(EGS) to climatic factors in different vegetation types.Five main vegetation types in study area were listed,including evergreen coniferous forest(ECF),evergreen broadleaved forest(EBF),deciduous broadleaved forest(DBF),shrub,and grassland.a) derivative and b) dynamic threshold method.1 day·sd?1denoted that an increase of 1 standard deviation (sd)in the climatic factors delayed or advanced the EGS by 1 day.Figure S8.The seasonal cycles of SIF and flux tower GPP from 2003 to 2010.Table S1.The change rate and significance of three climate factors in different vegetation areas from 2000 to 2018,including precipitation (Pre),maximum temperature (Tmax)and minimum temperature (Tmin).Five main biomes in this study area were listed,such as evergreen coniferous forest(ECF),evergreen broadleaved forest(EBF),deciduous broadleaved forest(DBF),shrub,and grassland.

    Acknowledgements

    The authors would like to acknowledge the National Tibetan Plateau Data Center,National Aeronautics and Space Administration (NASA) and ChinaFLUX for providing the data.We thank the two anonymous reviewers for helping us improve the manuscript.

    Authors’ contributions

    PR acquired the data,analyzed the data,and drafted the manuscript.PL analyzed the data and conceived and designed the study.PL,ZL,XZ,CP,JX,SW and XL revised the important intellectual content critically.The authors read and approved the final manuscript.

    Funding

    This study was financially supported by the National Natural Science Foundation of China (Grant No.41901117),Natural Science Foundation of Hunan Province,China (Grant No.2020JJ5362),the Outstanding Youth Project of Hu’nan Provincial Education Department (No.18B001),and the Natural Sciences and Engineering Research Council of Canada (NSERC)Discover Grant.

    Availability of data and materials

    The GOSIF product is available from https://globalecology.unh.edu/data.html.The MODIS EVI product analysed during the current study is available in the NASA,https://ladsweb.modaps.eosdis.nasa.gov/search/.The meteorological data is available from https://data.tpdc.ac.cn/zh-hans/.The flux data is available in the national flux network of China,http://www.cern.ac.cn/0index/index.asp.Other datasets used during the current study are available from the authors on reasonable request.

    Declarations

    Ethics approval and consent to participate

    Not applicable.

    Consent for publication

    Not applicable.

    Competing interests

    The authors declare that they have no competing interests.

    Author details

    1College of Geographic Science,Hunan Normal University,Changsha 410081,China.2Department of Biology Sciences,Institute of Environment Sciences,University of Quebec at Montreal,C.P.8888,Succ.Centre-Ville,Montreal H3C 3P8,Canada.3Earth Systems Research Center,Institute for the Study of Earth,Oceans,and Space,University of New Hampshire,Durham,NH 03824,USA.4International Institute for Earth System Sciences,Nanjing University,Nanjing 210023,China.5Jiangsu Provincial Key Laboratory of Geographic Information Technology,Key Laboratory for Land Satellite Remote Sensing Applications of Ministry of Natural Resources,School of Geography and Ocean Science,Nanjing University,Nanjing 210023,China.

    Received:4 January 2021Accepted:27 April 2021

    e午夜精品久久久久久久| 99国产精品99久久久久| 变态另类成人亚洲欧美熟女| 天天躁狠狠躁夜夜躁狠狠躁| 变态另类丝袜制服| 亚洲五月婷婷丁香| 美女扒开内裤让男人捅视频| 欧美一区二区国产精品久久精品 | 老司机午夜十八禁免费视频| 脱女人内裤的视频| 一二三四在线观看免费中文在| 日韩成人在线观看一区二区三区| 舔av片在线| 最近在线观看免费完整版| 色尼玛亚洲综合影院| 欧美性长视频在线观看| 国产在线观看jvid| 亚洲精品中文字幕一二三四区| 三级毛片av免费| 欧美最黄视频在线播放免费| 日韩中文字幕欧美一区二区| 亚洲成人国产一区在线观看| 久久亚洲真实| 国产高清videossex| 又黄又爽又免费观看的视频| 亚洲片人在线观看| 一个人观看的视频www高清免费观看 | 男女视频在线观看网站免费 | 欧美午夜高清在线| 欧美黑人巨大hd| 中出人妻视频一区二区| 老熟妇乱子伦视频在线观看| 日本 av在线| 两个人视频免费观看高清| 国产69精品久久久久777片 | 亚洲精品一区av在线观看| 国产三级黄色录像| 国产麻豆成人av免费视频| 一级毛片女人18水好多| 每晚都被弄得嗷嗷叫到高潮| 亚洲精品国产精品久久久不卡| 两性夫妻黄色片| 精品久久久久久成人av| 777久久人妻少妇嫩草av网站| 色在线成人网| 黄频高清免费视频| 亚洲欧美精品综合久久99| 国产亚洲精品一区二区www| 国产精品影院久久| 极品教师在线免费播放| 美女 人体艺术 gogo| 啦啦啦韩国在线观看视频| 999精品在线视频| 亚洲18禁久久av| 嫁个100分男人电影在线观看| 丰满人妻熟妇乱又伦精品不卡| 国产精品一区二区三区四区免费观看 | 又黄又粗又硬又大视频| 国产一区在线观看成人免费| 色播亚洲综合网| 特大巨黑吊av在线直播| 激情在线观看视频在线高清| 草草在线视频免费看| 最近最新中文字幕大全免费视频| 亚洲第一欧美日韩一区二区三区| 亚洲午夜精品一区,二区,三区| 五月伊人婷婷丁香| 免费搜索国产男女视频| 国产精品久久电影中文字幕| 国产精品亚洲一级av第二区| www.999成人在线观看| 嫁个100分男人电影在线观看| 久久人妻av系列| 全区人妻精品视频| 91麻豆精品激情在线观看国产| 久久天躁狠狠躁夜夜2o2o| 欧美中文综合在线视频| 一级黄色大片毛片| 88av欧美| 国产成人aa在线观看| 真人做人爱边吃奶动态| 亚洲真实伦在线观看| 丝袜美腿诱惑在线| 人妻久久中文字幕网| 亚洲一区中文字幕在线| 一进一出好大好爽视频| 好男人在线观看高清免费视频| 脱女人内裤的视频| 精品人妻1区二区| 极品教师在线免费播放| 精品久久久久久成人av| 亚洲精品一卡2卡三卡4卡5卡| 草草在线视频免费看| 久久久久国内视频| 啦啦啦韩国在线观看视频| 日韩免费av在线播放| 国产黄a三级三级三级人| 精品熟女少妇八av免费久了| 国产亚洲精品av在线| 欧美av亚洲av综合av国产av| av超薄肉色丝袜交足视频| 国产av一区在线观看免费| 国产精品免费视频内射| 亚洲精品美女久久久久99蜜臀| 给我免费播放毛片高清在线观看| 一级毛片精品| 欧美日本亚洲视频在线播放| 一个人观看的视频www高清免费观看 | 国产v大片淫在线免费观看| 亚洲国产精品合色在线| 久久中文看片网| 禁无遮挡网站| 日韩欧美在线乱码| 在线观看一区二区三区| 给我免费播放毛片高清在线观看| netflix在线观看网站| 免费看a级黄色片| 午夜亚洲福利在线播放| 亚洲第一电影网av| 我的老师免费观看完整版| 九色成人免费人妻av| 三级毛片av免费| 国产精品亚洲美女久久久| 99re在线观看精品视频| 欧美黑人巨大hd| 两个人看的免费小视频| 欧美+亚洲+日韩+国产| 久久久精品欧美日韩精品| 欧美3d第一页| 老司机在亚洲福利影院| 亚洲av成人av| 国产97色在线日韩免费| 久久精品亚洲精品国产色婷小说| 国产av一区二区精品久久| 午夜免费激情av| www.www免费av| 国产又黄又爽又无遮挡在线| 天天躁夜夜躁狠狠躁躁| 99国产综合亚洲精品| 精品午夜福利视频在线观看一区| 亚洲成人免费电影在线观看| a级毛片a级免费在线| 中出人妻视频一区二区| 日本黄大片高清| 女生性感内裤真人,穿戴方法视频| av福利片在线| 熟妇人妻久久中文字幕3abv| 宅男免费午夜| 精品无人区乱码1区二区| 99国产精品一区二区三区| 国产av一区二区精品久久| 久久午夜亚洲精品久久| 国产熟女xx| av国产免费在线观看| 丝袜人妻中文字幕| netflix在线观看网站| www日本黄色视频网| 中国美女看黄片| 99精品在免费线老司机午夜| 五月伊人婷婷丁香| 欧美人与性动交α欧美精品济南到| 精品高清国产在线一区| 黄色女人牲交| 99精品久久久久人妻精品| 搡老熟女国产l中国老女人| 男女下面进入的视频免费午夜| 亚洲五月婷婷丁香| 搡老熟女国产l中国老女人| 一区福利在线观看| 亚洲人与动物交配视频| 国产黄片美女视频| av有码第一页| 国产精品久久电影中文字幕| 999久久久国产精品视频| 亚洲九九香蕉| 男女那种视频在线观看| 亚洲乱码一区二区免费版| 人成视频在线观看免费观看| 天堂av国产一区二区熟女人妻 | 日韩有码中文字幕| 国产精品一区二区三区四区免费观看 | 亚洲成人免费电影在线观看| 日韩欧美免费精品| 精品第一国产精品| 丁香欧美五月| 国产真实乱freesex| 国产亚洲精品av在线| 日本 av在线| 1024手机看黄色片| 亚洲av五月六月丁香网| 色综合站精品国产| 法律面前人人平等表现在哪些方面| 人妻夜夜爽99麻豆av| 亚洲九九香蕉| 人妻丰满熟妇av一区二区三区| 精品久久久久久,| 麻豆av在线久日| 精品国内亚洲2022精品成人| 国产不卡一卡二| 欧美乱妇无乱码| 好看av亚洲va欧美ⅴa在| 国产主播在线观看一区二区| 精品无人区乱码1区二区| 婷婷精品国产亚洲av| 女生性感内裤真人,穿戴方法视频| 成人三级做爰电影| 欧美一区二区精品小视频在线| 一夜夜www| 成人18禁高潮啪啪吃奶动态图| 99久久久亚洲精品蜜臀av| 男女之事视频高清在线观看| 国产在线精品亚洲第一网站| 丁香欧美五月| 热99re8久久精品国产| 国语自产精品视频在线第100页| 亚洲精品美女久久av网站| 老司机午夜福利在线观看视频| 国产午夜福利久久久久久| 久久中文字幕一级| 日韩大尺度精品在线看网址| 亚洲国产精品999在线| 一本久久中文字幕| 久久久久国产精品人妻aⅴ院| 亚洲精品美女久久av网站| 精品午夜福利视频在线观看一区| 好男人电影高清在线观看| 成人午夜高清在线视频| 人人妻人人看人人澡| 两个人免费观看高清视频| 国产成人aa在线观看| 国产黄片美女视频| 久久久国产精品麻豆| 国产成人啪精品午夜网站| 一区二区三区高清视频在线| 精品久久久久久久末码| 一本综合久久免费| 久久精品综合一区二区三区| 亚洲人成网站高清观看| 午夜a级毛片| 一进一出抽搐动态| 妹子高潮喷水视频| av超薄肉色丝袜交足视频| 在线国产一区二区在线| 亚洲av成人av| 亚洲av电影不卡..在线观看| 国产av在哪里看| 禁无遮挡网站| 久热爱精品视频在线9| 精品久久久久久久毛片微露脸| 中文字幕最新亚洲高清| 国产又色又爽无遮挡免费看| 97超级碰碰碰精品色视频在线观看| 亚洲 欧美 日韩 在线 免费| 人成视频在线观看免费观看| 又黄又粗又硬又大视频| 久久久久久大精品| 亚洲 欧美一区二区三区| 男女下面进入的视频免费午夜| 狂野欧美白嫩少妇大欣赏| 亚洲专区国产一区二区| 久久精品国产亚洲av高清一级| 国产1区2区3区精品| 国产欧美日韩一区二区三| e午夜精品久久久久久久| 妹子高潮喷水视频| www.www免费av| 不卡av一区二区三区| 国产一区二区激情短视频| 亚洲国产精品久久男人天堂| 色在线成人网| 免费在线观看亚洲国产| 国内少妇人妻偷人精品xxx网站 | 亚洲精品av麻豆狂野| 在线永久观看黄色视频| 成年女人毛片免费观看观看9| 少妇粗大呻吟视频| 18禁裸乳无遮挡免费网站照片| 又黄又粗又硬又大视频| av在线天堂中文字幕| 夜夜躁狠狠躁天天躁| 午夜a级毛片| e午夜精品久久久久久久| 成人av一区二区三区在线看| 久久九九热精品免费| 亚洲国产精品合色在线| 听说在线观看完整版免费高清| 夜夜躁狠狠躁天天躁| 少妇被粗大的猛进出69影院| 成熟少妇高潮喷水视频| 亚洲美女黄片视频| 欧美日韩福利视频一区二区| 妹子高潮喷水视频| 在线观看免费日韩欧美大片| 中文亚洲av片在线观看爽| 嫩草影视91久久| 青草久久国产| 亚洲午夜精品一区,二区,三区| 亚洲狠狠婷婷综合久久图片| 亚洲av成人一区二区三| 天堂动漫精品| 老司机在亚洲福利影院| aaaaa片日本免费| 精品不卡国产一区二区三区| 国产黄色小视频在线观看| 国产精品野战在线观看| 天堂√8在线中文| 在线播放国产精品三级| 成人三级黄色视频| 国产精品亚洲av一区麻豆| 又粗又爽又猛毛片免费看| 男插女下体视频免费在线播放| 老熟妇乱子伦视频在线观看| 中文字幕人妻丝袜一区二区| av片东京热男人的天堂| 人妻丰满熟妇av一区二区三区| 久久香蕉精品热| 欧美日韩乱码在线| 亚洲av电影不卡..在线观看| 身体一侧抽搐| 久久久水蜜桃国产精品网| 欧美日韩一级在线毛片| 男男h啪啪无遮挡| 18禁黄网站禁片免费观看直播| 免费人成视频x8x8入口观看| 人妻久久中文字幕网| 国内精品一区二区在线观看| 久久精品国产亚洲av香蕉五月| 亚洲五月婷婷丁香| 99国产精品99久久久久| 丁香欧美五月| 日本一二三区视频观看| 亚洲午夜理论影院| 亚洲国产看品久久| 国产精华一区二区三区| 亚洲午夜理论影院| 国产精品av视频在线免费观看| 嫁个100分男人电影在线观看| 国产激情偷乱视频一区二区| 日韩成人在线观看一区二区三区| 日韩大码丰满熟妇| 岛国在线观看网站| 久久久久久久精品吃奶| 欧美激情久久久久久爽电影| 国产成人精品无人区| 欧美 亚洲 国产 日韩一| 亚洲av电影在线进入| 成人欧美大片| 国产欧美日韩精品亚洲av| 两个人看的免费小视频| 国产精品,欧美在线| 欧洲精品卡2卡3卡4卡5卡区| 久久99热这里只有精品18| 亚洲avbb在线观看| 久久精品夜夜夜夜夜久久蜜豆 | 亚洲精品一卡2卡三卡4卡5卡| 亚洲国产精品成人综合色| 又黄又爽又免费观看的视频| 欧美黑人巨大hd| 免费无遮挡裸体视频| 亚洲一码二码三码区别大吗| 亚洲欧美精品综合一区二区三区| 这个男人来自地球电影免费观看| 午夜两性在线视频| 久久午夜亚洲精品久久| 成人av一区二区三区在线看| 在线观看免费午夜福利视频| 久久欧美精品欧美久久欧美| 首页视频小说图片口味搜索| www日本黄色视频网| 国产v大片淫在线免费观看| 成人三级黄色视频| 欧美日韩亚洲国产一区二区在线观看| 久久99热这里只有精品18| 日本三级黄在线观看| av在线播放免费不卡| 中出人妻视频一区二区| 一a级毛片在线观看| 欧美日韩亚洲国产一区二区在线观看| 久久这里只有精品中国| 国产精品 欧美亚洲| avwww免费| 成熟少妇高潮喷水视频| 99久久精品热视频| 久久人妻av系列| 岛国视频午夜一区免费看| 精品第一国产精品| 亚洲精品美女久久av网站| 啦啦啦韩国在线观看视频| 国产黄a三级三级三级人| 日日夜夜操网爽| 亚洲熟妇熟女久久| 亚洲中文日韩欧美视频| 少妇裸体淫交视频免费看高清 | 久久香蕉国产精品| 欧美高清成人免费视频www| 舔av片在线| 亚洲专区字幕在线| 亚洲av熟女| 亚洲午夜精品一区,二区,三区| 色播亚洲综合网| 国产精品98久久久久久宅男小说| 俺也久久电影网| 欧美日韩黄片免| 亚洲色图 男人天堂 中文字幕| 三级国产精品欧美在线观看 | 国产1区2区3区精品| 琪琪午夜伦伦电影理论片6080| 看黄色毛片网站| 亚洲人成77777在线视频| 欧美乱码精品一区二区三区| 在线国产一区二区在线| 国产精品爽爽va在线观看网站| 久久久久久久久中文| 90打野战视频偷拍视频| 9191精品国产免费久久| 叶爱在线成人免费视频播放| 曰老女人黄片| 国产亚洲av嫩草精品影院| 亚洲一区二区三区色噜噜| 欧美日韩福利视频一区二区| 午夜免费成人在线视频| 99久久综合精品五月天人人| 日韩欧美一区二区三区在线观看| 久久久久久久久免费视频了| 在线观看一区二区三区| 国产又黄又爽又无遮挡在线| 亚洲精品粉嫩美女一区| 香蕉久久夜色| 亚洲一码二码三码区别大吗| 又紧又爽又黄一区二区| tocl精华| 国产三级中文精品| 久久久国产成人精品二区| 国产精品美女特级片免费视频播放器 | 国产精品国产高清国产av| 亚洲五月天丁香| 日韩中文字幕欧美一区二区| 国产精品美女特级片免费视频播放器 | 日本免费一区二区三区高清不卡| 国产成人av教育| 日本五十路高清| 国产精品亚洲av一区麻豆| 午夜福利欧美成人| 久久热在线av| 久久久久久大精品| 成人高潮视频无遮挡免费网站| av超薄肉色丝袜交足视频| 夜夜看夜夜爽夜夜摸| 亚洲欧美日韩高清专用| av中文乱码字幕在线| 成熟少妇高潮喷水视频| 精品乱码久久久久久99久播| 国产激情欧美一区二区| 免费观看人在逋| 亚洲一区中文字幕在线| 亚洲五月婷婷丁香| 97碰自拍视频| 欧美另类亚洲清纯唯美| 88av欧美| 久久精品国产综合久久久| 最近最新中文字幕大全免费视频| 全区人妻精品视频| x7x7x7水蜜桃| 亚洲精品久久成人aⅴ小说| 久久国产精品人妻蜜桃| 国产亚洲欧美在线一区二区| 国产伦一二天堂av在线观看| 最近最新免费中文字幕在线| 亚洲成人久久性| 桃色一区二区三区在线观看| 男女下面进入的视频免费午夜| 观看免费一级毛片| 久久久国产精品麻豆| 日本熟妇午夜| 毛片女人毛片| 亚洲精品一区av在线观看| 亚洲精品国产一区二区精华液| 亚洲,欧美精品.| 99国产精品99久久久久| cao死你这个sao货| 色哟哟哟哟哟哟| 很黄的视频免费| av在线播放免费不卡| 亚洲人成网站高清观看| 黄色视频,在线免费观看| 日韩精品中文字幕看吧| 制服诱惑二区| av天堂在线播放| 亚洲avbb在线观看| 长腿黑丝高跟| 亚洲自偷自拍图片 自拍| 国模一区二区三区四区视频 | 亚洲精品美女久久av网站| 好看av亚洲va欧美ⅴa在| av超薄肉色丝袜交足视频| 人妻夜夜爽99麻豆av| 国产区一区二久久| 久久精品国产亚洲av香蕉五月| 人妻丰满熟妇av一区二区三区| 国内毛片毛片毛片毛片毛片| 久久香蕉精品热| 美女 人体艺术 gogo| avwww免费| 亚洲欧美日韩无卡精品| 国产精品永久免费网站| 亚洲av日韩精品久久久久久密| 亚洲国产精品合色在线| 最近在线观看免费完整版| 在线视频色国产色| 99热这里只有是精品50| 免费观看人在逋| 成人特级黄色片久久久久久久| 三级毛片av免费| 日本 av在线| 国产av麻豆久久久久久久| 日日干狠狠操夜夜爽| 在线观看www视频免费| 成年免费大片在线观看| 精品乱码久久久久久99久播| 亚洲一区二区三区不卡视频| 黄片大片在线免费观看| 亚洲成人精品中文字幕电影| 国产精品一区二区免费欧美| 99国产综合亚洲精品| 国产视频一区二区在线看| 免费在线观看亚洲国产| 1024香蕉在线观看| 国产亚洲精品一区二区www| a在线观看视频网站| 亚洲精品在线观看二区| 一本精品99久久精品77| 叶爱在线成人免费视频播放| 一区二区三区激情视频| 国产单亲对白刺激| 国产欧美日韩精品亚洲av| 欧美精品啪啪一区二区三区| 一级作爱视频免费观看| 亚洲精品久久成人aⅴ小说| 国产1区2区3区精品| 99久久无色码亚洲精品果冻| 亚洲全国av大片| 久久伊人香网站| 90打野战视频偷拍视频| 在线观看免费视频日本深夜| 国产区一区二久久| 日本精品一区二区三区蜜桃| 亚洲欧美日韩东京热| 亚洲一卡2卡3卡4卡5卡精品中文| av在线播放免费不卡| 国产av麻豆久久久久久久| 久久亚洲真实| 神马国产精品三级电影在线观看 | 久久久久久久午夜电影| 亚洲色图av天堂| 久久久精品国产亚洲av高清涩受| 欧美午夜高清在线| 国产伦一二天堂av在线观看| 后天国语完整版免费观看| 日韩 欧美 亚洲 中文字幕| 日日摸夜夜添夜夜添小说| a级毛片在线看网站| 久久久久亚洲av毛片大全| 久久久久精品国产欧美久久久| av在线天堂中文字幕| 国产午夜福利久久久久久| 精品午夜福利视频在线观看一区| 黄色视频不卡| 在线a可以看的网站| 免费在线观看成人毛片| 1024香蕉在线观看| 日韩国内少妇激情av| 成人av在线播放网站| 国内揄拍国产精品人妻在线| 亚洲av片天天在线观看| 禁无遮挡网站| 午夜福利欧美成人| 色综合婷婷激情| 亚洲五月婷婷丁香| 国产真人三级小视频在线观看| 国产精品野战在线观看| 在线观看午夜福利视频| 人人妻人人澡欧美一区二区| 高潮久久久久久久久久久不卡| 欧美在线一区亚洲| 神马国产精品三级电影在线观看 | 99国产综合亚洲精品| 国产一区二区三区视频了| 精品日产1卡2卡| 国产精品九九99| 国产视频一区二区在线看| 亚洲第一欧美日韩一区二区三区| 亚洲国产欧美人成| 欧美成人午夜精品| 亚洲欧美精品综合久久99| 视频区欧美日本亚洲| 女同久久另类99精品国产91| 成年女人毛片免费观看观看9| 久久香蕉国产精品| 一边摸一边做爽爽视频免费| 国产v大片淫在线免费观看| 视频区欧美日本亚洲| 亚洲 欧美 日韩 在线 免费| 人妻久久中文字幕网| 90打野战视频偷拍视频| 国产99白浆流出| 亚洲国产精品久久男人天堂| av在线天堂中文字幕| 亚洲中文字幕一区二区三区有码在线看 | www国产在线视频色| 免费在线观看影片大全网站| 可以免费在线观看a视频的电影网站| 国产一区二区激情短视频| 欧美日韩国产亚洲二区| 俺也久久电影网| 在线观看一区二区三区|