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

橋邊河大型底棲動物生境適宜性

2020-09-17 03:09:04粟一帆李衛(wèi)明李金京孫徐陽
生態(tài)學(xué)報 2020年16期
關(guān)鍵詞:物種模型

粟一帆,李衛(wèi)明,李金京,孫徐陽,胡 威

三峽大學(xué)水利與環(huán)境學(xué)院, 宜昌 443002

河流生態(tài)系統(tǒng)作為陸地生態(tài)系統(tǒng)和水生態(tài)系統(tǒng)在物質(zhì)、能量和信息交換中的紐帶,是自然生態(tài)系統(tǒng)的重要組成部分之一[1]。受水利水電工程建設(shè)和人類活動的影響,自然河流生態(tài)系統(tǒng)健康受到了嚴(yán)重威脅,引發(fā)了保護河流水生態(tài)水環(huán)境的熱潮[2- 3]。水生生物是河流生態(tài)系統(tǒng)的重要組成部分,其生境質(zhì)量狀況對流域水資源管理、河流生態(tài)恢復(fù)有重要的意義[4]。如何選擇合適的水生生物生境適宜性評價方法,建立準(zhǔn)確的生物與環(huán)境因子間的響應(yīng)關(guān)系,受到學(xué)術(shù)界的廣泛關(guān)注。

建立水生生物與環(huán)境因子之間的適宜性曲線,是水生生物生境適宜性研究最常用的方法[5]。最早運用適宜性曲線描述河流生境質(zhì)量的是專家評判法,利用專家知識和實踐經(jīng)驗對環(huán)境因子與生物選擇在[0, 1]之間進行評分。最常見的專家評判法有PHABSIM模型和模糊邏輯評價法[6- 7],當(dāng)數(shù)據(jù)集偏小且對評價結(jié)果的主客觀性要求不高時,專家評判法適用性較好;如果考慮的環(huán)境因子過多,數(shù)據(jù)量過大,專家評判法不再適用。數(shù)學(xué)統(tǒng)計方法、數(shù)學(xué)模型可以較好的解決這一問題,生境適宜性研究中常用的數(shù)學(xué)統(tǒng)計方法、模型有多項式擬合、廣義線性模型 (Generalize Linear Model,GLM)和廣義加性模型 (Generalized Additive Model,GAM)等[8- 10]。Martínez-Rincón等利用GAM模型和改進回歸樹模型對刺鲅生境適宜性進行擬合[11],Mouton等利用模糊邏輯模型和單變量適宜性曲線方法進行褐鱒魚生境適宜性評價[12]。鄭文浩等通過加權(quán)平均法研究了太子河流域大型底棲無脊椎動物生境適宜性[13],李若男等利用模糊邏輯模型研究了漓江光倒刺鲃的生境適宜性[14],易雨君等利用GAM模型對長江中華鱘進行了生境適宜性評價[15]。

水生生物生境適宜性研究方法雖然很多,但優(yōu)缺點各異,使用單一的方法很難滿足精度和適應(yīng)性的問題。為探尋精度高,適用性強的水生生物生境適宜性研究模型,本研究選取長江一級支流橋邊河為研究對象,利用典范對應(yīng)分析(Canonical correlation analysis,CCA)與獨立性分析探尋影響橋邊河大型底棲無脊椎動物生境變化的關(guān)鍵因子,采用多項式擬合和GAM模型分別擬合關(guān)鍵因子與底棲動物優(yōu)勢種的聯(lián)系,并進行擬合結(jié)果精度比較,得到精度高、適用性強的大型底棲動物生境適宜性研究模型。以期為河流生境質(zhì)量評價和生態(tài)修復(fù)提供理論基礎(chǔ)。

1 研究區(qū)域與方法

1.1 研究區(qū)域

