999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

四川省岷江冷杉對氣候變化的響應及其潛在分布格局

2022-06-27 14:05:34潘少安李旭華馮秋紅劉興良孫建新
生態學報 2022年10期
關鍵詞:物種模型

潘少安,李旭華,馮秋紅,劉興良,孫建新

1 北京林業大學生態與自然保護學院,北京 100083 2 四川省林業科學研究院,森林和濕地生態恢復與保育四川省重點實驗室,成都 610081 3 四川臥龍森林生態系統國家定位觀測研究站,成都 610081

多年前,科學家們就已證實,全球氣候變化已經無可爭辯[1]。當前的氣候變化趨勢和氣候預測均表現出溫度上升隨氣候變化的一致趨勢[2]。陸地生態系統對溫度變化高度敏感,且溫室氣體的增加使陸地生態系統面臨發生重大轉變的危險[3]。近幾十年來,氣候變暖導致亞高山針葉林分布向更高海拔轉變[4],而水分的減少限制了亞高山樹木的分布范圍,同時還可能導致林木數量的減少[5]。有研究表明峨眉冷杉在海拔3000 m處分布面積逐漸減少最終消失并被其它物種取代,而在3600 m處以分散分布出現,最終成為優勢物種[6]。亞高山森林分布的動態變化對氣候變化具有十分重要的指示作用[7],尤其是冷杉屬物種,對氣候變化十分敏感,對溫濕度穩定性有著嚴格要求[8],因此,研究岷江冷杉對氣候的響應對于理解全球變化具有極其重要的意義。

岷江冷杉(Abiesfaxoniana)為松科冷杉屬常綠喬木,是我國特有樹種,也是川西亞高山地區森林的優勢樹種之一,通常在海拔2800 m以上的區域形成純林,或與紅杉、方枝柏等形成混交林[9],是岷江上游地區陰坡林線的主要建群樹種[10]。岷江冷杉林是四川省面積最大的原始森林類型和大熊貓的重要棲息場所,在涵養水源、保持水土方面具有十分重要的功能,亦維持著大熊貓棲息地生態系統的穩定與安全[11—12]。目前關于岷江冷杉的研究較為豐富,涉及岷江冷杉林的群落動態變化[13]、凋落物分解[14]、年輪-氣候響應研究[15—17]、林地水源涵養[18]以及天然更新[19—20]等多個方面。在全球氣候變化背景下,岷江冷杉種群的潛在分布格局及其對未來氣候變化的響應方面的研究卻鮮有報道。

近年來,研究表明森林群落隨著氣候變暖而向上遷移,但很少提及這些變化的具體方式和具體方法,而植物在應對氣候變化的過程中,遷移或擴散能力是植物生存的先決條件[21]。我國學者在對目標物種的潛在分布區預測,以及物種分布對氣候變化的響應等方面已開展了諸多研究[22],而物種分布模型即為最常用的研究方法,該方法的基本原理是利用觀測到的物種分布數據及其對應的環境變量信息,關聯性構建二者的相關模型,進而確定目標物種的實際分布并對物種的潛在分布進行預測[23]。目前常用的物種分布模型有MaxEnt模型、Garp模型、Bioclim模型、Dumain模型、廣義線性模型、廣義可加模型等[24],但MaxEnt模型以其運行樣本量要求低、易操作、運行時間短、數據處理能力強、預算精度高等優點[25—26]而得到廣泛應用。因此,本研究采用物種分布模型—最大熵模型(MaxEnt)和ArcGIS軟件,對岷江冷杉在當代氣候條件和未來氣候變化情景下的潛在分布區和生境適宜性進行評估和分析,并確定影響岷江冷杉空間分布格局的主要影響因子和分布閾值,分析岷江冷杉對未來氣候變化的響應策略,為川西山地岷江冷杉天然更新、大熊貓棲息地恢復、退化森林質量精準提升等提供理論依據。

1 材料與方法

1.1 岷江冷杉分布數據的采集

