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

人為引入和氣候變化對灰喜鵲未來分布的影響

2024-01-17 02:52:36顏尉珂張新宇
生態學報 2023年24期
關鍵詞:物種模型

顏尉珂,雷 宇,王 磊,張新宇,劉 強,*

1 西南林業大學云南省高原濕地保護修復與生態服務重點實驗室,昆明 650224 2 西南林業大學國家高原濕地研究中心/濕地學院,昆明 650224

氣候是影響物種分布的重要驅動因素[1]。聯合國政府間氣候變化專門委員會(IPCC)第六次評估報告指出,在當前的溫室氣體排放情景下,到2100年全球溫升可達到2.2-3.5℃[2]。由于全球氣候變化的區域變異性以及物種本身的生態位特征,不同的物種會表現出不同的響應模式。目前,已有許多物種表現出向高緯度、高海拔地區擴散或轉移的現象[3-5]。但某些物種也會呈現出不同的變化趨勢,如我國新疆地區的鳥類向東擴張[6],以及部分物種的分布區萎縮[5, 7-8]。此外,由于寵物貿易、園林建設等人為引入的原因,也導致了外來物種的定殖和擴張[9-10]。

灰喜鵲(Cyanopicacyanus)隸屬于雀形目(PASSERIFORMES)鴉科(Corvidae)灰喜鵲屬,主要分布于東亞,其中我國北方是其主要分布區[11]。主要以鱗翅目(LEPIDOPTERA)幼蟲和直翅目(ORTHOPTERA)昆蟲為食,是重要的食蟲鳥類,在國內常被人為引種用以防治松毛蟲(Dendrolimus)等森林害蟲[9, 12-14],對保持森林生態系統的穩定性具有重要作用。由于人為引種、籠養逃逸及放生等原因,已有一定數量的灰喜鵲在其自然分布地以外的地方建立種群。1975年,灰喜鵲作為逃逸物種在香港被記錄[15];廣東廣州、廣西桂林和福建光澤分別在1985年、1986年和1992年引進灰喜鵲以防治森林蟲害[12-13, 16];1988年,昆明動物園從北京引進灰喜鵲以供游客觀賞,后于1991年釋放野外[17];2011年,王龍舞等人在海南海口發現本地籠養逃逸的灰喜鵲種群[15]。目前,灰喜鵲已在甘肅[18]、云南[17]、廣西[16]、廣東[13]、福建[12]、海南[15]、香港[15]及臺灣[19]等省份定殖并呈現出擴散的趨勢。

物種分布模型(Species Distribution Model,SDM)也常稱為生態位模型(Ecological Niche Model,ENM),是一種結合了物種分布數據和相關環境變量數據的數學模型[20-21],此類模型基于一定的算法,在給定的環境變量構成的生態空間內,根據分布數據及其所關聯的環境變量參數來計算物種對生態位的需求,并將此結果投射到給定的時間和空間范圍內,最后以概率的方式呈現物種對不同時空中生境的喜好程度[21-23]。近年來,物種分布模型被廣泛地應用于保護生物學、入侵生物學及氣候變化對物種分布影響的研究中[21-22, 24]。最大熵(Maximum Entropy,MaxEnt)模型是目前應用最為廣泛、性能可靠的物種分布模型之一[24-27]。

本研究基于氣候生態位理論,結合灰喜鵲在我國境內的分布數據及氣候變量,通過構建MaxEnt模型,旨在探討以下科學問題:(1)灰喜鵲自然分布地種群和引入地種群對氣候因子的響應情況;(2)在人為引入與氣候變化背景下,當前和未來時期灰喜鵲種群在我國境內的分布格局與變化趨勢。這將有助于理解全球氣候變化背景下物種的分布及變化規律,并預測灰喜鵲對生態系統的潛在影響。

1 數據來源與研究方法

1.1 分布數據獲取及處理

