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

凈多樣化速率和進化時間對虎耳草目科間物種多樣性差異的影響

2022-11-16 13:28:50顧嘉豪張粒毫張皓昱仵天晴黃林青程瑞靜王慶剛徐曉婷
廣西植物 2022年10期
關鍵詞:物種差異研究

顧嘉豪, 張粒毫, 張皓昱, 仵天晴, 黃林青, 程瑞靜,2, 徐 瑩,2, 王慶剛, 徐曉婷,2*

( 1. 四川大學 生命科學學院, 成都 610065; 2. 四川大學 生物資源與生態環境教育部重點實驗室, 成都 610065; 3. 四川大學 生物醫學工程學院, 成都 610065; 4. 中國農業大學 資源與環境學院, 北京 100083 )

生物類群之間的物種總數的差異是自然界中普遍存在的現象,但是這一現象的形成機制仍然是未解之謎(Ricklefs & Renner,1994; Barraclough & Savolainen, 2010; Rabosky, 2010; Wiens, 2011; Robosky & Adams, 2012; Mereau & Bell, 2013; Machac, 2020; Li & Wiens, 2019)。從生物類群的宏觀進化角度來看,進化時間假說(time-for speciation hypothesis)和多樣化速率假說(diversification rate hypothesis)是解釋類群間物種多樣性差異的兩個最為核心的假說(Scholl & Wiens, 2016; Li & Wiens, 2019)。進化時間假說認為,進化時間越長,物種能積累的物種多樣性越高,因此進化的時間差異是導致類群間物種多樣性差異的主要原因。換句話說,一個起源古老的類群,比新近起源的類群擁有更長的進化時間,因此可以積累更多的物種。多樣化速率假說認為,在生物演化歷史中,不同類群種化速率和滅絕速率的動態變化導致的凈多樣化速率差異是引起物種多樣性變化的主要原因(Scholl & Wiens, 2016)。也就是說,物種形成速率高而滅絕速率低的類群具有更高的凈多樣化速率,因此具有更多的物種。

目前,進化時間假說和多樣化速率假說在物種多樣性形成中的相對作用還存在很多爭議(Marin et al.,2018)。進化時間假說成立的前提是演化過程中凈多樣化速率在不同類群之間無差異或差異極小。然而越來越多的研究發現,凈多樣化速率在不同類群之間有差異。特別是經歷過大型滅絕事件的類群,物種數目會驟然下降,從而導致物種多樣性的突然喪失,物種數目減少。例如,Rabosky 等(2012)對多細胞真核生物1 397個主要分支的研究發現,分支的分化時間,即干齡(stem age),與分支所包含的物種數之間沒有顯著關系,因而否定了進化時間假說。Pyron和Wiens (2013)發現,許多古老兩棲動物分支具有較低的物種多樣性,但是一些年輕的分支具有較多的物種。但是該研究并沒有直接驗證物種多樣性與多樣化速率之間的關系。而McPeek和Brown (2007)對不同動物類群物種多樣性的研究卻發現分支年齡是物種多樣性的決定因素。Marin和Hedges (2016)對不同類群的冠齡(crown age)與物種豐富度關系的研究也表明兩棲動物、鳥類和哺乳動物的物種豐富度與其最近共同祖先的年齡(age)顯著正相關。為進一步驗證進化時間和多樣化速率假說,Scholl和Wiens (2016)使用生命之樹(tree of life)對不同生物類群不同分類學等級的物種多樣性的差異進行了研究。研究中發現,在不同的分類等級,例如,門、目和科等,物種多樣性主要是由多樣化速率決定,而與進化時間無顯著相關關系,甚至出現負相關關系。古老的類群物種多樣化速率低,而年輕類群的多樣化速率高可能是導致該現象的主要原因。并且,進化時間與物種多樣化速率的關系在不同分類等級和不同的生物類群中的差異可能導致不顯著的物種多樣化速率與物種多樣性的關系。綜上,進化時間假說和多樣化速率假說在物種多樣性形成中的作用仍然存在爭議。

