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

    A combination of climate,tree diversity and local human disturbance determine the stability of dry Afromontane forests

    2021-07-24 07:09:48HadguHisheLouisOosterlynckKidaneGidayWandaDeKeersmaeckerBenSomersandBartMuys
    Forest Ecosystems 2021年2期

    Hadgu Hishe,Louis Oosterlynck,Kidane Giday,Wanda De Keersmaecker,3,Ben Somers and Bart Muys

    Abstract

    Keywords:Climate,Dryland,Disturbance,Restoration,Tigray,Growth stability,Biodiversity function

    Introduction

    A significant area of the globe(41%)is covered with drylands,and a large part of the human population(35%)resides in them(Safriel and Adeel 2008).Among dryland ecosystems,the dry forest biome covers an estimated 1079 million ha(Bastin et al.2017),accounting for almost half of the(sub)tropical forests(Aide et al.2013).Dryland forests are very important for biodiversity conservation,as they are known for their high level of endemism (Myers et al.2000);also for deep aquifer recharge as they show high infiltration rates in a lacking water environment(Bargués-Tobella et al.2020),and for moderating high temperatures.Dryland forests are some of the most threatened by human degradation and therefore,maintaining the remnant forests is crucial for a sustainable environment,and a seed source for possible restoration(Safriel et al.2005;Díaz et al.2018).

    Dry forests are among the most threatened ecosystems(Bognounou et al.2010)as they are found in regions of low productivity,supporting population with one of the fastest birth rates,where poverty prevails(Safriel and Adeel 2008).Dry forests have high conversion rates to other land use,and the remaining parts are degraded and fragmented(Sánchez-Azofeifa et al.2005).

    Due to climate change and other anthropogenic causes,desertification is widespread in drylands and is impacting the overall well-being of dwellers(Yan et al.2011).Climate change-induced prolonged dryness could change the vegetation composition of dryland forests,which might further complicate the socioeconomic situation in these areas(Huang et al.2016).Local disturbance factors such as illegal logging,uncontrolled browsing and grazing,and fire incidences are adding on to,and are possibly interfering with,the effect of global climate change on dryland forests(Lloret et al.2007;Jacob et al.2014;Abrha and Adhana 2019;Hishe et al.2020).Understanding how forests respond to increasing climate change and local human pressure is crucial to keep a sustained flow of the ecosystem services,ecosystem stability(Jactel et al.2006;Bauhus et al.2017;Duffy et al.2017)and should be an essential component of forest management(Huang et al.2016).This is important as not all forests respond in the same way to global and local disturbances.Their responses are modulated by local landscape characteristics such as species composition,altitude,slope and edaphic factors.Diverse versus monoculture stands,for example,are reported to respond differently to disturbance(Johnson et al.1996;Van Ruijven and Berendse 2007;De Keersmaecker et al.2018).While a number of studies reported that tree diversity has a positive effect on production,health and stability of forests,other have reported either neutral or negative effect of diversity which indicates for a need of further study(Waide et al.1999;McCann 2000).As a consequence restoration planning protocols will need context specific information.

    Different metrics have been proposed to define and quantify the responses of forests to disturbances(Webb 2007;Yan et al.2011).Among these,growth stability,resilience and resistance have been used widely(Verbesselt et al.2016;De Keersmaecker et al.2018).Many definitions are given to the mentioned stability concepts(Nikinmaa et al.2020).The resilience is defined as the recovery rate after a disturbance(Dakos et al.2012).Resistance,on the other hand,is the capacity of the forest to remain unchanged regardless of disturbances(Grimm and Wissel 1997).Growth stability is considered as a steady continuity of growth irrespective of external disturbance(Chen et al.2019).

    Ecosystem stability is affected by different factors,such as climate,topography and species diversity,among others(Yan et al.2011;Hutchison et al.2018).Insight in the response of the ecosystem to change in these factors is valuable for management and restoration purposes.In the absence of long-term ecological experiments,remote sensing data analysis is providing an opportunity to monitor long term forest dynamics(Wang et al.2004).Typically,vegetation indices based on the ratio between the reflectance in red and near-infrared(NIR)bands,such as the Normalized Different Vegetation Index(NDVI)(Kogan 1995),are used to characterize vegetation properties(Lu et al.2016).NDVI time series thus provide valuable information on forest dynamics and their response to external pressures(Lhermitte et al.2011;Verbesselt et al.2016;De Keersmaecker et al.2018).

    Forest stability metrics can be derived by applying statistical analysis to the entire NDVI time series,holistic approach,to take the possible recurrent stochastic perturbation events such as drought and other environmental variations in an open environment into consideration(Verbesselt et al.2016;Hutchison et al.2018).Within the holistic approach,temporal autocorrelation(TAC)(Verbesselt et al.2016),the depth of the anomalies(De Keersmaecker et al.2014)and the standard deviation of the anomalies(Pimm 1984)from a decomposed time series are commonly used as an indicator of forest resilience and resistance,respectively.TAC is based on the assumption that forests with lower resilience will recover more slowly,and growth progress is dependent on previous performances(Verbesselt et al.2016).Hence,higher TAC values indicate a slow forest response to these perturbations,showing lower recovery rate of the system.TAC is thus a measure of the slowness of forest response after disturbances and a direct indicator of resilience(Verbesselt et al.2016).TAC can be used to assess how close a system is close to a critical transition point(CTP) to another stable system,the higher the autocorrelation(close to one)the closer the system is to the CTP(Leemput et al.2018).Subtracting the TAC from one,on the other hand,indicates how close a system is to its prior disturbance state(the recovery rate in its broad sense),which could be considered as the resilience of the system(Verbesselt et al.2016).

    Similarly,as resistance is defined as the ability to withstand external shocks where highly resistant forests will deviate less than forests with low resistance during perturbations,the depth of the deviation is considered as an indicator of resistance(De Keersmaecker et al.2014).In addition,growth stability can be measured by calculating the area under the curve of the undecomposed NDVI at a yearly basis and is measured by the inverse of the coefficient of variation(mean divided by the standard deviation)of the respective years of the time series(Isbell et al.2009).

    Apart from quantifying the degree of stability of forests to disturbances,understanding and predicting the effect of environmental factors strengthening or weakening forest stability is little explored(Yan et al.2011).Therefore,this research aims at quantifying the effect of different explanatory variables describing tree species diversity,local degradation indicators and climate on forest resilience,resistance and growth stability over time using MODIS NDVI time series.Such information will be crucial for planning a successful restoration and forest management(Anjos and De Toledo 2018).With this respect,the study strives to test the following hypotheses:1)precipitation and temperature play a vital role in the stability of dry forests,2)topographic and edaphic factors and local land degradation indicators further modulate the difference in the stability of forests,3)stands with multispecies composition have more growth stability resistance and resilience under climate fluctuation and human disturbances than monocultures.

    Methods

    Study area description

    The study was carried out in Desa’a Forest,a large degraded dry Afromontane forest situated in the Tigray and Afar regions in the north of Ethiopia,for which an ambitious restoration plan is ongoing.The altitudes range from 900 m in Afar lowlands to 3000 m in the highlands of Tigray(Fig.1).Due to the large difference in topography and long north-south extension along the escarpment,the geologic formation of the forest area is diverse(Asrat 2002).The bedrock in Desa’a Forest is mainly made up of a Precambrian basement in the northern part and the Hintalo limestone dotted with Adigrat Sandstone in the southern landscape(Williams 2016).

    The precipitation pattern of the study area is influenced by topography and rain-bearing winds and is dominated by a large inter-annual variability(Nyssen et al.2005).Data from a nearby meteo-station and Worldclim (http://worldclim.org/version2)(Fick and Hijmans 2017) indicate that the average annual temperature and precipitation of the study area ranges between 13°C to 25°C and 400 to 700 mm respectively.Drought has a long history in the area,and caused regular famines,including in recent times.Recent droughts have been recorded for 2000,2002,2004 and 2009(Gebrehiwot and van der Veen 2013).In a recent study,2012 and 2013 were added among the driest years in the region(Tefera et al.2019).

    Fig.1 Location of Desa’a Forest in Ethiopia,and the position of the sampling points

    Desa’a Forest is most often classified as a dry Afromontane forest with a long dry season,where Juniperus procera Hochst.ex Endl.and Olea europaea ssuubbsspp..cuspidata(Wall.ex G.Don)Cif.are the dominant species(Friis et al.2010)in the canopy and understory,respectively.In Aynekulu et al.(2012),dry Afromontane forest(Juniper-Olea-Tarchonanthus group),semi-deciduous shrubland(Cadia-Acacia group),open acacia woodland and semidesert shrubland(Balanites group)was identified from top to bottom along the altitude gradient.The forest is under strong degradation pressure by livestock and overcutting and is undergoing fast species composition change(Aynekulu et al.2011))with a 500 m upward shift in the tree line for juniper and olive species so far(Aynekulu et al.2011).Desa’a forest covers an area of 150,000 ha.

    Data collection

    Environmental factors

    The ground data were collected by systematic sampling,based on a 2 km by 2 km grid.At each corner of the grids,303 plots of 400 m2were established on which all woody species,shrubs and trees,were identified following the nomenclature of Ethiopian flora(Tesemma 2007)and counted.For each tree,diameter at breast height(DBH)at 1.3 m above ground was measured using a calliper.For shrubs,diameter at stump height(DSH)at 30 cm above ground was measured.Trees with at least 5 cm in DBH and shrubs with at least 1 cm in DSH were considered.Only plots with a vegetation cover above 10%following the FAO definition of forest,131 plots were used(FAO 2010).For the shrub and tree layers,canopy cover was estimated by a group of three experts and an average was recorded.

    For each plot,slope,aspect and altitude were extracted from the 30 m spatial resolution ASTER Digital Elevation Model.The 19 standard Bioclimatic variables for 30 years(1970–2000)were extracted at 1 km resolution from the WorldClim WebPortal (http://worldclim.org/version2)(Fick and Hijmans 2017).The definition and nature of the bioclimatic variables are well documented in Fick and Hijmans(2017)and(O’Donnell and Ignizio 2012).

    Distance to nearby settlements and roads were extracted from a Euclidean distance raster constructed from a digitized road and settlement shapefiles.The shapefiles were obtained from a combination of data digitized from Google earth,and GPS tracked major and feeder roads,towns and centre of encompassing villages.

    In every plot,local disturbance indicators such as fire incidence,grazing and logging severity were estimated(see appendix 1,supplementary material)following Aynekulu et al.(2011).In each of the diversity inventory plots,soil depth was measured by penetrating a metal rod until the bedrock is reached.The thickness of the forest floor(ectorganic humus layer)was measured after cutting a profile with a spade(Eriksson and Holmgren 1996)(Table 1).

    Satellite imagery

    Moderate Resolution Imaging Spectroradiometer(MODIS)satellite data,i.e.the global MOD13Q1 data product with a temporal resolution of 16 days and a spatial resolution of 250 m,was used.MODIS NDVI time series from 2001 to 2018 were downloaded from Google Earth Engine(Hird et al.2017).Upon downloading,low data quality observations such as pixels covered by clouds were masked(Hird et al.2017).NDVI values were extracted for the pixels covering each inventory plot for every scene as a matrix of bimonthly NDVI over the 18 years in R-software.

    Data analysis

    Time series decomposition

    The time series were decomposed into trend,seasonality and remainder(anomalies)components using Seasonal-Trend decomposition using Loess(STL)(Abbes et al.2018)in R software.The trend component indicates long-term forest development,while the seasonal component depicts annual growth variations(Quan et al.2016).The remainder is the difference obtained when the trend and seasonality are subtracted from the original time series(Verbesselt et al.2016)(Fig.2).

    Deriving ecosystem stability metrics from the NDVI time seriesThree stability metrics were used to describe forest dynamics:resilience,resistance,and growth stability.While resilience and resistance were based on the anomalies of the NDVI time series(De Keersmaecker et al.2014),growth stability was based on the integrals of the undecomposed NDVI time series(Isbell et al.2009).

    Table 1 Categorical environmental factors collected in the field(Lower rank indicates better forest condition and higher values indicate bad forest condition;while soil depth,humus depth and erosion status were assessed into five ranks,grazing,cutting and fire incidence were ranked into four)

    Fig.2 An example of an NDVItime series of Desa’a forest,study area,decomposed using the STL algorithm

    ResilienceResilience(Fig.3)was computed using the temporal auto-correlation(TAC)of the anomaly.TAC and resilience are given in the following formula(Dakos et al.2012),Eqs.1 and 2,respectively.Highly correlated events(=high TAC)represent a slow recovery rate(=low resilience).

    where TAC is the temporal autocorrelation at lag 1,Xtstands for the observation at time t and n equals the total number of observations.

    ResistanceThe resistance was calculated as the lowest 5th percentile of the remainder(anomalies)per year(De Keersmaecker et al.2014)(Fig.3).Small values for the resistance metric represent highly resistant forests,i.e.forests that will deviate to a small extent during perturbations.

    Growth stabilityThe growth stability was calculated from the integral of the undecomposed NDVI time series(Yin et al.2012).The area under the curve of yearly based NDVI time series was considered as a good proxy for the net primary production(growth)of the forest.This area under the curve was obtained based on the top 75%of the yearly NDVI response to avoid the possible effect of seasonal variation in vegetation properties such as leaf sheds(Fig.4).The growth stability was then calculated as the inverse of the coefficient of variation(i.e.a ratio of mean to standard deviation)of the area under the curve.

    Fig.3 The concept of resilience(A)resistance(B)as used in this study on the remainder of the time series decomposition.Resilience is the recovery rate of the community,the resistance if the net change in the community

    Tree diversity

    Basal area(BA)based on species diversity was derived using the Shannon-Wiener diversity index(H′)and evenness index(J)equations(Shannon 1948),Eqs.3 and 4,respectively.

    where H′is the Shannon-Wiener diversity index,J is Shannon-Wiener evenness index,and BAiis the BA proportion(n/N)of individuals of the abundance of the ithspecies(one particular species)found(n)divided by the total number of individuals found(N)(species richness),and S is the number of species.These diversity indices were later used as explanatory variables in the regression analysis.

    Statistical analyses

    The four forest stability metrics were modeled against climate,tree species diversity,edaphic and topographic variables and land degradation indicators.Boosted Regression Trees(BRT)was applied as a regression model(Elith et al.2008)for each metric to explain the dynamics of the forest as a system and identify the most important factors predicting each metric.

    BRT allows handling of complex interactions while allowing simplicity for ecological interpretation(Elith et al.2008;Aertsen et al.2012).BRT combines the power of regression trees and boosting.It continuously partitions the data into homogeneous parts and fits a specific model to each partition.This avoids the loss of unexplained data if a single regression model could be fitted into such complex interactions.In R-environment,BRT was run using the gbm.step function developed by Elith et al.(2008)which as an extension of the“gbm”package(Ridgeway 2007),and explanatory variables could be simplified to concentrate on the most meaningful and important ones using the gbm.simplify to boost the power of the model(Elith et al.2008).

    The different variables used in the analyses were checked for multi-collinearity using the variation inflation factor(VIF)and Pearson correlation.Variables with higher VIF(>5)and Pearson correlation(>0.7)between predictors were not included in the reported outputs(Aertsen et al.2012).BRT was run for the different stability metrics by varying the learning rates(0.001–0.05),tree complexity(1–5)and bag fraction (0.50–0.75).Model performance was measured using R-squared,AIC and root mean square error(RMSE).In the BRT,the cross-validation(CV)statistic is the most important measure to evaluate the results(Elith et al.2008).The cross-validation correlation is the mean correlation of the predicted data iteratively based on the number of folds(Elith et al.2008).The higher the correlation,the higher the predictive power of the model.Because the algorithm is of a stochastic nature,based on the bag fraction used(the default is 75%),a portion of the data(here 50%was used)is used to train the model and the remaining for prediction capability test.Variable importance is determined by averaging the number of times,a variable is selected in the iterative division(splitting)of data weighted by the squared improvement to the BRT model(Gu et al.2019).Variables that are above the median of the group in the model value are highly important(significant),and those that are below are less important variables in the model(Gu et al.2019).Results were also supported by partial dependence plots to ease ecological interpretability of the effect trend of the factors considered.

    Fig.4 Fraction of the yearly NDVI(75%)used to extract growth stability for Desa’a forest

    To generate wall to a wall map of stability metrics over the forest,a kriging interpolation in ArcMap10.6 was applied to the stability metrics obtained on a plot level.Similarly,the stability matrices were summarized on an annual basis to show the stability status of the forest over the study period.A summary of the methodological approach is presented in the flow chart below(Fig.5).

    Results

    Stability status of Desa’a forest and correlation of the metrics

    The resilience,resistance,and growth stability of Desa’a forest from 2001 to 2018 depict a similar trend(Figs.5 and 6,Table 2).The resilience index showed lows in the years 2001,2007 and 2015(Fig.6).The resistance showed minima in 2004,2008,2009 and 2015.The growth stability,however,was declining throughout the study period except for a sudden rise in 2016(Fig.7).Additionally,the spatial distribution of the four metrics showed similar patterns(Fig.8),where vegetation in the south was more stable while in the center of the study area it was less stable.In the north,however,it was more stable except for the resilience metric.

    The correlation between the stability metrics used shows that resilience(r=0.56)and resistance(r=0.46)correlated significantly with growth stability.However,the correlation between resistance and resilience was weak(0.23).The correlation among resilience,resistance and productivity was positive.

    Drivers of stability

    Drivers of resilience

    Resilience was influenced by a combination of biophysical and climatic factors.In general,precipitation of the wettest month,species evenness,distance from the settlement and slope were the most effective variables explaining the resilience of Desa’a forest.The other factors had a similar share of influence(Table 2).

    The partial dependencies of the variables in the model indicated that three main types of responses could be observed.First,the precipitation of the wettest month,annual precipitation,annual temperature,Shannon diversity,distance to settlement,and annual temperature range showed a similar trend.Their influence was increasing up to a certain optimal condition and ceiled afterwards.In all except the precipitation of the wettest month,visible reductions in resilience were observed before an ultimate increment was recorded.Second,the effect of both species evenness and slope showed a unimodal shape,high at the mid values and lower at the two ends.Third,temperature seasonality and stoniness showed a negative effect on the resilience of the forest(Fig.9).

    Drivers of resistance

    Temperature seasonality and temperature of the driest quarter,forest floor thickness and precipitation of the wettest month were the variables that influenced the resistance of the forest most,with a total contribution of 53.6%(Table 3).

    Fig.5 Methodological flowchart showing the workflow used in this paper.DIis Diversity indices,VIF is variance inflation factor,TAC is temporal autocorrelation,SD is the standard deviation,and BRT is boosted regression tree.Same colours show similar work stages;blue represents final input variables and green data sources,for example

    Fig.6 The NDVIderived resilience and resistance of Desa’a Forest between 2001 to 2018.The solid line is the average of each metrics of all plots in a particular year and the broken line is the linear trendline of each metric

    The partial dependency plots revealed that the important variables affecting resistance had two general effect trends.First,the influence of temperature seasonality ended up in a decreasing trend though they showed different responses in the process.The resistance of the forest was lower in areas where temperature seasonality was lower than 180(1.8°C),the optimal size of temperature seasonality and got pick at around 220(2.2°C)above which an increase in temperature seasonality resulted in reduced resistance of forest communities.Second,the effect of the mean temperature of the driest quarter,humus depth and precipitation of the wettest month followed a positive trend.Around 185 mm precipitation of the wettest month is optimal to keep a resistant forest in the dry Afromontane environment(Fig.10).

    Drivers of growth stability

    Growth stability was governed dominantly by precipitation of the wettest month,taking about 44%of the total effect.Annual temperature range,precipitation of the warmest quarter and distance to settlement had similar effect strength accounting for 56%of the total(Table 4).

    The partial dependencies of the factors influencing growth stability(Fig.11)show that the stability of the forest has been increasing with all the important factors.However,the increment rate was different across the factors.The growth stability remained low up to around 155 mm of precipitation of the wettest month,and it exponential increased and ultimately ceiled at 180 mm(Fig.11).

    Model strength of the different stability metrics

    The performance of the model fit to the different stability metrics is given in Table 5.Modelling growth stability with the variables used was difficult compared to the other response variables,resulting in the lowest performance for all goodness-of-fit criteria used(Table 5).

    Discussion

    Resilience,resistance and growth stability status of Desa’a forest

    Over the study period,Desa’a Forest remained more or less resistant but not resilient,with a significant decrease in resilience in 2001,2007 and 2015.A slight drop below the average resistance was also observed in 2004,2008,2009,and 2015.The frequent,and acute drought occurrences might explain these drops in both resilience and resistance in the region.In the study period,reported droughts occurred in 2000,2002,2004(Gebrehiwot and van der Veen 2013),2012 and 2013(Tefera et al.2019),and 2015(Ahmed et al.2017).The resilience range of Desa’a forest(0.3–0.6)is incomparably lower than that of other African tropical forests(0.7–1.0)reported by Verbesselt et al.(2016)which might explain the severe and repetitive anthropogenic pressure the forest is facing(Aynekulu et al.2011).The growth stability,however,was continuously decreasing over the study period,which might be linked to continuous degradation in the forest that could be explained by the dieback of the dominant species,olive and juniper trees(Aynekulu et al.2011),browsing and lopping of various species(Giday et al.2018).The frequent drought occurrences that were linked to the declined resilience of the forest might also be a reasonable explanation for the decreased yield stability.A clear increment of growth was,however,observed in 2016.This might be attributed to the increased rainfall recorded in 2016(Berhane et al.2020).Because there was an acute drought in 2015(Ahmed et al.2017)and significant increment in precipitation in 2016,growth might have positively affected the biomass production in the forest.

    Table 2 The relative influence of the variables determining resilience in Desa’a forest(in bold are significant factors)

    Fig.7 Growth stability in Desa’a Forest,2001 to 2018.The solid line is the average growth stability of all plots in a particular year and the broken line is the linear trend of the growth stability

    Fig.8 Spatial distribution of resilience,resistance,and growth stability in Desa’a Forest

    Fig.9 Partial dependencies of factors affecting resilience in Desa’a forest.The relative importance of variables in the model(%out of 100)is given in brackets.Fitted functions are centred around the mean of the resilience and plotted on a common scale.Rug plots(ticks in X-axis)show the distribution of sample measurements.PWem stands for precipitation of the wettest month,DiSet for distance from the settlement,AP for annual precipitation,ShanIfor Shannon index,TS for temperature seasonality and TAR for temperature annual range

    Among the determinants of resilience,resistance,and growth stability,those above the median in the contribution of the factors are considered important(significant)factors(Gu et al.2019)and are discussed.

    Drivers of forest resilience

    Precipitation of the wettest month was the most important factor associated with resilience.Although dry forests in the tropics are generally considered more resilient,their recovery is heavily dependent on the amount of precipitation(álvarez-Yépiz et al.2018),which is in line with the results of this study.A similar result was also reported in a wide range of tropical forest ecosystems where extended drought and low precipitation slows the recovery of forests in different continents(Verbesselt et al.2016)and Amazon mountain forests(Nobre and Borma 2009).

    Table 3 The relative influence of the variables determining resistance in Desa’a Forest(in bold are significant factors)

    Fig.10 Partial dependencies of factors affecting resistance in Desa’a forest.The relative importance of variables in the model(%out of 100)is given in brackets.Fitted functions are centred around the mean of the resilience and plotted on a common scale.Rug plots(ticks in X-axis)show the distribution of sample measurements.TS stands for temperature seasonality,MTDQ for a mean temperature of the driest quarter,HumusDh for humus depth,PWeM for precipitation of the wettest month,and TAR for temperature annual range

    Generally,tree diversity was associated with resilience,yet the Shannon and evenness indicators had a different impact.In the literature,there are contradicting findings on the effect of diversity on stability,where positive effect of species diversity has been reported in grasslands(Tilman et al.2006;Van Ruijven and Berendse 2010),and in forests across Europe(Guyot et al.2016,Sousa-Silva et al.2018,Vannoppen et al.2019),while others argue that there is no true positive diversity effect found so far on resilience(Bauhus et al.2017).We found a positive association of Shannon diversity with resilience,but saturating eventually.The positive effect of diversity on resilience might be explained by the insurance effect where different species respond differently to disturbances stabilizing the overall resilience as a system regardless of the lowered performance of certain member species(Loreau 2004).However,the effect of evenness was unimodal,with the highest evenness values resulting in a lower forest resilience.In this forest,dominant species might be needed to some extent to keep the forest community more resilient.Such species could have particular functional traits that play a significant role in the stability of the forest community(Yan et al.2011).However,diversity indices lack information to indicate the functional role of species(Yan et al.2011)and limit the identification of the species that are disadvantaged when sites get more even.In Desa’a Forest,such late successional species could be those that are less competitive such as juniper tree(Alshahrani 2008),which are disadvantaged when they grow in even proportion to others,reducing the total resilience of the forest community.

    Table 4 The relative influence of the variables determining growth stability in Desa’a forest(in bold are significant factors)

    Fig.11 Partial dependencies of factors affecting growth stability in Desa’a forest.The relative importance of variables in the model(%out of 100)is given in brackets.Fitted functions are centred around the mean of the growth stability and plotted on a common scale.Rug plots(ticks in Xaxis)show the distribution of sample measurements.PWeM stands for precipitation of the wettest month,TAR for temperature annual range,PWaQ for precipitation of the warmest quarter and DiSet for distance to settlement

    Proximity to a settlement increases the probability of anthropogenic disturbance such as grazing and cutting,which are predominant in the forest(Giday et al.2018).Our results confirm that the resilience of the vegetation located further than 5 km from settlements was considerably increased.The anthropogenic disturbance could affect resilience by affecting species composition,which might introduce an artificial dominance of a certain tree species and reduce species richness.That could have a direct impact on the resilience of the forest(Hillebrand et al.2008).

    The negative effect of slope on the resilience might be linked to its effect on soil depth,moisture content and susceptibility to degradation where steep slopes and exposed rocky areas have a little medium for plant growthdue to erosion(Zhang et al.2015),and when disturbances prevail,they are more affected than those in good soil conditions and gentle slopes.In general,in line to our hypotheses,the combination of tree diversity,local human impact,topographic position and climate(mainly precipitation)controlled resilience in dry Afromontane forest.

    Table 5 Stability metrics and their model performance(TDC is training data correlation,CVC is cross-validation correlation)

    Drivers of resistance

    While temperature seasonality was negatively associated with resistance,mean temperature of the driest quarter,humus thickness and precipitation of the wettest month was positively associated.In contrary to resilience,the resistance of forests is dependent more on their productivity before a disturbance(Wang et al.2007;Van Ruijven and Berendse 2010).Therefore,forest communities growing in productive sites,having favourable environmental conditions,are expected to show higher resistance(Wang et al.2007).In line with this argument,our results indicated that vegetation growing in sites with thicker humus and more stony sites had higher and lowered resistance,respectively.The negative effect of increased temperature seasonality on forest resistance might be a general attribute to the tropical forests which have developed themselves under relatively stable climatic conditions(Blach-Overgaard et al.2010).Therefore,in response to their narrow climatic tolerance,as the seasonality of temperature increases,forests might lose the capacity to rearrange(to adapt quickly)themselves so reducing their resilience capability(Blach-Overgaard et al.2010).Our results indicate that higher temperature seasonality and annual temperature range were associated with lower resistance.In the highland parts of Desa’a Forest,where it is relatively colder and dominated by climax species,a negative correlation between temperature and growth of juniper and olive trees was reported(Mokria et al.2017;Siyum et al.2019).Temperature seasonality between 1.8°C and 2.2°C and an annual temperature range between 21°C and 22°C were associated with higher resilience.Increased temperature seasonality and annual temperature range prolongs the disturbance and slows the recovery and break the resistance(Anjos and De Toledo 2018)due to increased fluctuation and excessive evapotranspiration(Schroth et al.2009).

    In contrast to the resilience indicator,no association between biodiversity and resistance could be found.This is in line with the findings of Van Ruijven and Berendse(2010)who reported the positive effect of biodiversity on community resilience after a drought,but there was no association found with resistance.This is another strong evidence that resistance to disturbance depends on a prior forest condition(production,health,etc.).In contrast,the post-disturbance response of the forest could be supported by its constituents,such as diversity(Van Ruijven and Berendse 2010).While our hypothesis on the positive effect of climate and good edaphic properties on resistance holds true,the effect of tree diversity was not supported by our results.

    Drivers of growth stability

    The growth stability was mainly controlled by climate,the precipitation of the wettest month.The effect of tree diversity was not observed,and only the distance to settlement as an indicator of human impact indicator was detected though not significant.In dry forests,precipitation is the most important factor for the growth of trees and increased biomass(Hiltner et al.2016).Dry forests are affected by high evapotranspiration due to the high temperature and low precipitation(Souza et al.2016),and when precipitation gets higher,the growth of the forests is positively affected.The results are in line with the findings from different tropical forests;subtropical forest in China(Gu et al.2019),dry tropical montane forests of Ethiopia(Hiltner et al.2016),and in the dry Afromontane forests(Gebru et al.2020).Effect of anthropogenic disturbances can be mediated and suppressed by the effect of precipitation which initiates more growth and system repair in forests(Rito et al.2017),which could be the reason for the non-significant effect of disturbance on the growth stability of this forest.

    The relationship among resilience,resistance and growth stability in Desa’a forest

    Forest stability was successfully characterized using resilience and resistance from remotely sensed imagery in different forests(Sousa-Silva et al.2018;Frazier et al.2018).In Desa’a,a dry tropical Afromontane forest,the three stability metrics were modeled.The correlation analysis between the metrics showed that the correlation between resilience and resistance was very weak but positive.This is in line with the concept of DeRose and Long(2014),who argued that resistance and resilience act upon ecosystems differently.While resilience is related to the influence of disturbance on the structure and composition of the ecosystem,resistance is related to the influence of the structure and composition of an ecosystem on disturbance.In support of our results,Gazol et al.(2018)reported low resistant forests to be more resilient across different biomes.Against our findings,a negative correlation was found between resistance and recovery rate from another tropical dry forest(Bhaskar et al.2018).The difference in the correlation results might be due to the difference in the interaction of climate and local degradation factors(Bhaskar et al.2018).

    Conclusion and recommendation

    The dry Afromontane forest of Desa’a was generally resistant but less resilient experiencing a continuous decline growth in stability in the last two decades.Climate variability played a pivotal role in the resilience,and resistance of the forest.While the precipitation of the wettest month is the most important factor in all the stability metrics,an inter-annual variation above 2°C is was enough to degrade the resilience and resistance of the forest.Furthermore,tree species diversity was important to enhance the resilience of the dry Afromontane forest,but no evidence of tree diversity effects was found for resistance and growth stability.We found a threshold(0.7),above which tree species evenness leads to less resilience.Experimental research might be important to investigate into what extent of evenness species identity is important to promote resilience in the dry forests.Moreover,distance to the settlement,which is an indicator of degradation and slope were also important to promote resilience.Climate,both precipitation and temperature,edaphic factors,local human disturbance indicators and tree diversity were important for one or all of the stability metrics investigated in the dry Afromontane forest.

    Supplementary Information

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

    Additional file 1:Appendix 1.Methodological protocol for humus,soil depth,and local human disturbance indicators assessment.

    Acknowledgements

    Forest inventory data were obtained from WeForest,an international nonprofit non-government organization,which is working on the restoration of Desa’a forest in collaboration with different national and international institutes(https://www.weforest.org/project/ethiopia-desaa).

    Authors’contributions

    Conceptualization:HH&BM;Methodology:HH&BM;Formal Analysis:HH&LO;Writing–original draft:HH;Writing review&editing:BM,BS,WD&KG;Supervision:BM.The author(s)read and approved the final manuscript.

    Funding

    PhD IRO grant from KU Leuven and WeForest Ethiopia supported the data collection.It is one of the chapters of a PhD research and there is no specific ID attached to the funds.

    Availability of data and materials

    The datasets generated during and/or analyzed during the current study are available in the KU Leuven repository,and are accessible according to the regulation of the University.

    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

    1KU Leuven,Department of Earth and Environmental Sciences,Division Forest,Nature and Landscape,Celestijnenlaan 200E,P.O.Box 2411,3001 Leuven,Belgium.2Department of Land Resource Management and

    Environmental Protection,Mekelle University,College of Dryland Agriculture and Natural Resources,P.O.Box 231,Mekelle,Tigray,Ethiopia.3Laboratory of Geo-Information Science and Remote Sensing,Wageningen University,6708 PB Wageningen,The Netherlands.

    Received:4 October 2020 Accepted:18 January 2021

    欧美bdsm另类| 老师上课跳d突然被开到最大视频| a级一级毛片免费在线观看| 99久久九九国产精品国产免费| 日韩人妻高清精品专区| 国产一区二区亚洲精品在线观看| 能在线免费观看的黄片| 国产亚洲精品综合一区在线观看| 夜夜看夜夜爽夜夜摸| 国产精品av视频在线免费观看| 我要搜黄色片| 久久草成人影院| 亚洲av中文字字幕乱码综合| ponron亚洲| 男人舔奶头视频| 日韩中文字幕欧美一区二区| 亚洲最大成人手机在线| 精品久久久久久久末码| 欧美另类亚洲清纯唯美| 特级一级黄色大片| 久久精品国产自在天天线| 国产精品女同一区二区软件 | 不卡一级毛片| 天堂动漫精品| 人妻少妇偷人精品九色| 99久久精品一区二区三区| 身体一侧抽搐| 久久精品影院6| 男人和女人高潮做爰伦理| 午夜久久久久精精品| 亚洲精品一卡2卡三卡4卡5卡| 亚洲精品成人久久久久久| 免费看美女性在线毛片视频| 色av中文字幕| 国产毛片a区久久久久| 亚洲精品日韩av片在线观看| 国产黄色小视频在线观看| 欧美色欧美亚洲另类二区| 成人精品一区二区免费| 人人妻人人澡欧美一区二区| 午夜亚洲福利在线播放| 中出人妻视频一区二区| 欧美日韩国产亚洲二区| avwww免费| 丝袜美腿在线中文| 九九爱精品视频在线观看| 国产毛片a区久久久久| 精品久久久久久久久亚洲 | 色综合色国产| 免费av毛片视频| 亚洲综合色惰| 国产亚洲精品av在线| 欧美潮喷喷水| 国产 一区 欧美 日韩| 国产色爽女视频免费观看| 直男gayav资源| 日本在线视频免费播放| 国产探花在线观看一区二区| 免费在线观看日本一区| 国产午夜精品论理片| 国产精品国产高清国产av| 亚洲va日本ⅴa欧美va伊人久久| 波野结衣二区三区在线| 999久久久精品免费观看国产| 一进一出抽搐gif免费好疼| 亚洲真实伦在线观看| 亚洲精品日韩av片在线观看| 国产精品,欧美在线| 成年女人永久免费观看视频| 亚洲aⅴ乱码一区二区在线播放| 精品一区二区三区视频在线| av在线亚洲专区| 亚洲自拍偷在线| 亚洲av一区综合| 中出人妻视频一区二区| 女同久久另类99精品国产91| 亚洲一级一片aⅴ在线观看| 嫁个100分男人电影在线观看| 天堂影院成人在线观看| 小说图片视频综合网站| 亚洲国产欧美人成| 又爽又黄无遮挡网站| 亚洲va日本ⅴa欧美va伊人久久| av国产免费在线观看| 欧美成人性av电影在线观看| 久久人妻av系列| 亚洲国产高清在线一区二区三| 亚洲自偷自拍三级| 狂野欧美白嫩少妇大欣赏| 人人妻人人看人人澡| 欧美日本亚洲视频在线播放| 日本免费a在线| 成人二区视频| 国产精品1区2区在线观看.| 久久婷婷人人爽人人干人人爱| 久久午夜福利片| 国产综合懂色| 精品人妻偷拍中文字幕| 我要搜黄色片| 亚洲av电影不卡..在线观看| 亚洲人成网站在线播放欧美日韩| АⅤ资源中文在线天堂| 国产精品美女特级片免费视频播放器| 亚洲精华国产精华液的使用体验 | 两个人视频免费观看高清| 又爽又黄无遮挡网站| 欧美日韩瑟瑟在线播放| 88av欧美| 欧美丝袜亚洲另类 | 中亚洲国语对白在线视频| 国产在视频线在精品| 亚洲电影在线观看av| 国产成人福利小说| 国产精品久久久久久久久免| 国产精品亚洲一级av第二区| 国产av不卡久久| 精品一区二区三区视频在线| 黄色欧美视频在线观看| 1000部很黄的大片| 一夜夜www| 色av中文字幕| 欧美三级亚洲精品| 国产精品爽爽va在线观看网站| 免费看av在线观看网站| 国内精品一区二区在线观看| 亚洲av免费在线观看| 熟女人妻精品中文字幕| 国产亚洲精品av在线| 美女大奶头视频| 午夜福利在线观看吧| 观看免费一级毛片| 人妻夜夜爽99麻豆av| 亚洲最大成人av| 欧美xxxx性猛交bbbb| 色5月婷婷丁香| 午夜福利欧美成人| 最近在线观看免费完整版| 亚洲中文日韩欧美视频| 久久亚洲精品不卡| 午夜福利视频1000在线观看| 婷婷丁香在线五月| 干丝袜人妻中文字幕| 亚洲精品粉嫩美女一区| 免费观看人在逋| 成人无遮挡网站| 狂野欧美白嫩少妇大欣赏| 在线播放国产精品三级| av在线亚洲专区| 人妻夜夜爽99麻豆av| 久久精品夜夜夜夜夜久久蜜豆| 国产精品三级大全| 97热精品久久久久久| 国产单亲对白刺激| 国产亚洲精品综合一区在线观看| 国产女主播在线喷水免费视频网站 | 一个人看的www免费观看视频| 一进一出好大好爽视频| 久久精品国产亚洲网站| 黄色日韩在线| 舔av片在线| 精品日产1卡2卡| 成人特级黄色片久久久久久久| 久久久久久国产a免费观看| 最新在线观看一区二区三区| 国产精品久久久久久久电影| 校园人妻丝袜中文字幕| 精品人妻一区二区三区麻豆 | 深夜a级毛片| 亚洲精华国产精华精| 国内精品久久久久精免费| 婷婷精品国产亚洲av在线| 69av精品久久久久久| 久久人妻av系列| 日韩欧美国产一区二区入口| a在线观看视频网站| 精品久久久久久久人妻蜜臀av| www日本黄色视频网| 俺也久久电影网| 在线国产一区二区在线| 精品乱码久久久久久99久播| 国产精品久久久久久久久免| 亚洲专区中文字幕在线| 久久久久国内视频| 欧美绝顶高潮抽搐喷水| 国产国拍精品亚洲av在线观看| 亚洲成a人片在线一区二区| 亚洲av二区三区四区| 国产欧美日韩精品亚洲av| 国产精品伦人一区二区| 老熟妇仑乱视频hdxx| 国产伦人伦偷精品视频| 亚洲在线观看片| 日韩精品青青久久久久久| 国产精品三级大全| 两个人的视频大全免费| 大型黄色视频在线免费观看| 草草在线视频免费看| 床上黄色一级片| 亚洲五月天丁香| 成年免费大片在线观看| 久久久久久久亚洲中文字幕| 白带黄色成豆腐渣| 国产精品无大码| 国国产精品蜜臀av免费| 久久九九热精品免费| 黄色丝袜av网址大全| 欧美精品啪啪一区二区三区| 国产男靠女视频免费网站| 麻豆国产97在线/欧美| 日本欧美国产在线视频| 观看免费一级毛片| 女同久久另类99精品国产91| 99久久中文字幕三级久久日本| 嫁个100分男人电影在线观看| 国产综合懂色| 亚洲成人免费电影在线观看| 国产精品三级大全| 亚洲天堂国产精品一区在线| 18禁黄网站禁片免费观看直播| 欧美色视频一区免费| 亚洲va日本ⅴa欧美va伊人久久| 深夜精品福利| 免费黄网站久久成人精品| 成人毛片a级毛片在线播放| а√天堂www在线а√下载| 久久久精品大字幕| 亚洲综合色惰| 色视频www国产| 日韩中文字幕欧美一区二区| 91久久精品国产一区二区成人| 国产精品电影一区二区三区| 日日撸夜夜添| 91麻豆av在线| 国产高清视频在线播放一区| 国模一区二区三区四区视频| 看免费成人av毛片| 国产蜜桃级精品一区二区三区| 成年女人永久免费观看视频| 少妇的逼好多水| 天美传媒精品一区二区| 亚洲成人久久爱视频| 老师上课跳d突然被开到最大视频| 少妇高潮的动态图| 国产精品综合久久久久久久免费| 国产成年人精品一区二区| 亚洲四区av| 精品久久久久久久久久免费视频| 精品一区二区三区人妻视频| 久久久久久国产a免费观看| 少妇人妻精品综合一区二区 | 人人妻,人人澡人人爽秒播| 乱人视频在线观看| 欧美3d第一页| 国产黄片美女视频| 国产精品免费一区二区三区在线| 麻豆成人午夜福利视频| 日韩一区二区视频免费看| 欧美bdsm另类| 成年女人毛片免费观看观看9| 我的老师免费观看完整版| 夜夜爽天天搞| 国产精品人妻久久久久久| 国产精品国产三级国产av玫瑰| 久久草成人影院| 日韩人妻高清精品专区| 久久久国产成人免费| 三级男女做爰猛烈吃奶摸视频| 黄色视频,在线免费观看| 国内精品久久久久精免费| 日韩精品有码人妻一区| 女的被弄到高潮叫床怎么办 | 特大巨黑吊av在线直播| 亚洲精品一区av在线观看| 在现免费观看毛片| 午夜日韩欧美国产| 中文亚洲av片在线观看爽| 18+在线观看网站| 国产三级中文精品| 丰满人妻一区二区三区视频av| 国产精华一区二区三区| 日韩一区二区视频免费看| 亚洲精品一区av在线观看| 久久6这里有精品| 日本三级黄在线观看| 中文亚洲av片在线观看爽| 精品久久久久久久久av| 日本撒尿小便嘘嘘汇集6| 国产精品一区二区三区四区免费观看 | 色综合婷婷激情| 日本色播在线视频| 国产精华一区二区三区| 亚洲中文日韩欧美视频| 精品久久久久久久久亚洲 | 久久精品久久久久久噜噜老黄 | 精品一区二区免费观看| 特大巨黑吊av在线直播| 国产av在哪里看| 1024手机看黄色片| 国产一区二区三区视频了| 欧美xxxx黑人xx丫x性爽| 久久久久久久久久成人| 我要看日韩黄色一级片| 免费在线观看日本一区| 久久中文看片网| 精品久久久久久成人av| 亚洲欧美日韩无卡精品| x7x7x7水蜜桃| 男人舔奶头视频| 99在线视频只有这里精品首页| 在线a可以看的网站| 亚洲欧美日韩卡通动漫| 国产毛片a区久久久久| 日韩欧美免费精品| 国产免费男女视频| 午夜亚洲福利在线播放| 国产探花极品一区二区| 听说在线观看完整版免费高清| 免费一级毛片在线播放高清视频| 99在线人妻在线中文字幕| av黄色大香蕉| 少妇丰满av| 免费看日本二区| 国产男靠女视频免费网站| 日韩 亚洲 欧美在线| 国产一级毛片七仙女欲春2| 国产精品久久视频播放| 免费av不卡在线播放| 亚洲国产精品成人综合色| 国产精品亚洲美女久久久| 噜噜噜噜噜久久久久久91| 啦啦啦观看免费观看视频高清| 国产私拍福利视频在线观看| 十八禁国产超污无遮挡网站| 美女黄网站色视频| 99热这里只有是精品50| 十八禁网站免费在线| 老师上课跳d突然被开到最大视频| 精品国产三级普通话版| 亚洲第一电影网av| 男女做爰动态图高潮gif福利片| 亚洲精品456在线播放app | 免费看av在线观看网站| 亚洲精品一卡2卡三卡4卡5卡| 成人永久免费在线观看视频| 国产 一区精品| 日韩欧美精品免费久久| 国产色婷婷99| 亚洲精品一区av在线观看| 久久久久久久精品吃奶| 一本精品99久久精品77| 97碰自拍视频| 天天躁日日操中文字幕| 成人欧美大片| avwww免费| 欧美xxxx性猛交bbbb| 国产精品98久久久久久宅男小说| 有码 亚洲区| 精品99又大又爽又粗少妇毛片 | 岛国在线免费视频观看| 中文在线观看免费www的网站| 男人舔奶头视频| 亚洲精品久久国产高清桃花| 亚洲精华国产精华液的使用体验 | 尾随美女入室| 欧美黑人欧美精品刺激| 午夜精品一区二区三区免费看| 亚洲18禁久久av| 国产高清三级在线| 亚洲18禁久久av| 国内精品久久久久精免费| 波多野结衣高清作品| 国内精品久久久久精免费| 欧美精品国产亚洲| 精品久久久噜噜| 国产欧美日韩精品亚洲av| 大又大粗又爽又黄少妇毛片口| 波多野结衣高清作品| 中文字幕久久专区| 久久久久久国产a免费观看| 午夜视频国产福利| 美女xxoo啪啪120秒动态图| 国产伦一二天堂av在线观看| 深爱激情五月婷婷| 欧美激情久久久久久爽电影| 精品久久久久久久末码| 精品一区二区三区视频在线观看免费| 国产淫片久久久久久久久| 日韩欧美 国产精品| 亚洲成人久久爱视频| 91午夜精品亚洲一区二区三区 | 热99在线观看视频| 最近最新中文字幕大全电影3| 亚州av有码| 人妻夜夜爽99麻豆av| 欧美黑人巨大hd| 变态另类成人亚洲欧美熟女| 亚洲av免费高清在线观看| 亚洲中文字幕一区二区三区有码在线看| 男人和女人高潮做爰伦理| 一区二区三区激情视频| 久久午夜亚洲精品久久| 22中文网久久字幕| aaaaa片日本免费| 亚洲自拍偷在线| 亚洲欧美精品综合久久99| 精华霜和精华液先用哪个| 亚洲男人的天堂狠狠| 久久久成人免费电影| 尤物成人国产欧美一区二区三区| 国产精品一区二区三区四区久久| 久久久久久久久中文| 我要看日韩黄色一级片| 久久久久国产精品人妻aⅴ院| 嫩草影院新地址| 又黄又爽又免费观看的视频| 99riav亚洲国产免费| 欧美日韩黄片免| 简卡轻食公司| 午夜影院日韩av| 婷婷丁香在线五月| 91麻豆精品激情在线观看国产| 精品一区二区三区av网在线观看| 国产一区二区亚洲精品在线观看| 国产成人一区二区在线| 国产伦人伦偷精品视频| 三级男女做爰猛烈吃奶摸视频| 午夜福利视频1000在线观看| 男女啪啪激烈高潮av片| 91精品国产九色| 久久精品久久久久久噜噜老黄 | a级毛片免费高清观看在线播放| 韩国av一区二区三区四区| 日本爱情动作片www.在线观看 | 免费av毛片视频| 国产成人a区在线观看| 国产一区二区激情短视频| 中文资源天堂在线| 欧美黑人巨大hd| 亚洲电影在线观看av| 国产精品一区www在线观看 | 一本精品99久久精品77| 有码 亚洲区| 久久久久免费精品人妻一区二区| 久久久久久久精品吃奶| 亚洲色图av天堂| 亚洲性夜色夜夜综合| 成年免费大片在线观看| 精品久久久久久成人av| 亚洲av中文av极速乱 | 中文字幕高清在线视频| 日日撸夜夜添| 男人狂女人下面高潮的视频| 久久欧美精品欧美久久欧美| 久久婷婷人人爽人人干人人爱| 午夜激情欧美在线| 亚洲一区二区三区色噜噜| 亚洲aⅴ乱码一区二区在线播放| 中文字幕高清在线视频| 免费高清视频大片| 精品久久久噜噜| 免费电影在线观看免费观看| 别揉我奶头 嗯啊视频| 欧美bdsm另类| 亚洲va在线va天堂va国产| 亚洲av不卡在线观看| 久久中文看片网| 色视频www国产| 一进一出抽搐gif免费好疼| 久久天躁狠狠躁夜夜2o2o| 国内精品一区二区在线观看| 欧美性感艳星| 中文字幕久久专区| 精品久久国产蜜桃| 日韩一区二区视频免费看| 久久久久久国产a免费观看| 免费av观看视频| 男女边吃奶边做爰视频| 欧美一级a爱片免费观看看| 色综合亚洲欧美另类图片| 久久精品国产鲁丝片午夜精品 | 美女高潮的动态| 日韩欧美国产一区二区入口| 中文字幕精品亚洲无线码一区| 大型黄色视频在线免费观看| 精品人妻熟女av久视频| 男女啪啪激烈高潮av片| 亚洲成人久久爱视频| 男女啪啪激烈高潮av片| 免费在线观看成人毛片| 91午夜精品亚洲一区二区三区 | 国产乱人视频| 乱系列少妇在线播放| 国产乱人视频| 51国产日韩欧美| 22中文网久久字幕| 欧美xxxx性猛交bbbb| 成人国产一区最新在线观看| bbb黄色大片| 999久久久精品免费观看国产| 男人和女人高潮做爰伦理| 日本黄大片高清| 99精品久久久久人妻精品| 日本免费一区二区三区高清不卡| 国产精品一区www在线观看 | 久久午夜亚洲精品久久| 国内精品美女久久久久久| 精品欧美国产一区二区三| 人妻夜夜爽99麻豆av| 真人一进一出gif抽搐免费| 老熟妇仑乱视频hdxx| 很黄的视频免费| 99国产极品粉嫩在线观看| 国产视频内射| 欧美日韩精品成人综合77777| 亚洲av不卡在线观看| 一a级毛片在线观看| 亚洲av电影不卡..在线观看| 搡女人真爽免费视频火全软件 | 亚洲欧美激情综合另类| 亚洲成人精品中文字幕电影| 亚洲欧美精品综合久久99| 3wmmmm亚洲av在线观看| 欧美另类亚洲清纯唯美| 草草在线视频免费看| 老女人水多毛片| 精品久久久久久久久av| 啦啦啦韩国在线观看视频| 十八禁网站免费在线| 看片在线看免费视频| 极品教师在线免费播放| 91狼人影院| 两个人视频免费观看高清| 国产黄片美女视频| 久久久久国内视频| 午夜福利在线观看免费完整高清在 | 日日摸夜夜添夜夜添av毛片 | 伦理电影大哥的女人| 中亚洲国语对白在线视频| 成年女人看的毛片在线观看| 日韩亚洲欧美综合| 国产精品综合久久久久久久免费| 国产麻豆成人av免费视频| 国产三级中文精品| 少妇人妻精品综合一区二区 | 国产精品嫩草影院av在线观看 | 露出奶头的视频| 色5月婷婷丁香| 欧美在线一区亚洲| 国产高清不卡午夜福利| 午夜福利视频1000在线观看| 色哟哟哟哟哟哟| 三级毛片av免费| 亚洲自拍偷在线| 国产v大片淫在线免费观看| 久久这里只有精品中国| 尤物成人国产欧美一区二区三区| 色吧在线观看| 欧美最黄视频在线播放免费| 女生性感内裤真人,穿戴方法视频| 国产精品电影一区二区三区| 久久久久久大精品| 一个人看视频在线观看www免费| а√天堂www在线а√下载| 国产精品国产高清国产av| 美女xxoo啪啪120秒动态图| 久久久久免费精品人妻一区二区| 成人午夜高清在线视频| 伦理电影大哥的女人| 最后的刺客免费高清国语| 我的老师免费观看完整版| 午夜激情欧美在线| 亚洲国产精品sss在线观看| 国产久久久一区二区三区| 美女高潮喷水抽搐中文字幕| 欧美潮喷喷水| 黄色视频,在线免费观看| 欧美激情久久久久久爽电影| 成人无遮挡网站| 国产高潮美女av| 午夜激情福利司机影院| 亚洲自拍偷在线| 欧洲精品卡2卡3卡4卡5卡区| 深夜精品福利| 国产在线男女| 搡老妇女老女人老熟妇| 婷婷精品国产亚洲av| 春色校园在线视频观看| 三级国产精品欧美在线观看| 中文字幕精品亚洲无线码一区| 91av网一区二区| 直男gayav资源| 天堂√8在线中文| 午夜免费男女啪啪视频观看 | 久久午夜福利片| 人妻夜夜爽99麻豆av| 桃色一区二区三区在线观看| 一a级毛片在线观看| 久久国产精品人妻蜜桃| 精品久久国产蜜桃| 国产精品99久久久久久久久| 欧美激情在线99| 欧美日韩综合久久久久久 | 国产亚洲精品av在线| 亚洲黑人精品在线| 在线a可以看的网站| 好男人在线观看高清免费视频| 午夜免费成人在线视频| 99久国产av精品| 色吧在线观看| 人人妻人人看人人澡| 欧美在线一区亚洲| av在线蜜桃| 麻豆精品久久久久久蜜桃| 婷婷六月久久综合丁香| 国产高清三级在线|