橋邊河發(fā)源于宜昌市點軍與長陽交界的土城紅巖灣,于點軍社區(qū)注入長江,系長江右岸的一級支流。流域面積295 km2,主河道長40 km,最高海拔568 m,最低海拔55 m,是宜昌市重要的引用水源地之一。課題組于2019年1月和4月對橋邊河進行采樣調(diào)查,主要包括3個區(qū)域(圖1):點軍—橋邊(S1—S3),橋邊—土城(S4—S6),土城—源頭(S7—S8)。

圖1 橋邊河采樣點分布Fig.1 Distribution of sampling points in Qiaobian River

1.2 數(shù)據(jù)收集

1.2.1底棲動物采集

依據(jù)HJ 710.8—2014《生物多樣性觀測技術(shù)導(dǎo)則淡水底棲大型無脊椎動物》的規(guī)定,考慮河流形態(tài)、寬度、底質(zhì)類型等因素,采用500 μm D形網(wǎng)在每個采樣點兩岸上下游100 m內(nèi)的可涉水區(qū)域各采集3個分樣,每個分樣點采樣面積0.3 m2,將6個分樣點的樣品匯聚為一個樣品,每個采樣點采樣總面積均為1.8 m2。大型底棲無脊椎動物樣本經(jīng)0.5 mm鋼篩網(wǎng)篩洗后裝入500 mL標(biāo)本瓶,樣本加入10%甲醛溶液固定[16],在底部印有5 cm×5 cm網(wǎng)格線的白色搪瓷盤中進行底棲動物人工挑揀取樣,室內(nèi)挑揀并鑒定所有個體至可行的最低分類單元,其中寡毛類鑒定至綱,搖蚊類鑒定至亞科。樣品在分類計數(shù)時,若標(biāo)本損壞則只統(tǒng)計頭部,計算密度(個/m2)。稱重時,先用濾紙將樣品表面的水分吸干,再用萬分之一電子天平稱重,計算生物量(g/m2)[17]。

1.2.2生境因子采集

1.3 統(tǒng)計分析方法

1.3.1大型底棲無脊椎動物分布特征

運用優(yōu)勢度指數(shù)選取橋邊河底棲動物優(yōu)勢種[18],優(yōu)勢度指數(shù)計算見式(1):

P=(ni/N)fi

(1)

式中,ni為第i種的總個體數(shù),N為所有物種的總個體數(shù),fi為第i種在各站位出現(xiàn)的頻率,以P>0.02作為優(yōu)勢種判定閾值。

1.3.2生境因子選擇

利用束縛型排序(Constrained ordination)確定大型底棲無脊椎動物生境的水環(huán)境驅(qū)動因子[19]。大型底棲無脊椎動物群落的除趨勢對應(yīng)分析(Detrended correspondence analysis, DCA)最大梯度長度為3.6,單峰模型的典范對應(yīng)分析(Canonical correspondence analysis,CCA) 適用本研究[20]。為減少實驗分析誤差,選擇在3個或以上樣品中出現(xiàn),且相對豐度>1%的底棲動物種群進行CCA分析[21],為保證CCA分析時物種數(shù)大于環(huán)境因子數(shù),本文選取相對豐度排名前10的物種進行分析。采用Pearson相關(guān)性分析法與Spearman相關(guān)性分析法剔除相關(guān)性高的生境因子。運用CANOCO 5進行DCA與CCA分析,運用SPSS 22進行相關(guān)性分析。

1.3.3底棲動物生境擬合

以優(yōu)勢種在點位出現(xiàn)的個數(shù)為自變量,以關(guān)鍵生境因子為因變量,利用多項式擬合和廣義加性模型構(gòu)建大型底棲無脊椎動物與生境因子的隸屬度模型。

多項式擬合的一般公式[22]為

f(c0+c1x+...+cnxn)n=0,1,2…

(2)

式中,cn為多項式擬合系數(shù)。

廣義加性模型(GAM)的一般公式[23]為

g(μ(Y))=β0+f1(x1)+...+fm(xm)

(3)

式中,g()函數(shù)為連接函數(shù),μ(Y)為Y的期望,β0為截距,fi()是非參數(shù)光滑函數(shù)。