岷江冷杉的地理分布數據主要從3個方面獲取:包括岷江冷杉樣地調查數據、公開發表文獻資料、權威的標本信息網絡數據庫,如全球生物多樣性信息網GBIF(https://www.gbif.org),中國數字植物標本館等(http://www.cvh.org.cn)。對搜集到的物種分布數據進行整理,剔除重復的和研究區域范圍外的坐標點,共得到93個參與模型運算的物種分布點,并按要求將其另存為csv格式。

1.2 環境數據的收集與預處理

本文所用的環境數據包括氣候數據、土壤數據和地形數據,其中氣候數據來源于世界氣候數據庫WorldClim(https://www.worldclim.org/),分別選取當代(1970—2000年)氣候數據和未來兩個時期(2050s和2070s)3種氣候變化情景RCP2.6,RCP4.5和RCP8.5的生物氣候變量,3種氣候情景分別表示未來時期CO2排放濃度低、中、高的3種可能變化路徑[27],每個時期的氣候數據變量均為19個。土壤數據從寒區旱區科學數據中心提供的世界土壤數據庫(Harmonized World Soil Database,HWSD)中提取出研究區域與岷江冷杉生長相關的土壤因子,如表層土壤容重、pH值、沙粒含量、粉粒含量、粘土含量、碎石百分比、有機碳含量、基本飽和度、碳酸鹽含量等。地形數據包括海拔、坡度和坡向,其中海拔數據(空間分辨率1 km)從全球海拔數據庫(https://globalmaps.github.io/)中下載并提取,坡度和坡向則分別根據海拔數據運用ArcGIS中的空間分析工具計算得到。對上述所有環境因子進行研究區掩膜提取和空間處理,使其具有統一的空間坐標系(GCS-WGS84)和相同的空間分辨率(30 arc sec×30 arc sec,約1 km)。

1.3 環境因子的相關性分析

研究指出環境因子間的多重共線性會導致物種分布模型的過度擬合[28—29],因此在運行模型前,對環境因子分別進行Pearson相關性分析,進而篩選出相關性較低且對物種分布有重要影響的環境因子(相關系數|r|<0.8)來參與模型運算。最終共選取了16個環境因子參與岷江冷杉分布區預測與生境評價,包括6個氣候因子,7個土壤因子和3個地形因子。具體指標情況見表1。

表1 用于MaxEnt模型的環境變量

1.4 模型運行與準確性

將符合MaxEnt模型格式要求的岷江冷杉物種分布數據,分別與當前和未來兩個時期不同氣候情景下的環境變量數據導入模型中,隨機選取25%的分布點作為測試集,75%的物種分布點作為訓練集,模擬共得到7個岷江冷杉適生分布區的預測結果。模型預測結果的精確度則通過受試者工作特征曲線(Receiver Operating Characteristic curve,ROC曲線)下面積(Area Under Curve,AUC)進行評價[30],該值大小與模型的預測精度呈正相關[31],當AUC值小于0.6時模型預測失敗,當AUC值介于0.8—0.9表明模型預測精度較好,達到0.9—1.0時則模型預測精度很高[32]。同時,模型通過貢獻率、置換重要值和刀切法檢驗(Jackknife test)3個方面來綜合評估個各環境變量在影響岷江冷杉地理分布中的重要性。

1.5 適宜分布區劃分及未來適生區空間變化分析

模型的預測結果是通過物種在待預測地區的存在概率P(取值0—1)來表示物種的分布適宜性[33],將模型運行得到的分布區預測結果進行適宜性等級劃分和可視化表達,結合岷江冷杉的生態學特點,采用重分類中的人工(Manual)分級方法,將岷江冷杉的潛在分布區劃分為4個等級,分別是:適宜性指數在0—0.1為非適生區,0.1—0.3為低適生區,0.3—0.5為中適生區,0.5—1.0為高適生區。為了進一步分析岷江冷杉在未來氣候變化情景下適宜分布區的空間變化格局,對比當代適生區分布結果,我們分別提取了新增適生區、喪失適生區、保留適生區和非適生區4種岷江冷杉適生區的變化類型。

2 結果與分析

2.1 模型的預測精度

圖1 岷江冷杉分布區預測的受試者工作曲線 Fig.1 Receiver operating characteristic (ROC) of Abies faxoniana AUC: 曲線下面積 Area under curve

模型運行結果顯示,訓練集和測試集的AUC值分別為0.955和0.853(圖1),達到了模型預測精度要求,由此表明該模型對岷江冷杉潛在分布區預測結果的可信度較高。

2.2 當代氣候條件下岷江冷杉的潛在分布區預測

對預測的岷江冷杉潛在適生區等級劃分的結果如圖2所示,經空間分析模塊統計各適生區的分布面積,結果表明岷江冷杉中適生分布區和高適生分布區的面積約為33176.86 km2和31416.84 km2,分別占全省總面積的5.91%和5.60%,低適生區占全省總面積的12.30%,而非適生區的占比最高為76.19%。岷江冷杉的中、高適生區集中分布于阿壩藏族羌族自治州汶川縣、理縣、黑水縣、松潘縣等地;雅安市的天全縣、寶興縣、綿陽市的平武縣、北川縣等地,除此之外在成都市大邑縣、甘孜藏族自治州的丹巴縣、康定縣等地的零星的面積也是岷江冷杉的潛在適生分布區。

圖2 當代氣候條件下岷江冷杉適生等級地理分布格局Fig.2 Geographic distribution pattern of suitable grade for Abies faxoniana in contemporary China

2.3 環境變量貢獻率及主導因子篩選

在影響岷江冷杉地理分布的環境因子中,氣候因子的累計貢獻率為74.16%,地形因子和土壤因子的累計貢獻率分別為13.33%和12.51%,而這三類因子對應的置換重要性則分別為84.79%,12.25%和2.96%(表2),由此表明岷江冷杉地理分布主要受氣候因子的影響,而地形對岷江冷杉地理分布的影響次之。

具體來講,對岷江冷杉地理分布貢獻率較高的前4個環境變量分別是降水季節性變異系數、氣溫年變化幅度、晝夜溫差月均值和海拔,其累計貢獻率為80.77%;置換重要值較高的環境因子分別是氣溫年變化幅度、年降水量、海拔和降水季節性變異系數,其累計值達到了84.99%(表2)。此外,根據刀切法檢驗結果,在僅使用單變量時的正規化訓練增益最高的3個變量分別是氣溫年變化幅度、年降水量、海拔,而在使用除單個環境變量外的其他變量時,氣溫年變化幅度和降水季節性變異系數的增益值降低幅度最大(表2),由此可見,這二者具有其他環境變量中沒有的有效信息。綜合上述3類評估結果,影響岷江冷杉地理分布的氣候因子主要是降水季節性變異系數、氣溫年變化幅度和年降水量,地形因子則主要是海拔。

2.4 主要貢獻率環境變量閾值分析

根據環境因子響應曲線,計算影響岷江冷杉地理分布的各主導因子的閾值(即物種存在概率>0.5的范圍)。結果顯示(圖3),岷江冷杉適宜生長的海拔范圍是2310—3757 m;氣溫年變化幅度的適宜范圍為29.3—32.5 ℃;降水季節性變異系數的適宜范圍為76.5—85.4,最適宜為81.7左右;而適宜岷江冷杉生長的年降水量的范圍是694.1—791.7 mm。

表2 影響岷江冷杉分布的環境因子貢獻率

圖3 主要環境因子的響應曲線Fig.3 Response curves of four major environmental factors

2.5 未來氣候變化情景下岷江冷杉適生區及其空間格局變化

在未來兩個時期(2050s和2070s)的不同氣候變化情景下,岷江冷杉潛在中適生區和高適生區的面積較當代氣候條件下適生區面積均有所增加,且潛在高適生區面積的增加幅度總體上高于中適生區面積的增加幅度。此外,未來兩個時期的中適生區的面積占比均隨著排放濃度的增加而逐漸減少,而高適生區的面積占比則表現出先增加后降低的變化(表3),即在中等排放濃度時的潛在高適生區面積最高。

表3 不同氣候模式下岷江冷杉潛在分布區的面積占比/%

岷江冷杉適生分布區的空間變化格局顯示,在未來不同氣候變化情景下,岷江冷杉的潛在適生分布區總體向南移,新增適生區多分布于盆周山地,而喪失適生區集中分布于原適生區北緣和西北緣(圖4)。其中在2050s三種排放情景下,岷江冷杉潛在適生區的新增率在RCP8.5排放情景下最高,為34.27%,其次是RCP2.6排放情景下的34.20%(表4);而在2070s,岷江冷杉新增適生區面積在RCP8.5排放情景下最少,新增率為31.62%,在RCP4.5排放情景下的新增適生面積最大,新增率為39.79%。岷江冷杉潛在適生區在未來兩個時期不同氣候變化情景下,喪失面積最大的是2070s的RCP8.5排放情景,喪失面積最小的是在2070s的RCP4.5排放情景。

表4 未來氣候變化情景下岷江冷杉適生區的變化率/%

圖4 不同氣候變化情景下岷江冷杉適生區空間變化格局Fig.4 Spatial variation patterns of suitable areas of Abies faxoniana under different climate change scenariosRCP: 典型濃度路徑 representative concentration pathway

3 討論

3.1 模型模擬可靠性及岷江冷杉的潛在分布區分析

當前,利用物種分布模型來預測物種的地理分布已開展了很多研究,在眾多物種分布模型中,最大熵模型(MaxEnt)是使用最為廣泛且具有較好預測能力的模型[22]。本研究基于MaxEnt模型和ArcGIS空間分析,定量展示了當代氣候變化條件下岷江冷杉在四川省境內的適生分布區,預測結果顯示岷江冷杉潛在分布區要比實際已知的生長區域范圍大,主要分布于阿壩藏族羌族自治州、雅安市中部和北部以及甘孜藏族自治州丹巴縣、康定縣等地的零星地區,這些地區所在的森林分區為川西亞高山云冷杉林區,位于青藏高原東南緣,也是金沙江、大渡河、岷江、涪江等大江大河主要的水源涵養區或源頭。適生分布區的氣候受高原地形的決定性影響,形成冬寒夏涼的高山氣候特點,這也符合岷江冷杉耐陰性強,喜冷濕氣候的生物學特性。研究結果表明,當代氣候條件下,岷江冷杉潛在適生區(中適生區和高適生區)的面積占研究區總面積的11.51%,說明岷江冷杉的生態位較為狹窄。此外,模型模擬結果經ROC曲線精度檢驗,AUC值高于0.8,且預測的潛在分布區與實際調查分布基本一致,由此表明了MaxEnt模型對岷江冷杉的預測精度較好,較為真實地反映了岷江冷杉潛在適生區的分布情況,但模擬的精度尚未達到很高,這可能與物種分布數據量不足有關,模型的表現與物種分布數據的數量正相關[34]。

3.2 影響岷江冷杉地理分布的主導因子

在不同的環境尺度下, 影響物種分布的各類因素的作用程度是不同的。Soberon和Peterson[35]將影響物種分布的因素總結為4類:(1)非生物因素,包括氣候、土壤條件、地形、水文等;(2)生物因素,包括種間競爭、取食、互利共生等;(3)該地物種自身遷移能力和地理區域的特性;(4)物種或種群對新環境的適應能力。通常,非生物因素主要在大尺度空間影響物種的分布,這些因素很大程度上決定了物種的分布范圍和格局,包括生理制約、物種對氣候和生境梯度的響應和選擇等。劉曉彤等[22]文中指出,統計自2000年以來有關物種分布模擬的研究中,有50%以上的研究只考慮了氣候因素,涉及地形和土壤及其他因子的研究較少,本研究綜合考慮了氣候、地形及土壤因子對岷江冷杉潛在分布區的影響。結果表明,氣候因子對岷江冷杉地理分布預測的貢獻最大,土壤因子和地形對分布預測的貢獻率較小。后二者可能在小尺度范圍內影響植物的分布,但在區域尺度上的影響較小,這可能是由于地形和土壤因子對植物的影響是間接的[36]。影響岷江冷杉分布的降水因子主要是降水季節性變異系數,這一氣候因子反映了岷江冷杉在不同生長期對降水的需求不同,研究表明降水量的季節性變異系數在76.5—85.4的區域為岷江冷杉的適宜生長范圍。影響岷江冷杉分布的溫度因子是年均氣溫變化幅度,其適宜范圍是29.3—32.5 ℃,即在年均氣溫變幅較高的地區適合岷江冷杉分布。上述兩個氣候因子的主導性反映了氣候條件的變異性決定了岷江冷杉的分布格局。在地形因子中岷江冷杉受海拔的影響最大,在海拔2310—3757 m區間是岷江冷杉的適宜生長范圍,而《四川松杉植物植物地理》[37]指出岷江冷杉垂直分布于海拔2400—4000 m,坡度>25°的陰坡、半陰坡的高山地帶,本文結果符合岷江冷杉的實際生長環境;雖然海拔適宜范圍的最低閾值低于2400 m,說明在海拔2300—2400 m的亞高山地帶也是岷江冷杉的潛在適生分布區,在未來氣候變化或人工種植的情況下岷江冷杉是可以存活的。本研究中參與模型構建的環境數據只涉及了非生物因素,未考慮岷江冷杉自身適應能力、種間相互作用及人類干擾數據,因而預測結果與實際情況也會存在一定的偏差。

3.3 岷江冷杉地理分布對未來氣候變化的響應

在全球氣候變化背景下,預測物種在未來氣候變化下的分布變化趨勢,對物種多樣性保護具有重要意義[27]。岷江冷杉對氣候變化的響應表現為:在未來兩個時期的不同氣候變化情景下,其潛在適生區面積均較當代氣候條件下的適生區面積有所增加,但中等適生區的面積占比均隨著排放濃度(低、中、高)的增加而逐漸減少,而潛在高適生區面積在中等濃度時最高,表明不同生境中的岷江冷杉對不同CO2濃度以及溫度的變化產生不同的響應,在同一個時期內隨著CO2濃度的升高,一定溫度范圍內,溫度升高很大程度上將抑制低海拔岷江冷杉的生長,從而導致岷江冷杉部分生境喪失;但在高濃度排放情景下,高海拔地區的溫度升高明顯,使得岷江冷杉有向高海拔低溫區擴張的可能[38],高海拔地區升高的溫度有利于岷江冷杉形成層更早的活動,表明高海拔地區將由不適宜生境逐漸向適宜生境轉變;在未來更長的時間內,中等CO2濃度排放情景更有利于岷江冷杉林生境的擴張,低CO2濃度排放情景和高CO2濃度排放情景下,升溫幅度低于2 ℃和升溫幅度太高都對岷江冷杉的生長不利。

未來氣候變化情景下岷江冷杉的潛在適生分布區總體向南遷移,這與Liao[8]等人的研究結果一致,該研究表明在未來氣候變化情景下,岷江冷杉的潛在適生區會沿著岷江流域向南擴張至云南省的東北部。綜合各個排放情景,生境喪失區域主要分布在阿壩藏族羌族自治州的若爾蓋縣北部、阿壩縣,以及雅安市、樂山市、眉山市以及涼山彝族自治州交界的馬爾康縣、天全縣、滎經縣、瀘定縣、石棉縣等。根據生境適宜分布區預測結果(圖2),這些分布區主要位于非適生區和低適生區的交界處,處于岷江冷杉分布的林線下緣。全球變暖背景下,中低海拔岷江冷杉存在季節性水分脅迫,抑制岷江冷杉生長,尤其是對大齡樹木的影響更加強烈[15]。此外,物種生長過程中的局部環境的生存競爭壓力也發揮著不可忽視的作用。高海拔地區氣溫逐漸升高,原本不利生境條件得到改善,林線下緣岷江冷杉面臨更大的資源競爭壓力[39]。較大的生存壓力對岷江冷杉種群幼苗的存活極為不利,以至于岷江冷杉種群幼苗死亡率較高,不利于種群的更新與擴張[10]。四川省地形復雜多樣,海拔垂直梯度較大,復雜的地形地貌條件導致局部氣候變化多樣,四川省生物氣候類型未來轉變趨勢由溫度低的類型像溫度高的類型轉變[40],低海拔一些原本不適宜岷江冷杉分布的區域,可能轉化成潛在適宜分布區。

對于物種的新增適生區,表明未來氣候變化使得這些區域的氣候條件更加適宜該物種的生長,物種分布的可能性增加,但實際的分布格局還要取決于物種向新增分布區遷移的可能性[41]。未來氣候變化條件下的新增適生區能夠為岷江冷杉未來的更新與恢復提供指引,可以通過人為的方式在這些區域引種或培育岷江冷杉,增加其遷入的可能性。而對于喪失適生區,可以通過采取遷地保護、提高岷江冷杉幼苗繁育等措施來提高其應對氣候變化的能力。未來氣候變化條件下,保留適生區是物種應對氣候變化的安全地和避難所[42],因而更應該注重和加強該區域岷江冷杉的保護。可以看到,未來不同氣候變化情景下岷江冷杉適生區的保留率基本保持在80%左右,說明岷江冷杉應對氣候變化的能力較強。

4 結論

(1)最大熵模型MaxEnt能夠較好地預測岷江冷杉的潛在適生分布(AUC值>0.8)。在當代氣候條件下,岷江冷杉的中、高適生區集中分布于阿壩藏族羌族自治州、雅安市中部和北部、綿陽市的平武縣、北川縣等地,分別占全省總面積的5.91%和5.60%。

(2)影響岷江冷杉地理分布的氣候因子主要是降水季節性變異系數、氣溫年變化幅度、年降水量、海拔和土壤碎石體積百分比。其中氣候因子占主導作用,地形因子次之。適宜岷江冷杉分布的海拔范圍2310 —3757 m,氣溫年變化幅度的范圍為29.3—32.5 ℃;降水季節性變異系數的適宜范圍為76.5—85.4,以及年降水量的適宜范圍為694.1—791.7 mm。

(3)在未來氣候變化情景下,岷江冷杉的適生分布區南移,新增適生區多分布于盆周山地,而喪失適生區則主要分布于原適生區北緣和西北緣。其中在岷江冷杉潛在適生區在未來兩個時期不同氣候變化情景下,喪失面積最大的是2070s的RCP8.5排放情景,喪失面積最小的是在2070s的RCP4.5排放情景。

猜你喜歡
物種模型
物種大偵探
物種大偵探
一半模型
吃光入侵物種真的是解決之道嗎?
英語世界(2023年10期)2023-11-17 09:18:18
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
回首2018,這些新物種值得關注
電咖再造新物種
汽車觀察(2018年10期)2018-11-06 07:05:26
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲美女一级毛片| 中文字幕无线码一区| 国产一区二区精品高清在线观看| 免费观看成人久久网免费观看| 国产剧情伊人| 国内精品久久人妻无码大片高| 激情五月婷婷综合网| 国产成人精品一区二区| 国产精品短篇二区| 国产亚洲现在一区二区中文| 久久精品中文字幕免费| 亚洲国产无码有码| 欧美成人精品高清在线下载| 国产肉感大码AV无码| 1769国产精品视频免费观看| 亚洲色无码专线精品观看| 亚洲制服丝袜第一页| 久夜色精品国产噜噜| 国产又色又爽又黄| 欧美成人aⅴ| 国产精品视频第一专区| 亚洲第一天堂无码专区| 日韩不卡高清视频| 国内精品自在自线视频香蕉| 精品欧美视频| 日韩在线2020专区| 成年av福利永久免费观看| 欧美特级AAAAAA视频免费观看| 国产精欧美一区二区三区| 综合亚洲色图| 久久男人资源站| 在线视频一区二区三区不卡| 国产精品视频导航| 毛片手机在线看| 国产av无码日韩av无码网站| 亚洲一区二区无码视频| 国产精品无码作爱| 久久精品一品道久久精品| 日韩视频精品在线| 欧美精品影院| 91精品在线视频观看| 亚洲综合精品香蕉久久网| 国产精品自在线拍国产电影| 国产69囗曝护士吞精在线视频| 国产高清无码第一十页在线观看| 久久不卡精品| 99久久精品免费看国产免费软件| 成人免费视频一区| 女人18毛片水真多国产| 91精品国产91久无码网站| 91区国产福利在线观看午夜| 97精品久久久大香线焦| 大香网伊人久久综合网2020| 亚洲成a人片77777在线播放| 2018日日摸夜夜添狠狠躁| 91久久青青草原精品国产| 97视频在线精品国自产拍| 88av在线播放| 国产视频a| 国产97视频在线| 欧美爱爱网| 国产高清不卡视频| 精品自拍视频在线观看| 精品国产成人av免费| 亚洲综合九九| 69av免费视频| 在线看免费无码av天堂的| 国产永久无码观看在线| 天天综合网站| 日韩在线播放欧美字幕| 国产在线第二页| 欧美日韩91| 青青草一区二区免费精品| 伊人网址在线| 久久精品女人天堂aaa| 亚洲中文制服丝袜欧美精品| 午夜a视频| 日韩成人高清无码| 久久综合色视频| 午夜高清国产拍精品| 国产高潮视频在线观看| 国产精品性|