虎耳草目(Saxifragales)包含了15個科的共約3 000個種,包括了喬木、灌木、多年生或一年生草本植物、多肉植物和水生植物等多種生活型,廣泛分布于全球不同的生態系統中(Soltis et al., 2013)。虎耳草目內部科的分化時間從早白堊紀晚期至第三紀早期,且科內物種數目存在巨大差異,最大的景天科(Crassulaceae)包含了大約1 600個種,而鎖陽科(Cynomoriaceae)和四心木科(Tetracarpaeaceae)、隱瓣藤科(Aphanopetalaceae)卻只有不到10個物種。同時,虎耳草目的分子系統學研究基礎扎實,基于301個核基因建立的系統發育骨架,確定了各科之間的系統發育關系和分化時間(Folk et al., 2019)。Folk 等 (2019)還結合小片段構建了種級水平的系統發育樹,物種覆蓋率達當時公認物種的70%以上,為準確估計物種的多樣化速率奠定了基礎。因此,虎耳草目是研究科間物種豐富度差異、驗證進化時間和多樣化速率假說的理想類群。本研究采用系統發育廣義最小二乘模型(phylogenetic generalized least squares,PGLS)分析了虎耳草目科的物種多樣性與科的冠齡、干齡和多樣化速率的關系,發現了物種多樣性與多樣化速率有顯著的正相關關系,與進化時間沒有顯著關系,且進化時間與多樣化速率的綜合效應對這種物種豐富度差異有更強的解釋性。本研究結果發現了多樣化速率在虎耳草目物種多樣性形成中的重要作用,支持了物種多樣性格局的多樣化速率假說,同時也支持了進化時間和多樣化速率的綜合效應。

1 材料與方法

1.1 虎耳草目系統發育樹

本研究采用的系統發育樹主要基于Folk等(2019)發表的虎耳草目的系統發育樹和徐瑩等(2021)構建的虎耳草屬的系統發育樹。Folk等(2019)首先基于來自627個物種的301 個核基因序列,構建了虎耳草目主要分支的系統發育樹并通過化石和分子鐘定年確定了主要分支的分化時間。之后,基于GenBank中下載虎耳草目的核基因ITS片段和葉綠體基因組的matK等24個基因片段,建立了包含15個科1 455個物種的系統發育樹。徐瑩等(2021)構建的系統發育樹包含了353個虎耳草屬的物種,涵蓋了Folk等(2019)構建的虎耳草目系統發育樹中的所有虎耳草屬物種,虎耳草屬的分化時間也與Folk等(2019)一致。因此,本研究將兩棵系統發育樹進行整合,在R語言中將徐瑩等(2021)構建的虎耳草屬系統發育樹替換了Folk等(2019)虎耳草目系統發育樹中虎耳草屬分支(Gordon, 1986;Bininda-Emonds, 2004),最終整合得到了一棵包含1 539個種的虎耳草目系統發育樹,涵蓋虎耳草目下全部15個科,并且各科均具有較高覆蓋率(表1) 。

表 1 虎耳草目15個科有關數據統計Table 1 Relevant data statistics of 15 families of order Saxifragales

1.2 科的物種多樣性

科的物種多樣性數目按照生物物種名錄數據庫(COL,https://www.catalogueoflife.org/,訪問時間:2021年5月)中每個科接受的物種名稱進行統計,變種和亞種等種下等級不統計,雜種不統計。由于物種多樣性數目不符合正態分布,參考Rabosky 等(2012)的分析方法對科的物種多樣性進行了以10為底的對數轉換以改善擬合效果。

1.3 科的干齡和冠齡計算

冠齡與干齡是針對系統發育樹中分支年齡的兩個不同的概念(圖1)。前者代表該類群現存物種最近共同祖先的年齡,而后者代表該類群最近共同祖先與其姊妹類群的最近共同祖先的年齡。有學者研究認為在取樣較好的情況下,冠齡用來做數據統計分析更為合理,因為冠齡可以消除由于進化停滯或者滅絕事件導致的長枝的影響(Stadler et al., 2014; Sanchez-Reyes et al., 2017)。為了更加全面充分地評估進化時間與物種豐富度的關系,本研究同時分析了各個科干齡以及冠齡與其物種多樣性的關系。

圖 1 干齡和冠齡的示意圖(單位:百萬年)Fig. 1 Schematic diagram of stem age and crown age (Unit: Myr)

為了提取不同科的冠齡和干齡,本研究利用R語言ape程輯包(Paradis & Schliep, 2019;R Core Team, 2021)中getMRCA()函數從系統發育樹上找到每個科內的物種對應的共同祖先的節點,該節點的年齡即為該科的冠齡,該節點的父節點即為該科的干齡。虎耳草目各個節點的年齡使用R 程序中的ape程輯包中的branching.times ( )函數提取。各個科的冠齡節點的編號使用ape程輯包中的getMRCA()函數提取。根據每個科的冠齡節點和干齡節點編號(nodelable),分別獲取到15個科的冠齡和干齡。

1.4 科的多樣化速率

本研究利用宏觀進化貝葉斯分析法(Bayesian analysis of macroevolutionary mixtures, BAMM),使用R語言中的程輯包BAMMtools獲取科的多樣化速率(Rabosky et al., 2014)。馬爾可夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)運行長度為1 000萬代,每1 000代取樣一次。運行結束后,在R語言中使用程輯包coda檢查鏈的收斂性和有效樣本(effective sampling sizes, ESS)的大小。ESS檢查結果大于200,說明MCMC過程有足夠的獨立樣本,估計的結果比較穩定。