利用Origin 9 進行多項式擬合,利用R軟件(version 3.3.2)的“mgcv”工具包進行GAM分析,基礎(chǔ)數(shù)據(jù)處理與統(tǒng)計在Excel 2016中完成。

2 結(jié)果

2.1 底棲動物組成與生境因子

2.1.1底棲動物組成

調(diào)查共采集到底棲動物1152頭,隸屬于3門5綱19科。其中節(jié)肢動物種類最多,占總物種數(shù)的68.4%,其次為軟體動物,占總物種數(shù)的26.3%,環(huán)節(jié)動物最少,占總物種數(shù)的5.3%。優(yōu)勢度指數(shù)計算結(jié)果顯示,橋邊河主要優(yōu)勢物種有黃色羽搖蚊、溪蟹、方格短溝蜷、圓頂珠蚌、耳蘿卜螺、梨形環(huán)棱螺和河蜆。各優(yōu)勢種沿程分布見表1,黃色羽搖蚊、方格短溝蜷、圓頂珠蚌和耳蘿卜螺主要分布于點軍—橋邊,譚氏泥蟹主要分布于土城—源頭段。由于河蜆在河流各個區(qū)域均有出現(xiàn),本研究以河蜆為指示物種進行研究。

表1 橋邊河優(yōu)勢種分布

2.1.2生境因子

橋邊河水深等生境因子的統(tǒng)計結(jié)果見表2。總體來看,橋邊河流域生境狀況較好,TN、TP平均值分別為0.470、0.078 mg/L,未超過GB 3838—2002 Ⅱ類標(biāo)準(zhǔn)限值。但局部地區(qū)的生境狀況較差,TN的最大值 (1.064 mg/L) 超過GB 3838—2002Ⅱ類標(biāo)準(zhǔn)限值2倍以上,NO3-N與TP的最大值也超過了Ⅱ類標(biāo)準(zhǔn)限值。橋邊河沿程水深不大,最淺的地方為8號樣點,水深僅0.2 m。橋邊河河床底質(zhì)粒徑空間異質(zhì)性明顯,下游底質(zhì)多為細沙為主;靠近上游區(qū)域河床底質(zhì)以卵礫石為主,D50在50 mm左右。

2.2 化學(xué)因子與底棲動物的關(guān)聯(lián)性

底棲動物物種與化學(xué)因子的CCA分析見表3與圖2,前四軸的變異解釋率累積為47.72%,蒙特卡羅檢驗(n=499)對所有軸均有顯著性意義(P<0.05)。與第二排序軸相關(guān)性較高的化學(xué)因子有TN、NO3-N、NH4-N、CODMn、pH、DO。TN、NO3-N、NH4-N、CODMn與第二排序軸呈顯著負相關(guān),pH,DO與第二排序軸呈顯著正相關(guān)。優(yōu)勢種河蜆與DO、pH的夾角為銳角,河蜆與DO呈正相關(guān);河蜆與TN、NO3-N、NH4-N、CODMn的夾角為鈍角,河蜆與這些因子呈負相關(guān)。

表2 橋邊河生境因子值

表3 化學(xué)因子CCA相關(guān)性

2.3 物理因子與底棲動物的關(guān)聯(lián)性

底棲動物物種與物理因子的典范對應(yīng)分析見表4、圖3,前四軸的變異解釋率累積為67.4%,蒙特卡羅檢驗(n=499)對所有軸均有顯著性意義(P<0.05)。與第二排序軸相關(guān)性較高的因子有D50、Dep、Tur。Tem與第一軸的相關(guān)性較高,但第一軸解釋類群變化的貢獻相較第二軸較小,所以Tem的貢獻較小。優(yōu)勢種河蜆與D50、Dep、Tur的夾角為鈍角,表明河蜆與這些環(huán)境因子呈負相關(guān)關(guān)系。

表4 物理因子CCA相關(guān)性

2.4 生境因子與底棲動物的響應(yīng)關(guān)系

