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

不同氣候條件下沙冬青屬植物在我國的潛在分布
——基于生態位模型預測

2020-11-24 09:08:50段義忠王海濤杜忠毓賀一鳴柴乖強
生態學報 2020年21期
關鍵詞:物種植物模型

段義忠,王 馳,王海濤,杜忠毓,2,賀一鳴,柴乖強

1 榆林學院生命科學學院,陜西省陜北生態修復重點實驗室, 榆林 719000 2 寧夏大學西北退化生態系統恢復與重建教育部重點實驗室, 銀川 750021

西北干旱區荒漠植物的形成與古地中海退卻、喜馬拉雅及青藏高原的隆升、荒漠的形成演化及第四紀冰期的氣候變化等密切相關。第四紀冰期和間冰期反復交替引起的氣候波動對該區域植物地理分布格局產生了重大影響[1]。地質、氣候與生境的巨大變化對西北干旱區植物區系和生物多樣性產生了重要影響,更新世以來冰期與間冰期反復作用,進一步影響到該區生物類群的分化、遷移和擴散。西北干旱區植物區系種類相對貧乏,但其廣闊的面積、特殊的干旱地理和歷史成分,使其在中國植物區系乃至歐亞及北溫帶植物區系中都占有重要地位。近年來,全球平均溫度的上升,對植被的生長及其分布產生了巨大的影響,使全球分布的多個物種面臨滅絕的危險[2- 4]。未來氣溫的持續變化及降水格局改變,諸多生態系統和生物類群將受到氣候變化的影響[5]。因此研究氣候變化對該區域物種潛在分布、物種保護及資源利用等具有十分重要的意義[6]。

生態位模型是根據物種實際分布范圍和環境變量,通過一定的算法來預測物種能否在該區域生存的模型[7],并能預測植物在全球未來氣候變化下的潛在分布[8],越來越多的生態位模型被應用于物種歷史地理分布過程及未來氣候變化下分布趨勢的研究[9]。目前應用比較廣泛的生態位模型主要有最大熵值模型(MaxEnt),基于生物氣候數據的Bioclim和Domain模型。MaxEnt模型是通過物種“只存在”的分布數據信息以最大熵原理為基礎,統計已知研究物種的分布點,根據物種現實的分布信息和環境變量,推斷物種未知概率分布的數學計算方法,再得到目標物種的潛在分布[10]。Bioclim模型是一種框架生態位模型[11],它是通過獲取已有的物種分布數據和環境數據,再利用分布百分比算法對物種在一段時間內可能的分布區進行分析,計算每個網格的環境適生度[12]。Domain模型是利用Gower算法,通過計算環境空間中的點和點之間的最大相似程度,最終確定物種的分布范圍[13]。近年來,生態位模型對于物種生境的適宜性判斷運用較為廣泛,眾多學者已經通過生態位模型對多個物種的潛在適生區進行模擬,得到了很好的預測效果[14- 16]。

沙冬青屬(Ammopiptanthus)植物,隸屬于豆科(Leguminosae),是古地中海第三紀孑遺珍稀植物、亞洲荒漠區特有的瀕危植物、也是我國西北干旱地區唯一的超旱生常綠闊葉灌木[17- 19]。沙冬青屬植物具有防風固沙、改良土壤、水土保持、耐干旱、耐嚴寒、耐鹽堿化、耐貧瘠、荒漠防沙治沙、藥用和根系發達等特性,而備受研究人員關注[20- 23]。該屬包含兩個種,沙冬青(Ammopiptanthusmgolicus)和矮沙冬青(Ammopiptanthusnanus)。沙冬青又名蒙古沙冬青,主要分布在中國內蒙古阿拉善左旗、甘肅武威市、寧夏等地;矮沙冬青又名新疆沙冬青,主要分布在中國新疆喀什、烏恰縣等地區,在天山西部(相鄰的吉爾吉斯斯坦)也有分布[24- 26]。盡管國內外針對沙冬青屬植物在分子生物學、環境生態學、遺傳分子學、區系分布、瀕危原因、抗旱抗寒機理和利用價值等基礎部分進行了深入細致的研究[27-33],但目前未發現基于生態位模型對其在我國潛在分布區預測的相關研究報導。基于此,本研究基于MaxEnt、Bioclim和Domain 3種模型,結合ArcGIS軟件對沙冬青屬自然種群在末次間冰期(Last Interglacial)、末次冰盛期(Last Glacial Maximum)、當代(Current)和未來(2050s)四個不同時期氣候情景下在我國的潛在分布區進行預測,分析其與環境變量之間的聯系,揭示影響其潛在分布的主導因子,并為該屬植物今后的繁殖、資源利用及保護等提供一定的理論參考依據,同時本研究將對全球氣候變化下干旱半干旱地區物種多樣性和生態系統穩定性具有重大意義。