分布數據的來源包括全球生物多樣性信息數據庫(https://www.gbif.org/)、中國觀鳥記錄中心(http://www.birdreport.cn/)以及文獻記錄[12-13, 15]。共獲取到61026條分布數據,分布數據收集時間截止2022年7月。

首先刪除信息錯誤、缺失、模糊及重復的分布數據,然后將剩余數據在ArcMap 10.6中與使用的環境變量疊加,刪除在環境變量范圍外的分布數據。為了降低空間自相關對模型構建的干擾作用,使用R 4.1.3中的Wallace包對分布數據進行篩選,每個2.5 min×2.5 min的柵格里僅保留1個分布點[28-29]。最終保留有效分布點2458個。

以《中國動物圖譜》[30]中灰喜鵲分布界線為參照,將位于該區域內的分布點定義為自然分布點,共計2203條。將位于自然分布地外且有明確的證據表明灰喜鵲為引入物種的地點定義為引入地,最終確定新疆烏魯木齊、云南昆明、貴州貴陽、廣西桂林、廣東廣州、福建光澤、海南、香港和臺灣等9個省市為引入地[12-13, 15-17, 19],篩選出有效分布點129個(圖1)。

圖1 灰喜鵲分布區和分布點Fig.1 Distribution area and points of Azure-winged Magpie in China

1.2 環境變量獲取及篩選

1.2.1環境變量獲取

本研究共選取世界氣候數據庫(https://www.worldclim.org/)中的19種氣候變量用于模型分析(表1),分辨率為2.5 minutes。未來時期的環境變量選擇2041-2060年(2050s時期)和2061-2080年(2070s時期)兩個階段,氣候模式選擇北京氣候中心發布的BCC-CSM2-MR氣候系統模式,該模式對于中國的氣候變化趨勢具有較好的模擬效果[31],情景模式選擇ssp245,即中等輻射強迫與中等社會脆弱性的組合情景[32]。

表1 環境變量列表Table 1 List of environment variables

1.2.2環境變量篩選

環境變量在物種分布模型中的過度擬合問題受到了廣泛的關注,為了避免模型的過度擬合,剔除部分環境變量可以使模型的預測性能得到一定的提升[33-34]。學者們普遍采用Pearson相關系數分析法對模型使用的環境變量進行多重線性相關性分析,以環境變量間的相關性系數結合環境變量的貢獻率作為篩選的依據[22, 35]。使用ArcMap 10.6中的“Band Collection Statistics”工具對使用的19種環境變量進行相關性分析,相關性系數r表示兩兩變量之間線性關系的密切程度,本研究選擇|r|>0.8作為環境變量間有強相關性的標準[35-36];然后將19種環境變量導入到MaxEnt模型中試運行,使用模型結果中環境變量的貢獻值作為篩選的依據[27, 33, 37],環境變量的貢獻率表示環境變量的重要程度,貢獻率越高的環境變量對預測結果的貢獻越大[38];最后結合以上兩種方法的結果來篩選環境變量,當兩個環境變量之間的相關性系數r>0.8時,則保留貢獻率更高的環境變量。最終優選出bio3、bio4、bio10、bio11、bio12和bio15等6種環境變量用于綜合分布模型與自然分布地模型,bio2、bio3、bio6、bio7、bio8、bio13、bio14和bio15等8種環境變量用于引入地模型。

1.3 MaxEnt模型構建

1.3.1模型參數優化

使用R 4.1.3中的Wallace包分別對自然分布地模型、引入地模型及綜合分布模型進行計算,優選出每個模型最佳的調控倍頻(Regularization multiplier, RM)和特征組合(Feature combination, FC)。在Wallace中設置調控倍頻的范圍為0.5-4.0,間隔為0.5;特征組合項選擇:H、L、LQ、LQH、LQHP(H,hinge,片段化;L,linear,線性;Q,quadratic,二次型;P,product,乘積型)[29, 39-41]。delta.AICc是一種模型評估指標,在候選模型中delta.AICc值最低的模型被確定為最優模型[42],根據Wallace計算結果中最低的delta.AICc值選擇3個模型最優的RM和FC:自然分布地模型使用RM=0.5,FC=H;引入地模型使用RM=3,FC=LQHP;綜合分布模型使用RM=0.5,FC=LQHP。

MaxEnt模型的其他參數設置如下:選擇25%的分布數據作為測試集,75%的分布數據作為訓練集;設置最大背景點數為10000;模型重復次數為10次;重復運行模式選擇Crossvalidate,即交叉驗證法;最大迭代次數為1000次;選擇10 percentile training presence閾值規則作為適生區與非適生區的劃分依據;其余模型參數為默認參數。

1.3.2適生區劃分

采用10 percentile training presence閾值作為劃分適生區與非適生區的界限,適宜值低于該閾值的地方為非適生區,高于或等于該閾值的地方為適生區[35, 43],其中自然分布地模型的閾值為0.2629,引入地模型的閾值為0.1170,綜合分布模型的閾值為0.2579。然后將適生區按照適宜值由低到高三等分為低適生區、中適生區和高適生區[38]。

1.3.3模型預測精度評估

以AUC值作為模型預測精度評估的標準,值越大表示模型預測結果越有規律,環境變量與物種地理分布之間的相關性越大,模型預測效果越好[26, 38, 44]。當AUC值為0.5-0.6時,代表模型預測失敗;當AUC值為0.6-0.7時,代表模型預測效果較差;當AUC值為0.7-0.8時,代表模型預測效果一般;當AUC值為0.8-0.9時,代表模型預測效果良好;當AUC值為0.9-1.0時,代表模型預測效果優秀[45-46]。

2 結果與分析

2.1 模型評估結果

自然分布地模型、引入地模型和綜合分布模型的AUC值分別為0.851、0.962和0.826,為良好或優秀水平,表明模型預測結果較好。

2.2 環境變量對灰喜鵲分布的影響

對自然分布地種群分布影響較大的環境變量為bio12(年均降雨量)、bio4(溫度季節性變化)和bio11(最冷季度均溫)(圖2)。自然分布地種群較適宜(存在概率>0.5)的年均降雨量的范圍為507.45-760.84 mm和868.06-1345.63 mm,較適宜的溫度季節性變化標準差的范圍為853.59-1138.49,較適宜的最冷季度均溫的范圍為-3.80-6.29℃,自然分布地種群對于3種環境變量的響應曲線均呈現出先增加后減少的“單峰”曲線(圖3)。

圖2 環境變量的貢獻率Fig.2 Contribution rate of environmental variablesbio2:晝夜溫差月均值Mean Diurnal Range;bio3:等溫性Isothermality;bio4:溫度季節性變化Temperature Seasonality;bio6:最冷月份最低溫Min Temperature of Coldest Month;bio7:溫度全年波動范圍Temperature Annual Range;bio8:最濕季度均溫Mean Temperature of Wettest Quarter;bio10:最熱季度均溫Mean Temperature of Warmest Quarter;bio11:最冷季度均溫Mean Temperature of Coldest Quarter;bio12:年均降雨量Annual Precipitation;bio13:最濕月降雨量Precipitation of Wettest Month;bio14:最干月降雨量Precipitation of Driest Month;bio15:降雨量季節性變化Precipitation Seasonality

圖3 自然分布地模型中重要環境變量的響應曲線Fig.3 Response curve of important environmental variables in natural distribution model

對引入地種群分布影響較大的環境變量為bio6(最冷月份最低溫)和bio14(最干月降雨量)(圖2)。引入地種群較適宜的最冷月份最低溫的范圍為7.99-17.30℃,較適宜的最干月降雨量的范圍為16.65-202.32 mm,引入地種群的存在概率均隨著最冷月份最低溫和最干月降雨量的增加而升高(圖4)。

圖4 引入地模型中重要環境變量的響應曲線Fig.4 Response curve of important environmental variables in introduction model

對灰喜鵲種群分布影響較大的環境變量為bio12(年均降雨量)、bio11(最冷季度均溫)和bio4(溫度季節性變化)(圖2),表明灰喜鵲更偏好溫暖潮濕的東部季風區。灰喜鵲種群較適宜的年均降雨量的范圍為509.30-1347.55 mm和2310.55-4176.65 mm,較適宜的最冷季度均溫的范圍為-4.00-6.41℃和15.34-19.62℃,較適宜的溫度季節性變化標準差的范圍為849.26-1142.89,由于綜合分布模型的分布數據主要來自自然分布地,所以貢獻率較高的環境變量與自然分布地模型一致,但也由于引入地分布數據的影響,灰喜鵲種群對這3種環境變量的響應曲線均呈現出了先增加后減少、再增加再減少的“雙峰”曲線(圖5)。

圖5 綜合分布模型中重要環境變量的響應曲線Fig.5 Response curve of important environmental variables in comprehensive distribution model

2.3 潛在適生區的分布與變化趨勢

2.3.1自然分布地模型

通過對比當前和未來時期的適生區面積和分布發現:自然分布地種群適生區的面積和適宜值都將隨著氣候的變化而大幅增加,擴張的區域主要位于高緯度的東北地區和高海拔的青藏高原(圖6,圖7)。當前時期的適生區面積共175.46萬km2,占中國國土總面積的18.24%,主要分布于自然分布地內的華北、華中和華東地區,其中多數為中適生區,高適生區主要分布于湖北、安徽、江蘇及渤海周圍的陸地區域,自然分布地外,僅青藏高原存在較多的適生區,西藏的札達縣還出現了小面積的高適生區;2050s時期的適生區面積共252.24萬km2,占中國國土總面積的26.23%,相比當前時期增加了43.76%,新增的適生區主要位于東北地區和青藏高原。其中高適生區在河北、山東、河南、安徽、江蘇和湖北等地大幅增加,表明灰喜鵲不僅有明顯擴張的現象,自然分布地種群的密度可能也會隨適宜值的升高而增加;2070s時期的適生區面積共283.66萬km2,占中國國土總面積的29.49%,相比2050s時期增加了12.46%,適生區面積小幅增加,分布上的變化較小,僅新增了湖南和川東地區,表現出了輕微的南擴現象。質心分析表明灰喜鵲適生區總體向高緯度方向轉移,相比當前時期,2050s時期的質心朝北偏東23°方向轉移了253.38 km,2070s時期的質心朝北偏東25°方向轉移了229.16 km(圖7)。

圖6 不同時期灰喜鵲適生區面積Fig.6 Suitable areas of Azure-winged Magpie in different periods

圖7 自然分布地模型中不同時期灰喜鵲的適生區分布及質心轉移Fig.7 Distribution of suitable areas and centroid transfer of Azure-winged Magpie in different periods in natural distribution model

2.3.2引入地模型

通過對比當前和未來時期的適生區面積和分布發現:引入地種群適生區的面積和適宜值均隨著氣候的變化而小幅增加,分布上的變化較小(圖6,圖8)。當前時期的適生區面積共103.34萬km2,占中國國土總面積的10.74%,主要分布于自然分布地外的新疆北部、四川東部、云南、廣西、廣東、福建東南部、海南和臺灣,此外在自然分布地內的重慶、湖北、江蘇和上海也存在少量的適生區,其中中、高適生區主要分布于廣西、廣東、海南和臺灣等地,這表明引入地種群很可能在這些沿海地區擴張和定殖;2050s時期的適生區面積共125.82萬km2,占中國國土總面積的13.08%,相比當前時期增加了21.75%;2070s時期灰喜鵲的適生區面積共128.98萬km2,占中國國土總面積的13.41%,相比2050s時期增加了2.5%。質心分析表明灰喜鵲適生區總體向高緯度方向轉移,相比當前時期,2050s時期的質心朝北偏東37°方向轉移了88.85 km,2070s時期的質心朝北偏東46°方向轉移了67.69 km(圖8)。

圖8 引入地模型中不同時期灰喜鵲的適生區分布及質心轉移Fig.8 Distribution of suitable areas and centroid transfer of Azure-winged Magpie in different periods in introduction model

2.3.3綜合分布模型

通過對比當前和未來時期的適生區面積和分布發現:灰喜鵲種群適生區的面積和適宜值都將隨著氣候的變化而大量增加,擴張的區域主要位于高緯度的東北地區和新疆北部、高海拔的青藏高原以及南部的四川東部和兩廣地區(圖6,圖9)。當前時期的適生區面積共226.82萬km2,占中國國土總面積的23.58%,主要分布于自然分布地內的華北、華中和華東地區,與自然分布地模型類似,自然分布地內的低、中、高適生區的分布格局也與自然分布地模型基本一致,此外在自然分布地外的新疆北部、青藏高原、廣西南部、廣東南部、海南和臺灣也存在一定面積的適生區,但適宜程度普遍較低;2050s時期的適生區面積共378.97萬km2,占中國國土總面積的39.4%,相比當前時期增加了67.08%,適生區面積大幅增加,新增的適生區主要位于東北地區、新疆北部、青藏高原、四川東部和兩廣地區,新增的適生區中主要為中、低適生區,當前時期的中適生區大部分轉變為高適生區;2070s時期的適生區面積共431.53萬km2,占中國國土總面積的44.87%,相比2050s時期增加了13.87%,適生區面積和適宜值進一步增加,但分布上的變化較小。質心分析表明灰喜鵲適生區總體向高緯度方向轉移,相比當前時期,2050s時期的質心朝北偏西19°方向轉移了164.07 km,2070s時期的質心朝北偏西34°方向轉移了153.44 km(圖9)。

圖9 綜合分布模型中不同時期灰喜鵲的適生區分布及質心轉移Fig.9 Distribution of suitable areas and centroid transfer of Azure-winged Magpie in different periods in comprehensive distribution model

3 討論

3.1 影響灰喜鵲分布的主要環境變量

氣候因子對物種分布的影響與物種的適應性以及食物的可利用性有關[8, 47-48]。在自然分布地模型與綜合分布模型中,貢獻率最高的環境變量均為年均降雨量。降水主要通過影響植被和食物,從而間接影響鳥類的分布[48-50]。在降水量充沛的地區,鳥類的棲息地和昆蟲等食物資源也相對充足[50]。此外松毛蟲作為灰喜鵲的重要食物來源[12],松毛蟲的卵、幼蟲和蛹需要較濕的環境來保證其正常的生長發育,在干旱的環境中松毛蟲的密度會明顯下降[51]。所以在年均降雨量的響應曲線中,灰喜鵲的存在概率在較低的年均降雨量范圍內會隨著年均降雨量的升高而迅速增加。有研究表明降雨量與鳥類豐度呈正相關[50, 52],本文結果與之一致。在引入地模型中,年均降雨量因貢獻率較低而被最干月降雨量替代,其原因可能是引入地大多為降雨量充足的南方地區,在降雨量有保障的前提下,年均降雨量就不再是影響灰喜鵲分布的主要環境變量,而最干月降雨量則依然是影響引入地種群分布的主要環境變量。除年均降雨量外,最冷季度均溫(或最冷月份最低溫)在3種模型中也表現出了較高的貢獻率。最冷季度均溫過低會影響鳥類的身體機能[47],此外在最冷季度均溫較高的地區,鳥類可利用的食物資源也相對較充足[47, 50]。所以在環境變量響應曲線中,最冷季度均溫(或最冷月份最低溫)較高的地方,灰喜鵲的存在概率也較高。

比較3種模型的環境變量響應曲線的差異,發現引入地種群適宜的最冷季度均溫(或最冷月份最低溫)和年均降雨量比自然分布地種群更高,這和引入地多為溫暖濕潤的南方地區有關。表明灰喜鵲從自然分布地引入到引入地后,氣候生態位產生了重要的變化,證明灰喜鵲具有較強的適應能力,這在胡箭關于灰喜鵲在昆明城區的生態適應的研究中也得以驗證[17]。

3.2 灰喜鵲的分布格局與變化趨勢

本研究結果表明,當前時期灰喜鵲在我國境內的適生區主要分布于其自然分布地內的華北、華中和華東地區。在未來氣候變化的背景下,3種模型的結果均顯示灰喜鵲有顯著的擴張趨勢,主要表現為向高緯度的東北地區和高海拔的青藏高原地區擴散。許多國內外研究也表明,在全球氣候變化條件下,物種有向高緯度、高海拔地區遷移的趨勢[8],該現象在鳥類中也普遍存在[5];Hitch和Leberg研究了北美地區56種鳥類的擴散趨勢,發現分布在低緯度地區的鳥類平均每年向北移動2.35 km[53];Parmesan在對鳥類等物種的研究中發現,這些物種都表現出了向高海拔地區轉移的現象[54];國內的一些研究也證明了我國的鳥類有同樣的擴散趨勢:杜寅等分析了我國鳥類20年間的分布數據,發現共有120種鳥類的分布范圍已向北擴張[48];原分布于印度及東南亞地區的鉗嘴鸛(Anastomusoscitans),近年來已向北擴散至我國[55];李東明等在青海西寧發現了約60只白頭鵯(Pycnonotussinensis),表明白頭鵯不僅存在向北擴散的現象,還存在向高海拔地區擴散的趨勢[56];鮮莉莉在青藏高原東北部對灰喜鵲開展的繁殖研究中表示,灰喜鵲是近年來通過殖民的方式擴散至青藏高原東北部的[18],這也從另一方面驗證了本研究結果的正確性。此外,部分鳥類也會呈現出不同的變化趨勢:馬鳴研究發現由于人為因素、生境的可入侵性、全球氣候變化及生物地理上的空缺等原因,在我國新疆有至少42種鳥類存在向東擴散的現象[6];多數鳥類的分布范圍在逐漸縮小,特別是分布在高緯度、高海拔地區的鳥類,幾乎沒有可供其拓展的生存空間[5, 8]。

我國氣候的變化趨勢與全球氣候變化總體趨勢基本一致[48]。我國大部分地區的氣溫將呈現出增長的趨勢,且高緯度地區的氣溫增幅較大[8, 48-49]。我國的降水量總體上同樣處于增加的趨勢[8],有研究發現青藏高原的降水量正在變得越來越多[57-58],北半球中高緯度的降雨量也處于增長的趨勢[5]。本研究結果發現灰喜鵲更偏好溫暖潮濕的氣候,我國氣候總體上呈暖濕化趨勢,且高緯度、高海拔地區的溫度和降水增加較快,這恰好與本研究結果中灰喜鵲種群的擴散方向一致。所以本研究認為,我國氣候的暖濕化趨勢將導致灰喜鵲向高緯度和高海拔地區擴散以及其種群數量的增加。

由于人為引種、籠養逃逸及放生等原因,目前灰喜鵲已在其自然分布地外建立了多處可自我維持的種群,在自然分布地外的西南地區和華南地區也存在較大面積、較高適宜值的適生區。灰喜鵲之所以能在其自然分布地外的南方地區定殖,很大的原因在于灰喜鵲偏好暖濕的氣候,而南方地區相比其自然分布地更加溫暖潮濕。其次灰喜鵲較強的適應能力[17],也使其能在氣候條件完全不同的地方存活,引入地種群的氣候生態位可能已產生了重要變化。雖然引入灰喜鵲對控制松毛蟲等森林害蟲具有積極的意義,但其較強的適應力[17]和競爭力[18]也可能會對部分食蟲鳥類,尤其是同樣以松毛蟲為食的松鴉(Garrulusglandarius)、喜鵲(Picapica)、大杜鵑(Cuculuscanorus)及四聲杜鵑(Cuculusmicropterus)等鳥類造成影響[11, 51]。

猜你喜歡
物種模型
物種大偵探
物種大偵探
一半模型
吃光入侵物種真的是解決之道嗎?
英語世界(2023年10期)2023-11-17 09:18:18
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
回首2018,這些新物種值得關注
電咖再造新物種
汽車觀察(2018年10期)2018-11-06 07:05:26
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 激情五月婷婷综合网| 国产欧美综合在线观看第七页| 国产综合精品一区二区| 欧美日韩国产系列在线观看| 67194在线午夜亚洲 | 日本久久网站| 精品国产成人av免费| аⅴ资源中文在线天堂| 99热亚洲精品6码| av天堂最新版在线| 欧美h在线观看| 日本在线欧美在线| 亚洲天堂在线免费| 亚洲高清在线播放| WWW丫丫国产成人精品| 在线亚洲天堂| 欧美第二区| 青青操国产视频| 超清无码熟妇人妻AV在线绿巨人 | 国产丝袜第一页| 四虎亚洲国产成人久久精品| 欧美精品啪啪一区二区三区| 毛片网站免费在线观看| 久久狠狠色噜噜狠狠狠狠97视色| 污视频日本| 无码内射中文字幕岛国片| 日本午夜视频在线观看| 国产高清在线精品一区二区三区| 97在线观看视频免费| 无码日韩精品91超碰| 国产欧美日韩视频怡春院| 日韩国产黄色网站| 青草精品视频| 大乳丰满人妻中文字幕日本| 亚洲精品高清视频| 欧美a√在线| 色国产视频| 性欧美精品xxxx| 亚洲无码高清免费视频亚洲| 老色鬼欧美精品| 精品国产免费观看| 欧美一级视频免费| 视频在线观看一区二区| 91色综合综合热五月激情| 国产精品福利在线观看无码卡| 亚洲国产系列| 国产黑丝一区| 日韩高清无码免费| 亚洲天堂日韩av电影| 国产精品一区二区国产主播| 国产喷水视频| 日韩精品亚洲一区中文字幕| 婷婷五月在线| 亚洲AV无码一区二区三区牲色| 午夜精品久久久久久久2023| 亚洲欧美另类日本| 色哟哟国产精品| 国产成人精品一区二区秒拍1o| 国产成人91精品| 91啦中文字幕| 中文字幕免费在线视频| 小说 亚洲 无码 精品| av在线5g无码天天| 国产在线八区| www.youjizz.com久久| 午夜激情福利视频| 婷婷在线网站| 真实国产精品vr专区| 97视频精品全国在线观看 | 五月激情婷婷综合| 久久久噜噜噜| 国产精品流白浆在线观看| 91美女视频在线| 视频在线观看一区二区| 欧美精品高清| 久久精品66| 都市激情亚洲综合久久| 欧美一区精品| 中美日韩在线网免费毛片视频| 国产成人一级| 欧美伦理一区| 久草视频中文|