對生境因子進行獨立性分析,首先進行K-S正態(tài)檢驗,對符合正態(tài)分布的指標(biāo)進行Pearson相關(guān)分析,對不符合正態(tài)分布的指標(biāo)進行Spearman相關(guān)分析,以相關(guān)系數(shù)|R|>0.75為闕值[24],相關(guān)性分析結(jié)果見表5。水質(zhì)因子中,CODMn與NH4-N、TN與NO3-N、TN與NH4-N相關(guān)系數(shù)均超過0.75,具有較高的相關(guān)性。鑒于數(shù)據(jù)的易得性,剔除NH4-N與NO3-N;pH與DO存在較高的相關(guān)性,橋邊河沿程pH變異度不高,因此剔除pH因子。環(huán)境因子中,D50與Tur和Dep之間存在較高的相關(guān)性,因此剔除D50。

圖2 化學(xué)因子CCA排序圖Fig.2 Canonical correspond analysis of chemical factors SP1為Chironomus flaviplumus;SP2為Ilyoplax deschampsi;SP3為Semisulcospira cancellata;SP4為Unio douglasiae;SP5為Radix auricularia;SP6為Bellamya purificata;SP7為Corbicula fluminea SP8為Corophium acherusicum costa;SP9為Chinese white prawn;S10為Bellamya aeruginosa; TN:總氮,Total nitrogen;NO3-N:硝態(tài)氮,Nitrate;NH4-N:銨態(tài)氮,Ammonium;TP:總磷,Total phosphorus;PO4-P:正磷酸鹽,Orthophosphate;CODMn:化學(xué)需氧量,Chemical oxygen demand;DO:溶解氧,Dissolved oxygen;pH:氫離子濃度指數(shù),Hydrogen ion concentration

圖3 環(huán)物理子CCA排序圖Fig.3 Canonical correspond analysis of physical factors SP1為Chironomus flaviplumus;SP2為Ilyoplax deschampsi;SP3為Semisulcospira cancellata;SP4為Unio douglasiae;SP5為Radix auricularia;SP6為Bellamya purificata;SP7為Corbicula fluminea.;SP8為Corophium acherusicum costa;SP9為Chinese white prawn;S10為Bellamya aeruginosa..;Tur:濁度,Turbidity;Tem:溫度,Temperature;Vel:流速,Flow velocity;Dep:水深,Water depth;D50:河床中值粒徑,Median bed size

利用多項式擬合和廣義加性模型構(gòu)建底棲動物與生境因子的隸屬度關(guān)系,擬合結(jié)果見圖4。由于Tur數(shù)據(jù)變異度過大,采用lg(x+1)進行數(shù)據(jù)處理。河蜆的最適宜CODMn含量為1.228 mg/L,最適宜的TN含量為0.269 mg/L,最適宜的DO含量為11.170 mg/L,最適宜的Tur為1.130 NTU,最適宜的Dep為0.3 m。

圖4 橋邊河生境因子隸屬度Fig.4 Membership degree of habitat factors in Qiaobian River

3 討論

3.1 生境適宜性模擬對比

GAM廣義加性模型與傳統(tǒng)多項式擬合模型的擬合結(jié)果對比見表6。CODMn、TN、Tur、Dep的擬合結(jié)果顯示GAM廣義加性模型優(yōu)于傳統(tǒng)多項式擬合模型,非線性擬合情況下GAM廣義加性模型擬合結(jié)果較好;DO模擬結(jié)果顯示兩種方法的R2大致接近,多項式擬合結(jié)果較優(yōu)于GAM擬合結(jié)果,線性擬合情況下多項式擬合結(jié)果較好。TN擬合結(jié)果兩者接近,GAM擬合結(jié)果較優(yōu),多項式擬合結(jié)果出現(xiàn)過擬合情況,TN處于0.65—1.00 mg/L時,底棲動物適宜性出現(xiàn)負值,多項式擬合過擬合現(xiàn)象在Tur擬合中也有體現(xiàn)。多項式擬合是預(yù)測水生生物物種豐富度和棲息地適宜性的常用方法[25-26],但由于其處理非線性問題時常出現(xiàn)多重共線問題和過擬合等現(xiàn)象,因此,多項式擬合方法有很大的局限性。廣義加性模型因其靈活性和處理非線性問題的憂越性,被廣泛應(yīng)用于預(yù)測河流中魚類、大型底棲無脊椎動物、沉水植物的分布中[27-28]。Myers等[29]的研究表明,當(dāng)線性模型不能很好表征物種與各因子相關(guān)性時,GAM等模型可以進行相應(yīng)替代,但需注意GAM模型難以輸出具體擬合函數(shù)的問題。本研究顯示,線性擬合情況下兩種方法均適用;非線性擬合情況下GAM廣義加性模型較優(yōu)于傳統(tǒng)多項式擬合模型,且GAM模型在處理離散程度大的數(shù)據(jù)集時,可以很好的規(guī)避多項式擬合過程中出現(xiàn)的過擬合現(xiàn)象。