虎耳草目15個科的物種多樣化速率使用BAMMtools程輯包的getCladeRates( )函數提取。該方法提取了科內所有分支的枝長加權平均的種化速率(speciation rate)和滅絕速率(extinction rate)。凈多樣化速率等于種化速率與滅絕速率的差值。由于滅絕速率的估計存在較大的不確定性,本研究中同時分析了種化速率和凈多樣化速率與物種多樣性的關系。

此外,根據各個科是否為溫帶適應類群,利用ape程輯包中的drop.tip()函數將虎耳草目系統發育樹拆分為兩部分,即溫帶適應類群和常綠喬木類群[常綠喬木類群為圍盤樹科(Peridiscaceae)、虎皮楠科(Daphniphyllaceae)、鼠刺科(Iteaceae)和蕈樹科(Altingiaceae);虎耳草目下其余11個科為溫帶適應類群],利用BAMMtools程輯包中的plot.bammdata()函數繪制出凈多樣化速率隨時間的變化圖像。

1.5 物種多樣性與多樣化速率和進化時間的關系

本研究使用一元線性回歸模型(linear regression model,LM)和系統發育廣義最小二乘模型(PGLS)對物種豐富度與多樣化速率和進化時間的關系進行了分析。為驗證多樣化速率假說,本文同時建立了種化速率和凈多樣化速率與物種多樣性的一元線性回歸模型和PGLS模型。同時,也分別建立了冠齡和干齡與物種多樣性的一元線性回歸關系和PGLS模型來驗證進化時間假說。由于多樣化速率和進化時間可能同時對物種多樣性造成影響,因此,本研究也考慮了多樣化速率和進化時間的交互效應對物種多樣性的影響。交互效應中,進化時間由冠齡來代表。LM使用R語言中的基礎函數lm ( )完成。PGLS使用R程輯包中的 nlme的gls()函數完成。PGLS分析中使用了一個虎耳草目的科級系統發育樹。該系統發育樹是使用R語言的ape程輯包中的drop.tip()函數將物種數目超過1個的科中的物種進行去除,保證每個科只有1個代表物種。

2 結果與分析

2.1 虎耳草目科的干齡和冠齡與物種豐富度的關系

虎耳草目15個科的干齡和冠齡變化范圍差異很大。最年輕的連香樹科(Cercidiphyllaceae)冠齡僅為3百萬年(Myr),而最古老的圍盤樹科的冠齡則超過了8千萬年(表1)。這充分體現了虎耳草目內部不同科在進化時間尺度上的差異性(表1)。

LM分析的結果表明冠齡對物種豐富度的解釋率可達到0.261,但是干齡對物種豐富度的解釋率只有0.164(圖2)。PGLS分析的結果與LM一致,冠齡與物種豐富度的斜率為0.015,略低于干齡與物種豐富度的斜率(β=0.019)。但無論是冠齡還是干齡對物種豐富度的影響都不顯著(P>0.05)。該結果說明虎耳草目內15個科的物種豐富度與科的冠齡和干齡沒有顯著相關性,進化時間并不是影響科之間的物種多樣性差異的主要因素。

圖中右下角為LM模型的解釋率(R2)和顯著性檢驗的P值。The lower right corner are the LM model’s explanatory rate (R2) and P value of significance test.圖 2 虎兒草目15個科的冠齡與物種豐富度(A)及干齡與物種豐富度(B)的系統發育廣義最小二乘模型(PGLS)與線性回歸模型(LM)擬合Fig. 2 Phylogenetic generalized least squares (PGLS) and linear regression model (LM) for crown age and species richness (A) and stem age and species richness (B) of 15 families of order Saxifragales

2.2 虎耳草目科的物種多樣化速率與物種豐富度的關系