1 材料和方法

1.1 物種分布數據的獲取

沙冬青和矮沙冬青地理分布數據收集自中國數字植物標本館(CVH, http: //www.cvh.org.cn/)、《FloraofChina》、野外調查和公開發表的文獻。為保證其樣點信息的準確無誤,對所選取樣點進行篩選,最終得到沙冬青68個分布點,矮沙冬青57個分布點(圖1)。

圖1 沙冬青屬分布點Fig.1 The distribution points of Ammopiptanthus

1.2 地理和環境數據

地理數據在國家基礎地理信息系統(http: //nfgis.nsdi.gov.cn/)中下載1:400萬的中國行政區劃矢量圖作為分析的底圖。末次間冰期(Last Interglacial,130 ka BP)、末次冰盛期(Last Glacial Maximum,21 ka BP)、當代(1961—2000 年)、2050年(2041—2060年)環境數據從世界氣候數據庫(http: //www.worldclim.org/)下載,每個時期包括19個空間分辨率為2.5′(1km2)的環境氣候數據[34](表1)。為模擬未來2050年(2050s)氣候變化情景下物種的分布,下載CCSM4全球氣候模式中典型濃度路徑4.5(RCP 4.5)和典型濃度路徑8.5(RCP 8.5)兩種溫室氣體排放情景[35]的氣候數據,環境數據分辨率均為2.5′[36]。

表1 研究使用的環境數據

1.3 模型預測及評價

MaxEnt模型在MaxEnt 3.3.3k軟件中進行,分別導入沙冬青與矮沙冬青的地理分布數據和19個環境變量,選擇刀切法(Jackknife)計算變量的貢獻率和置換重要性,隨機選取25%的分布數據作為測試樣本,75%作為訓練樣本,設置10次重復,其他參數默認[37]。最終利用ArcGIS 10.0軟件進行作圖。Bioclim和Domain模型在DIVA-GIS軟件中運行。將10組訓練集和驗證集文件導入DIVA-GIS軟件,選擇CLASS_REP字段的全部類別及最大空間范圍[38],再將結果中AUC值最大的“.grd”文件,轉換為“.asc”柵格文件格式,導入ArcGIS 10.0軟件中進行可視化。基于MaxEnt、Domain和Bioclim模型得到沙冬青屬在我國分布區的柵格文件,利用自然斷點法進行重分類劃分。本研究中采用非閾值依賴判斷方法類中的受試者工作特征(receiver operating characteristic,ROC)曲線來評價預測效果[39- 40]。ROC評價曲線下的面積為AUC值,AUC值可用于判斷預測結果的準確性,取值范圍為0—1,AUC值小于0.6時認為模型預測結果失敗,0.6—0.7時,模型預測結果較差,0.7—0.8時,模型的預測效果一般,0.8—0.9時,其預測結果良好,當為0.9—1.0時,模型預測效果優秀。AUC值越接近1時,模型預測結果越可靠,反應物種的潛在分布能力越準確[41]。本研究劃分為非適生區、邊緣適生區、低適生區、中適生區、高適生區和最佳適生區6個等級。

2 結果與分析

2.1 模型的準確性評價及其預測沙冬青屬當代的分布區

2.1.1MaxEnt、Bioclim和Domain模型的準確性評價

基于物種分布數據和環境變量在MaxEnt、Bioclim和Domain模型中對沙冬青和矮沙冬青當代的潛在地理分布范圍進行了模擬,10次重復平均AUC值均大于0.80(圖2),沙冬青分別為0.991、0.925和0.967,矮沙冬青分別為0.991、0.829和0.957,綜合發現MaxEnt模型在預測沙冬青屬潛在地理分布區中表現最佳。此外,末次冰盛期、末次間冰期及未來氣候變化下(2050s)的模擬中,AUC值均高于0.80,表明預測結果良好準確。

圖2 MaxEnt、Bioclim和Domain模型對沙冬青和矮沙冬青潛在適生區預測結果的AUC值Fig.2 The AUC value of MaxEnt, Bioclim and Domain models for the prediction results of potential suitable area of A. mgolicus and A. nunas

2.1.2MaxEnt、Bioclim和Domain模型預測沙冬青屬當代的潛在分布區

利用MaxEnt、Bioclim和Domain模型對沙冬青屬在我國潛在分布區進行預測,發現MaxEnt和Bioclim模型預測結果與沙冬青屬當代在我國的實際分布數據最接近,而Domain模型偏差極大(圖3),表明Domain模型并不適用于沙冬青屬在我國的潛在分布。基于10次重復的平均AUC值,本研究中舍棄Domain模型,選擇效果較好的MaxEnt及Bioclim模型對沙冬青和矮沙冬青在我國潛在適生分布區進行預測。

圖3 沙冬青和矮沙冬青當代分布點和MaxEnt、Bioclim和Domain模型預測沙冬青和矮沙冬青當代的潛在分布區Fig.3 Distribution points of A.mgolicus and A. nunas and the predict potential distribution area by MaxEnt, Bioclim and Domain

2.2 沙冬青屬在我國的潛在適生分布區

MaxEnt模型預測結果表明,沙冬青在末次間冰期分布區面積占全國總面積的5.51%,主要分布在內蒙古鄂爾多斯市和阿拉善左旗、甘肅武威、寧夏銀川、陜西榆林、新疆吐魯番和烏魯木齊等區域。高適生區和最佳適生區主要分布在內蒙古鄂爾多斯、寧夏銀川和甘肅武威等區域,占全國總面積的1.15%;低適生區和中適生區主要分布在內蒙古鄂爾多斯市、阿拉善左旗、甘肅民樂縣、新疆烏魯木齊市以及河北張家口,占全國總面積的1.81%;邊緣適生區主要分布在陜西榆林、甘肅武威、山西太原、內蒙古鄂爾多斯和河北石家莊等區域。末次冰盛期與末次間冰期相比,沙冬青適生分布區面積增加0.96%,占全國總面積的6.47%,增加部分主要集中在邊緣適生分布區,但最佳適生區和高適生區面積減少0.02%,主要是陜西榆林市。沙冬青在我國當代時期的潛在分布區總面積與末次冰盛期相比減少2.46%,但其最佳適生區和高適生區分布面積增加1.65%,主要集中在內蒙古鄂爾多斯、阿拉善、寧夏石嘴山等區域。邊緣適生區分布面積減少2.46%,低適生區和中適生區分布面積減少為1.11%,占全國總面積的0.68%(圖4)。

圖4 MaxEnt模型預測的不同時期沙冬青在我國的潛在適生分布區Fig.4 Potential suitable area of A. mgolicus in China predicted by MaxEnt model in different period

與當代時期沙冬青在我國的適生分布區相比,其未來(2050s)在我國潛在適生分布區面積均有小幅度的增加。RCP4.5氣候情景下,沙冬青總潛在適生分布區面積占全國面積的4.03%,較當代適生分布面積僅增加0.02%;其中最佳適生區及高適生區面積占全國總面積的2.90%,低適生區及中適生區面積占全國總面積的0.56%,邊緣適生區面積占全國總面積的0.57%。RCP8.5氣候情景模式下,沙冬青總潛在適生分布區面積占全國總面積的4.15%,較當代時期總潛在適生分布面積增加0.14%。最佳及高適生分布區面積占全國總面積的2.88%,低適生區及中適生區面積占全國總面積0.67%,邊緣適生區面積占全國總面積0.60%。兩種氣候情景下沙冬青在我國潛在適生分布區域均以寧夏賀蘭山為中心向內蒙古西部地區和甘肅東部地區擴散(圖4)。

基于MaxEnt模型預測矮沙冬青在我國的潛在適生分布區發現,隨時間變化有擴張趨勢(圖5)。其在末次間冰期適生分布區占全國總面積的0.78%,高適及最佳適生區主要集中分布在新疆中西部烏恰縣及喀什地區,占我國總面積的0.30%。中適及低適生分布區僅僅占我國總面積的0.18%,零散分布在我國西藏自治區、青海省,邊緣適生分布于我國的甘肅省、陜西省和山西省。與末次間冰期相比,末次冰盛期總分布面積增加了0.39%。主要增加新疆塔城、西藏阿里及青海共和縣等中低適生分布區,有向我國東部遷移趨勢,如陜西省、山西省、河北省等地出現邊緣適生區。矮沙冬青當代在我國潛在適生分布區占全國總面積的2.91%,較末次冰盛期增加1.79%。主要是高適及最佳適生區有所擴張,占全國總面積的2.23%,主要分布在新疆烏恰縣、甘肅張掖、青海格爾木一帶。分布區范圍多集中在36°—44°N,75°—101°E。

未來(2050s)不同氣候情景下,矮沙冬青的潛在適生分布區域有所不同。RCP 4.5情景下矮沙冬青潛在適生區域占全國總面積的3.01%,與當代適生分布區相比總適生區面積增加0.05%,高適生區主要集中在新疆烏恰縣南部、烏魯木齊北部及寧夏等周邊地區。但在RCP 8.5情景下矮沙冬青適生分布區占全國總面積的2.86%,與當代適生分布區相比總適生區面積減少0.10%(圖5)。未來(2050s)兩種氣候情景下矮沙冬青潛在適生分布區均有向北遷移趨勢。

圖5 MaxEnt模型預測的不同時期矮沙冬青在我國的潛在適生分布區Fig.5 Potential suitable area of A. nunas in China predicted by MaxEnt model in different period

Bioclim模型預測沙冬青末次間冰期在我國的潛在分布區占全國總面積的4.48%(圖6)。高適及最佳適生區面積約占全國總面積的0.58%。與末次間冰期相比,沙冬青在我國的潛在適生區面積增加4.75%,主要分布在我國西北地區,占全國總面積的9.23%。沙冬青當代在我國的潛在分布區主要分布在以內蒙古、甘肅、寧夏、陜西交界地區,占全國總面積的8.30%。未來RCP4.5情景下,沙冬青在我國潛在適生分布區總面積占我國總面積的8.08%,與當代相比減少0.22%。未來RCP8.5情景下,沙冬青在我國潛在適生分布區總面積占我國總面積的7.63%,與沙冬青當代在我國潛在適生面積相比減少0.67%。兩種不同情景下,沙冬青在我國潛在適生分布區面積均有所減少,未來RCP8.5情景下沙冬青在我國潛在適生分布區縮減面積更大(圖6)。

圖6 Bioclim模型預測的不同時期沙冬青在我國的潛在適生分布區Fig.6 Potential suitable area of A. mgolicus in China predicted by Bioclim model in different period

與MaxEnt模型預測矮沙冬青在我國潛在分布區相似,Bioclim模型預測矮沙冬青末次間冰期在我國潛在適生分布區面積占全國總面積的2.19%,高適及最佳適生區面積僅為全國總面積的0.04%。與末次間冰期相比,矮沙冬青末次冰盛期適生分布區面積增加0.65%,占全國總面積的2.84%,主要分布在我國新疆地區,總體分布區向我國北部遷移。矮沙冬青當代在我國潛在適生分布區占全國總面積的4.39%,主要分布在我國新疆烏恰縣。未來RCP4.5情景下,矮沙冬青在我國潛在適生區分布面積占全國總面積的3.42%。與當代潛在適生分布區面積相比,矮沙冬青在我國潛在適生區分布面積減少0.97%。未來RCP8.5情景下,矮沙冬青在我國潛在適生區分布面積占全國總面積的4.29%(圖7)。

圖7 Bioclim模型預測的不同時期矮沙冬青在我國的潛在適生分布區Fig.7 Potential suitable area of A. nunas in China predicted by Bioclim model in different period

2.3 影響沙冬青屬潛在分布的主導環境變量

利用Jackknife刀切法分析各環境變量對沙冬青屬分布適宜度(圖8),結合各環境變量對沙冬青屬分布影響的貢獻率(表2)表明:溫度是影響沙冬青屬分布的主要因素。具體而言,年均溫(Bio1)的貢獻率14.4%、溫度年較差(Bio7)貢獻率10.5%、最冷月的最低溫(Bio6)貢獻率9.7%、最濕季度均溫(Bio8)貢獻率9.1%、降水量季節性變化(Bio15)貢獻率8.7%以及最冷季度的平均降水量(Bio19)貢獻率8.5%,其累積貢獻率達61.5%。說明上述6個環境變量是影響沙冬青適宜分布生境的主導環境變量。

表2 基于MaxEnt模型預測的影響沙冬青和矮沙冬青分布的環境變量的貢獻率

圖8 基于刀切法的潛在環境變量對沙冬青和矮沙冬青分布的貢獻率Fig.8 Contribution rate of potential environmental variables based on Jackknife to the A. mgolicus and A. nanus Bio1: 年均溫Annual mean temperature; Bio2: 平均日較差Mean diurnal range; Bio3: 等溫性Isothermality; Bio4: 溫度季節性變化Variation of temperature seasonality; Bio5: 最熱月最高溫The maximum temperature of the warmest month; Bio6: 最冷月最低溫The minimum temperature of the coldest month; Bio7: 溫度年較差Temperature annual range; Bio8: 最濕季度均溫Mean temperature of the wettest quarter; Bio9: 最干季度均溫Mean temperature of the driest quarter; Bio 10: 最熱季度均溫Mean temperature of the warmest quarter; Bio11: 最冷季度均溫Mean temperature of the coldest quarter; Bio12: 年降水量Annual precipitation; Bio13: 最濕月降水量Precipitation of the wettest month; Bio14: 最干月降水量Precipitation of the driest month; Bio15: 降水量季節性變化Variation of precipitation seasonality; Bio16: 最濕季度降水量Precipitation of the wettest quarter; Bio17: 最干季度降水量Precipitation of the driest quarter; Bio18: 最熱季度降水量Precipitation of the warmest quarter; Bio19: 最冷季度降水量Precipitation of the coldest quarter

影響矮沙冬青適宜分布生境的關鍵環境變量有溫度季節性變化(Bio4)、等溫性(Bio3)、最冷季度的平均降水量(Bio19)、最濕月降水量(Bio13)、平均日較差(Bio2)和溫度年較差(Bio7),貢獻率分別為17.9%、16.8%、9.8%、6.3%、5.2%和5.2%,累積貢獻率達61.2%。說明矮沙冬青適宜的分布生境主要由上述6個環境變量控制。

3 討論

隨著全球物種分布數據的共享和快速發展的空間分析技術,生態位模型在生物多樣性保護的多個領域得到開拓和應用,其中相關性方案的生態位模型利用物種已知的分布數據和相關環境變量,根據一定的算法運算來構建模型,判斷物種的生態需求,并將運算結果投射至不同的時間和空間中以預測物種的實際和潛在分布,被越來越多地應用在入侵生物學、保護生物學、全球氣候變化對物種分布的影響,生物地理學及傳染病空間傳播的研究中。例如在入侵生態學研究上,岳茂峰等[42]利用MaxEnt模型對刺軸含羞草入侵植物在我國的潛在分布進行了模擬,得出刺軸含羞草在我國的適生區主要分布在云南、海南、廣東西南部以及中國臺灣地區,到2050年,在溫室氣體A1排放模式下其在全球適生區面積與當代相似,但在我國的適生區略有減少;在保護生物學研究上,張文秀等[43]運用生態位模型對5個時期瀕危植物白豆杉的潛在分布區進行重建,表明在未來氣候變暖條件下,白豆杉適宜分布區減少,該物種往高海拔山地收縮,低海拔(約600 m以下)的適宜區基本消失。

本研究基于三種生態位模型,MaxEnt、Bioclim和Domain在空間尺度上和時間尺度上對沙冬青屬植物在我國的潛在適生區進行預測。結果表明,MaxEnt模型、Bioclim模型和Domain模型相比較,其MaxEnt模型對沙冬青屬的分布區預測較為準確,更有利于對該物種的長期保護和研究。MaxEnt模型預測出沙冬青最適生區主要分布在我國內蒙古鄂爾多斯市、阿拉善左旗、甘肅中東部及寧夏北部等干旱地區;矮沙冬青最適生區主要分布在我國新疆烏恰縣和烏魯木齊地區,這與目前對沙冬青屬植物的研究分布保持一致[44]。MaxEnt模型預測沙冬青屬的潛在分布區可用于對該物種的核心分布區進行研究,這對于沙冬青和矮沙冬青兩種瀕危植物的資源保護及其推廣種植具有重要的意義。

氣候環境是物種分布的決定性因素,而物種的分布變化也是氣候變化的最直接,最明顯的反應。全球變暖已經改變了陸地生態系統的結構,這反過來又改變了物種的棲息地和地理分布函數[45]。如氣候、地形、土壤、人為干擾和空間限制因素對不同空間尺度的分布具有重要意義[1]。預測結果表明沙冬青和矮沙冬青的空間分布差異明顯,并且呈現不連續分布,與先前的研究結果保持一致[46]。影響沙冬青潛在地理分布的最主要環境變量是年均溫,影響矮沙冬青潛在分布的最主要環境變量是溫度季節性變化,如在最濕月份和最干月份等,這可能是因為沙冬青所處海拔要低于矮沙冬青所處海拔,且沙冬青地理環境處于荒漠帶與草原帶交匯地,也是多種植物區系交匯的地區,且由于其位于荒漠帶東部的水分條件較好,造成了兩者之間的差異[47]。沙冬青屬植物的分布范圍會受氣候變化的影響,隨著時間推移,該屬的兩個物種分布會進行不同程度的減少或擴張[48]。沙冬青屬植物在2050年(RCP4.5和RCP8.5)的最佳適生區及高適生分布相比,在高等CO2濃度下(RCP8.5)最佳適生區及高適生范圍有所縮減,呈現破碎化分布,這可能受我國氣候整體呈現變暖、變濕的影響[49];也進一步說明了沙冬青屬植物未來分布是受到其主導環境變量水分及溫度的影響。

在調查中我們發現,僅利用環境變量預測沙冬青屬植物的潛在分布還存在一定的誤差,實際上在我國西北干旱區人為活動、土壤養分等也會影響沙冬青屬植物的生長與分布;目前在人煙稀少的沙地和黃河灘前平原分布較多,且長勢較好,群落更新較快;然而在靠近礦區、鐵路及放牧區等被破壞較為嚴重,僅剩年齡較大的沙冬青屬植物存在,這說明人為破壞因素成為沙冬青屬植物分布區日漸縮減的主要因素之一。通過3種模型預測,沙冬青屬植物目前的分布主要在我國西北干旱半干旱地區,其中沙冬青在我國內蒙古西部及寧夏北部地區分布最為集中,矮沙冬青僅在我國新疆喀什地區南部昆侖山與帕米爾交界的狹長地帶分布集中[50],研究結果對沙冬青屬自然保護區建立及其應用治理環境問題具有十分重要的意義。

4 結論

研究結果表明當代沙冬青最適生區主要集中在內蒙古中部、寧夏和甘肅等地;未來沙冬青最適生區在內蒙古中部、寧夏北部和甘肅北部呈現向外擴張的趨勢。矮沙冬青最適生區主要集中在新疆南部;未來矮沙冬青最適生區向新疆烏恰縣南部、烏魯木齊北部移動和擴大。未來氣候環境變化將極大的影響植物的分布,在2050年CO2濃度為RCP4.5和RCP8.5情景下,沙冬青屬植物的潛在適生區面積都有所增加,但在高濃度CO2(RCP8.5)情景下,沙冬青屬植物在最佳及高適生區將減少,建議對沙冬青屬植物最佳及高適生區建立自然保護區。影響沙冬青分布的主要環境變量是年均溫、最冷月最低溫和最冷季度的平均降水量;影響矮沙冬青分布的主要環境變量是最冷季度的平均降水量、最濕月降水量和溫度季節性變化。

猜你喜歡
物種植物模型
一半模型
吃光入侵物種真的是解決之道嗎?
英語世界(2023年10期)2023-11-17 09:18:18
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
回首2018,這些新物種值得關注
電咖再造新物種
汽車觀察(2018年10期)2018-11-06 07:05:26
哦,不怕,不怕
將植物穿身上
3D打印中的模型分割與打包
植物罷工啦?
主站蜘蛛池模板: 国产日韩AV高潮在线| 国产精品极品美女自在线| 亚洲日本韩在线观看| 91综合色区亚洲熟妇p| 尤物午夜福利视频| 特级精品毛片免费观看| 欧美亚洲另类在线观看| 噜噜噜综合亚洲| 亚洲天堂视频在线免费观看| 无码又爽又刺激的高潮视频| 99re视频在线| 四虎成人精品| 99re在线视频观看| 色婷婷在线影院| 免费人成视频在线观看网站| 天堂av综合网| 99精品影院| 亚洲资源站av无码网址| 一级毛片无毒不卡直接观看| 亚洲最大在线观看| AV在线麻免费观看网站| 国产xxxxx免费视频| 人妻出轨无码中文一区二区| 狠狠干欧美| 波多野结衣一区二区三区四区| 国内精自线i品一区202| 精品久久人人爽人人玩人人妻| 自慰网址在线观看| 熟女视频91| 精品三级网站| 国产精品爆乳99久久| 午夜视频免费一区二区在线看| 美女无遮挡免费视频网站| 国产成人免费手机在线观看视频| 精品小视频在线观看| 成人午夜在线播放| 亚洲天堂首页| 亚洲天堂精品在线| 国产三级国产精品国产普男人 | 手机在线国产精品| 久久久无码人妻精品无码| 九九这里只有精品视频| 一级毛片在线播放免费| 久久精品欧美一区二区| 人妻中文久热无码丝袜| 免费看一级毛片波多结衣| 福利在线不卡| 国产精品99在线观看| 88av在线| 亚洲中文字幕av无码区| 欧美日韩北条麻妃一区二区| 国产精品亚洲αv天堂无码| 欧美精品高清| 亚洲欧洲天堂色AV| 精品久久久久久久久久久| 久久国产精品影院| 欧美成人一级| 国产精品亚洲专区一区| 91色在线观看| 99人妻碰碰碰久久久久禁片| 精品视频福利| 亚洲色欲色欲www在线观看| 无码AV高清毛片中国一级毛片 | 91精品啪在线观看国产60岁| 午夜一级做a爰片久久毛片| 欧美激情第一区| 拍国产真实乱人偷精品| 91在线免费公开视频| 色亚洲成人| 热re99久久精品国99热| A级毛片无码久久精品免费| 欧美不卡视频在线观看| 成人va亚洲va欧美天堂| 中文字幕波多野不卡一区| 日韩国产无码一区| 免费国产不卡午夜福在线观看| 久久性视频| 全午夜免费一级毛片| 亚洲精品国产综合99久久夜夜嗨| 国产精品无码久久久久AV| 亚洲欧洲日产无码AV| 97综合久久|