3.2 生境因子與大型底棲無脊椎動物聯(lián)系

CCA排序與Pearson相關(guān)性分析結(jié)果顯示,影響橋邊河大型底棲無脊椎動物優(yōu)勢種生境質(zhì)量的主要因子為CODMn、TN、DO、Tur、Dep。CODMn、TN、Tur和Dep與河蜆的分布成負相關(guān)關(guān)系,DO與河蜆的分布成顯著正相關(guān)關(guān)系,分析結(jié)果與段學(xué)花等人的研究結(jié)果類似[30-32]。

CODMn是測定河流有機物含量的重要指標(biāo),主要應(yīng)用于河流污染評估和工業(yè)廢水性質(zhì)的研究以及廢水處理廠的運行管理[33]。廣義加性模型與多項式擬合模型結(jié)果顯示,橋邊河流域河蜆的最適宜CODMn含量為1.228 mg/L,多項式擬合結(jié)果顯示河蜆生境適宜性隨著CODMn升高而降低,當(dāng)CODMn超過2.159 mg/L時,河蜆生境適宜性降低為0。河蜆多分布河流近岸帶,靠近工業(yè)廢水排污口、污水處理廠排污口,水體CODMn較高易造成水體缺氧和富營養(yǎng)化[34],因此隨著橋邊河水域CODMn增高,河蜆生境適宜性顯著降低。廣義加性模型顯示,當(dāng)CODMn升至1.6—1.8 mg/L時,河蜆生境適宜性下降趨勢有一定減緩,分析發(fā)現(xiàn)這可能與河蜆是一種中度耐污種有一定關(guān)系。大型底棲無脊椎動物耐污值研究結(jié)果顯示,蜆類耐污值約為6—8(總分10分),可以在一定的污染環(huán)境中生存[35],且一定濃度的污水可以為河蜆提供生產(chǎn)所需有機物,因此,適宜性模擬結(jié)果出現(xiàn)較小波動。然而較高的CODMn依舊對河蜆的生產(chǎn)生活呈顯著的抑制作用。

表5 相關(guān)分析結(jié)果

表6 隸屬度對比分析

TN是判別河流富營養(yǎng)化的主要因子之一[36]。廣義加性模型與多項式擬合模型結(jié)果顯示,橋邊河流域河蜆的最適宜TN含量為0.269 mg/L,多項式擬合結(jié)果顯示隨著TP升高河蜆生境適宜性急劇降低。水體中N等營養(yǎng)物質(zhì)含量過多將導(dǎo)致水體富營養(yǎng)化,河流生態(tài)系統(tǒng)失衡,藻類等單一物種迅速繁殖,嚴(yán)重擠占河流中大型底棲無脊椎動物生態(tài)位,抑制底棲動物的生長繁殖[37]。水體中大量營養(yǎng)物質(zhì)經(jīng)水域微生物分解,會消耗水中溶氧,也會抑制底棲動物生長[38]。多項式擬合結(jié)果當(dāng)TN升至0.8—1.0 mg/L時,河蜆生境適宜性出現(xiàn)一定程度回升,這可能與底棲動物食物鏈有一定聯(lián)系。當(dāng)水體有機物增多時,水體中的藻類能夠利用環(huán)境中的氮素等顯著提高其初級生產(chǎn)力,為河蜆提供充足食物,河蜆生境適宜性有一定提升。