虎耳草目凈多樣化速率較高的科以溫帶和高山類群為主,集中在虎耳草目中以草本植物為主的草本分支中,例如茶藨子科(Grossulariaceae)、小二仙草科(Haloragaceae)和景天科等科的種化速率和凈多樣化速率都比較高(表1,圖3),并且以1 500萬年以來增長較為明顯(圖3,圖4)。物種主要為常綠喬木的科的種化速率和凈多樣化速率相對較低 (表1,圖3, 圖4)。

綠色標注的科是常綠喬木類群。The families shown in green are evergreen trees.圖 3 虎耳草目凈多樣化速率-時間圖Fig. 3 Net diversification rate-time map of order Saxifragales

圖 4 虎耳草目中常綠喬木和溫帶適應類群隨時間的凈多樣化速率變化曲線Fig. 4 Net diversification rate curves over time for evergreen trees and temperate adaptive taxonomic groups of order Saxifragales

LM和PGLS分析發現,虎耳草目科的種化速率和凈多樣化速率都與物種豐富度正相關(圖5)。但種化速率與物種豐富度的關系在LM中顯著(R2=0.280,β=7.319,P<0.05),而在PGLS模型中接近顯著(β=6.425,P<0.05),而這兩種方法發現凈多樣化速率與物種豐富度均有顯著正相關性(PGLS:β=9.790,P=0.031;LM:R2=0.380,β=11.791,P<0.05)。

陰影代表LM模型的95%的置信區間。圖中右下角為LM模型的解釋率(R2)和顯著性檢驗的P值。下同。The shadow represents 95% confidence interval of the LM model. In the lower right corner are the LM model’s explanatory rate (R2) and P value of significance test. The same below.圖 5 虎兒草目15個科的種化速率(A)和凈多樣化速率(B)與物種豐富度的系統發育廣義最小二乘模型(PGLS)與線性回歸模型(LM)擬合Fig. 5 Phylogenetic generalized least squares model (PGLS) and linear regression model (LM) fitted speciation rates (A) and net diversification rate (B) with the species richness of 15 families of order Saxifragales

2.3 進化時間與凈多樣化速率的綜合效應與科物種豐富度的關系

進化時間與凈多樣化速率對物種豐富度的共同作用在LM和PGLS模型中都顯著(LM:R2=0.603,β=0.297,P<0.05;PGLS:β=0.287,P<0.05) (圖6)。LM模型中,進化時間和凈多樣化速率的共同作用解釋了物種多樣性變化的60.3%,PGLS模型的斜率也比單獨使用進化時間和凈多樣化速率高。該結果支持了進化時間與凈多樣化速率的綜合效應對于虎耳草目科間物種豐富度差異的影響,且比凈多樣化速率這一單一變量對其的影響解釋率更高。

圖 6 虎兒草目15個科進化時間和凈多樣化速率的綜合效應與物種豐富度的系統發育廣義最小二乘模型(PGLS)與線性回歸模型(LM)擬合Fig. 6 Phylogenetic generalized least squares model (PGLS) and linear regression model (LM) fitted the interaction effects of evolutionary time and net diversification rate with the species richness of 15 families of order Saxifragales

3 討論與結論

在本研究中,PGLS 分析發現虎耳草目下15個科的物種多樣性與凈多樣化速率呈顯著的正相關關系,但是與科的冠齡正相關關系不顯著。LM中,凈多樣化速率對物種豐富度的解釋率達0.380。由于 Scholl和Wiens (2016)發現凈多樣化速率與年齡的交互作用可能是導致物種多樣性差異的原因,因此本研究也檢驗了物種凈多樣化速率與冠齡的交互作用,并發現交互作用對物種豐富度的解釋率達到0.603。因此,本研究認為凈多樣化速率的差異是主導物種多樣性差異的主要因素,進化時間與多樣化速率的共同作用在虎耳草目物種多樣性差異中起到了次要作用。