DO是水生生物賴以生存的必要生命因子之一,水體中DO的多寡極大影響水生生物的攝食、繁殖等生命過程。廣義加性模型與傳統(tǒng)多項式擬合模型結(jié)果顯示,橋邊河流域河蜆的最適宜DO含量為11.170 mg/L,橋邊河流域DO處于有益于河蜆生長的范圍,且一定范圍內(nèi)DO的增高,河蜆生境適宜性有進一步提升的趨勢。河蜆等原生底棲動物物種可以直接利用水中的溶解氧進行生長繁殖,DO是河蜆生境分布的直接驅(qū)動因素之一,因此河蜆生境適宜性與DO呈顯著正相關(guān);另一方面,DO的增加可以抑制水體的反硝化過程[39],亞硝酸鹽、硝酸鹽等不易在反硝化細菌作用下生成氨氮,減弱了河蜆等的氨氮中毒風(fēng)險[40]。同時溶解氧的增加將提高河蜆的食物同化率,促進河蜆的增長。

水深是研究大型底棲無脊椎動物生境適宜性的常用因子之一,一般認為隨著水深的增加,底棲動物群落密度等會相應(yīng)遞減[41]。橋邊河流域大型底棲無脊椎動物生境適宜性研究顯示,水深因子與河蜆生境質(zhì)量成負相關(guān),這與Buss等的研究結(jié)果類似[42]。廣義加性模型與傳統(tǒng)多項式擬合模型結(jié)果顯示,橋邊河流域河蜆的最適宜水深為0.3 mg/L。水深因子主要通過影響其他生境因子的變化,來作用于大型底棲無脊椎動物生境質(zhì)量的變化。隨著水深的增加,深層水流的光照驟減[43],初級生產(chǎn)力降低,河蜆食物來源減少,生物量受影響。同時深水區(qū)DO濃度較小,厭氧細菌的強烈活動增加反硝化過程,抑制河蜆的生長繁殖[44]。

河流中泥沙、黏土、有機物等懸浮物質(zhì)含量過多將造成河流Tur增高,Tur因子是影響大型底棲無脊椎動物生境質(zhì)量的重要因子之一[45]。廣義加性模型與傳統(tǒng)多項式擬合模型結(jié)果顯示,橋邊河流域河蜆的最適宜Tur為1.130 NTU,在一定范圍類Tur的升高導(dǎo)致河蜆棲息地適宜性的降低,但當(dāng)Tur在3—5 NTU范圍內(nèi)時,河蜆生境適宜性出現(xiàn)短暫的回升現(xiàn)象。濁度較高的水體河床底質(zhì)大多為淤泥,河床基質(zhì)穩(wěn)定性較差,大型底棲無脊椎動物生境質(zhì)量較低[46]。同時濁度較高的河流水體透明度較低,水生藻類等水生植物光合作用受到抑制,水體DO含量降低,抑制底棲動物生長。而河床泥沙再懸浮過程會釋放沉積物中蘊藏的大量有機物[47],小范圍的河床擾動給河蜆的生長提供了充足的營養(yǎng)物質(zhì),短期內(nèi)可能造成河蜆生境適宜性出現(xiàn)上升現(xiàn)象。

4 結(jié)論

(1) 影響橋邊河大型底棲無脊椎動物優(yōu)勢種生境質(zhì)量的主要因子為CODMn、TN、DO、Tur、Dep;其中CODMn、TN、Tur、Dep與其分布成負相關(guān),DO與其分布成正相關(guān)。