本研究的結果不支持進化時間假說,這與Rabosky等(2012)和Scholl和Wiens (2016)得到的結論一致,即不同類群間物種豐富度的差異并非是由進化時間的差異引起。盡管Pyron和Wiens (2013)研究發現熱帶兩棲類物種類群之間的物種多樣性差異與冠齡有顯著的正相關關系,但是在該研究中也同時發現物種多樣性與凈多樣化速率有顯著的正相關關系,該研究還發現物種多樣化速率與凈第一性生產力和面積存在正相關關系。熱帶地區的科的物種的生態位演化快, 種化速率高,滅絕速率低,古老類群能夠保存,是該熱帶類群物種多樣性高的重要原因。在虎耳草目中,草本和灌木等能夠適應溫帶環境的科有著更高的物種形成速率和凈多樣化速率(Smith & Donoghue,2008),數據分析也顯示這可能是虎耳草目中物種豐富度與凈多樣化速率正相關關系形成的主要原因。Folk等(2019)對虎耳草目生態位祖先重建的結果表明景天科、虎耳草科和茶藨子科等溫帶適應的科的共同祖先,在8千萬年前就已經能夠適應較為干旱和寒冷的環境。在漸新世全球氣候變冷時,這些類群的物種多樣化速率就已經高于虎耳草目木本分支的蕈樹科等。中新世后,隨著全球氣候的進一步變冷,這些溫帶類群的物種多樣化速率快速上升。相比喜濕熱環境的蕈樹科、虎皮楠科等類群,虎耳草目中的草本和灌木類群對寒冷環境的良好適應可能為該類群在全球變冷背景下的快速多樣化(Jian et al.,2008)提供了契機。例如,小二仙草科在始新世—漸新世的全球氣候變冷和干旱化背景下,分布范圍不斷擴張,多樣化速率經歷了大幅度提升(Moody & Garcia, 2021)。

地質歷史事件不僅會重新塑造原有的地形地貌還會改變區域氣候條件,從而對其間的生物演化過程產生影響(唐志堯和方精云,2004; Antonelli et al., 2018)。山脈的隆升直接導致的種群隔離分化以及間接通過生境的復雜化等導致的種群的生殖隔離和豐富的生態位促進了高山植物的快速演化,這也與虎耳草目多個科的快速物種多樣化有著密切的聯系 (Zhang et al., 2014a; 劉杰等, 2017; Moody & Garcia, 2021)。山地草本類群對寒冷生境的快速適應和分化,也是中國青藏高原、北美西海岸山脈和歐洲南部阿爾卑斯山脈成為虎耳草目中高山類群物種多樣性中心的重要原因之一(Zhang et al., 2008; Zhang et al., 2014b; Ebersbach et al., 2017; Xing & Ree, 2017; Stubbs et al., 2020)。隨著青藏高原的隆升,許多植物類群出現了快速多樣化,其中包括虎耳草目中適應寒冷的類群,虎耳草屬(SaxifragaL.)、紅景天屬(RhodiolaL.)等類群的多樣化速率快速升高(Zhang et al., 2014a; Ebersbach et al., 2017; Folk et al., 2019)。因此,虎耳草目的高山類群物種多樣性快速增加,物種多樣性與凈多樣化速率呈現正相關關系。

對類薔薇目(Rosids)的研究發現,包括虎耳草目在內的其他目的多樣化速率也隨著全球變冷而上升,并且多樣化速率隨著緯度升高而升高(Sun et al., 2020)。Igea和Tanentzap(2020)針對被子植物的研究也發現類似的現象,即植物類群在溫帶和寒帶的多樣化速率常常高于熱帶。這些研究的結果表明溫帶類群在全球變冷的過程中物種多樣性在不斷地積累。盡管空間上看,由于高緯度受冰期影響嚴重,物種多樣性較低,但是從類群上看,溫帶類群的物種多樣性的增加速率高于熱帶類群。因此,除了虎耳草目外,在被子植物的科中也可能存在類似的現象,即科間的物種多樣性差異與凈多樣化速率可能存在顯著的正相關關系。

有研究者提出物種多樣性差異很可能由多因素共同導致,因此分支年齡和多樣化速率都可能在其中發揮著作用(Valente et al., 2011; Bloom et al., 2014; Yan et al., 2018)。近年來,也有學者提出該問題的討論應當換一種思考模式,由驗證哪種假說更正確轉向研究這些影響因素之間的互作關系,并建立相關模型來更加數據化與統計學化地看待物種豐富度差異的形成機制(Pontarp et al., 2019)。盡管本研究中發現了凈多樣化速率與年齡的交互作用對物種多樣性差異有很高的解釋率,但是目前該結論是否適用于其他類群,還有待驗證??傊?,迄今為止該問題仍然留有許多值得討論的空間,利用更為復雜的統計模型方法以及更為全面的動植物數據集去探討該問題可能是未來科研工作者們的研究方向。