(2) 橋邊河流域河蜆最適宜CODMn濃度為1.228 mg/L,CODMn升高導(dǎo)致河蜆生境適宜性降低;最適宜TN濃度為0.269 mg/L,TN升高導(dǎo)致河蜆生境適宜性急劇下降;最適宜DO濃度為11.170 mg/L,DO與大型底棲無脊椎動物境適宜性呈顯著正相關(guān)。最適宜Dep為0.3 m,Dep因子與河蜆生境質(zhì)量成負相關(guān)。最適宜Tur為1.130 NTU,一定范圍內(nèi)Tur的升高導(dǎo)致河蜆生境適宜性的降低,但超過某一限定值時,河蜆生境適宜性存在回升的趨勢。

(3) GAM廣義加性模型與傳統(tǒng)多項式擬合模型的擬合結(jié)果顯示,線性擬合情況下兩種方法均適用;非線性擬合情況下GAM廣義加性模型較優(yōu)于傳統(tǒng)多項式擬合模型,且GAM模型在處理離散程度大的數(shù)據(jù)集時,可以很好的規(guī)避多項式擬合過程中出現(xiàn)的過擬合現(xiàn)象。

猜你喜歡
物種模型
物種大偵探
物種大偵探
一半模型
吃光入侵物種真的是解決之道嗎?
英語世界(2023年10期)2023-11-17 09:18:18
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
回首2018,這些新物種值得關(guān)注
電咖再造新物種
汽車觀察(2018年10期)2018-11-06 07:05:26
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 国产免费观看av大片的网站| 综合亚洲色图| 国产精品福利尤物youwu | 日本爱爱精品一区二区| 成人在线亚洲| 91久久性奴调教国产免费| 国内毛片视频| 亚洲品质国产精品无码| 2021国产精品自产拍在线| 免费a级毛片视频| 亚洲精品国产精品乱码不卞 | 亚洲一级毛片在线观| 国产一级在线播放| 国产av剧情无码精品色午夜| 无码精品福利一区二区三区| 国产区福利小视频在线观看尤物| 国产黄在线免费观看| 久久精品丝袜高跟鞋| 中文字幕乱码二三区免费| 不卡网亚洲无码| 色综合天天综合| 久久动漫精品| 国产成+人+综合+亚洲欧美 | 亚洲综合精品香蕉久久网| 国禁国产you女视频网站| 一本一道波多野结衣一区二区| 日本高清视频在线www色| 亚洲一区二区视频在线观看| 久久美女精品| 亚洲an第二区国产精品| 国产jizz| 蜜桃臀无码内射一区二区三区| 在线观看免费国产| 欧美亚洲欧美区| 丁香婷婷久久| 国产精品无码一区二区桃花视频| 丁香婷婷久久| 在线欧美a| 麻豆精品视频在线原创| 中文字幕永久在线看| 亚洲精品天堂自在久久77| av一区二区三区高清久久| 欧美激情视频二区三区| 99精品福利视频| 一本大道香蕉久中文在线播放| 午夜老司机永久免费看片 | 日韩高清无码免费| 成人蜜桃网| 国产视频欧美| 婷婷成人综合| 亚洲精品在线影院| 18禁不卡免费网站| www.亚洲一区| 91小视频在线播放| 蜜桃臀无码内射一区二区三区| 午夜不卡视频| 亚洲一级毛片在线观| 国产精品网址在线观看你懂的| 免费毛片全部不收费的| 亚洲国产中文欧美在线人成大黄瓜| 午夜无码一区二区三区| 中文字幕人妻无码系列第三区| 国产精品亚洲欧美日韩久久| 国产乱视频网站| 国产玖玖视频| 二级特黄绝大片免费视频大片| 亚洲综合婷婷激情| 在线精品视频成人网| 亚洲伊人天堂| 一本综合久久| 美女视频黄频a免费高清不卡| 中文字幕久久亚洲一区| 亚洲综合色在线| 日本欧美一二三区色视频| 亚洲永久色| 尤物精品国产福利网站| 四虎永久免费网站| 国产女人18水真多毛片18精品 | 国产福利免费视频| 国产午夜无码片在线观看网站| 精品国产毛片| 永久免费精品视频|