核基因的同義突變率比葉綠體基因高5倍,比線粒體基因高20倍,且核基因擁有多個獨立位點、雙親遺傳等特征(Small et al., 2004),因此,采用核基因構建的系統發育樹與以往小片段建樹相比,能夠更加全面和客觀地反映分支的系統發育關系。本研究中使用的系統發育樹是Folk等(2019)基于301個核基因構建的虎耳草目主要分支的系統發育樹及分支時間的基礎上,結合葉綠體基因組小片段和ITS序列所構建。本研究結果中各個科的關系與APG IV (APG et al., 2016) 一致。然而,由于物種采樣困難以及大量的核基因測序時間和經費比較多,因此核基因組或者大量核基因建樹的物種覆蓋率都比較低。而取樣率低對物種多樣化過程分析的準確性有著非常大的影響(Chang et al., 2020),因此基因組結合小片段建樹是提高物種覆蓋率的解決方案之一。然而,由于小片段的信息位點不足,可能無法十分準確地估計近緣物種之間的系統發育關系,也可能對系統發育樹的枝長估計產生影響,從而影響到多樣化速率的估計。隨著葉綠體基因組測序價格的下降和分析方法的成熟,增加葉綠體基因組的物種覆蓋率可能是解決近緣物種系統發育關系的較為經濟快捷的方法。

猜你喜歡
物種差異研究
吃光入侵物種真的是解決之道嗎?
英語世界(2023年10期)2023-11-17 09:18:18
相似與差異
音樂探索(2022年2期)2022-05-30 21:01:37
FMS與YBT相關性的實證研究
遼代千人邑研究述論
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
找句子差異
EMA伺服控制系統研究
回首2018,這些新物種值得關注
電咖再造新物種
汽車觀察(2018年10期)2018-11-06 07:05:26
生物為什么會有差異?
主站蜘蛛池模板: 中文国产成人精品久久| 国产成人亚洲精品无码电影| 日韩大片免费观看视频播放| 久久精品人人做人人爽97| 国产呦精品一区二区三区下载| 蝴蝶伊人久久中文娱乐网| 欧美在线国产| 四虎永久在线| 欧美综合中文字幕久久| 欧美一级在线| 91在线中文| 蜜桃视频一区二区三区| 亚洲swag精品自拍一区| 人人91人人澡人人妻人人爽 | 日本伊人色综合网| 91破解版在线亚洲| 91精品国产91欠久久久久| 精品国产电影久久九九| 国内丰满少妇猛烈精品播| 亚洲一区免费看| 成人午夜视频在线| 一级一毛片a级毛片| 欧美区国产区| 五月婷婷欧美| 日韩天堂在线观看| 国产成人综合亚洲欧美在| 亚洲欧洲日韩久久狠狠爱| 一级爱做片免费观看久久 | 一本色道久久88综合日韩精品| 欧美日韩激情| 亚洲第一视频网| 毛片在线看网站| 亚洲成肉网| 亚洲AV免费一区二区三区| 中文字幕在线永久在线视频2020| 亚洲av无码牛牛影视在线二区| 都市激情亚洲综合久久| 一级毛片无毒不卡直接观看| 亚洲日韩欧美在线观看| a级毛片免费看| AV不卡无码免费一区二区三区| 午夜综合网| 欧美成人午夜影院| 国产欧美日韩91| 秋霞一区二区三区| 亚洲色中色| www欧美在线观看| 99在线国产| 国产精品大白天新婚身材| 凹凸国产分类在线观看| 国产精品久久久久久影院| 免费不卡在线观看av| 国产区在线观看视频| 国产精品久久久久久久伊一| 久久婷婷六月| 国产不卡一级毛片视频| 91破解版在线亚洲| 国产精品视频第一专区| 99re在线视频观看| 日韩小视频在线播放| 日韩av电影一区二区三区四区| jizz国产在线| 国产福利拍拍拍| 成人久久精品一区二区三区| 精品久久香蕉国产线看观看gif | 91人妻日韩人妻无码专区精品| av在线5g无码天天| 日韩不卡免费视频| 3p叠罗汉国产精品久久| 美女潮喷出白浆在线观看视频| 思思热在线视频精品| 日韩精品无码免费一区二区三区| 久久香蕉国产线看观| 日韩中文精品亚洲第三区| 久视频免费精品6| 午夜福利网址| 夜夜高潮夜夜爽国产伦精品| 日韩在线2020专区| 国产成人精品男人的天堂下载 | 72种姿势欧美久久久大黄蕉| 国产精品视频a| 99国产精品一区二区|