• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于溫帶和熱帶玉米群體全基因組FST和XP-EHH的 選擇信號檢測

    2019-03-18 06:07:38楊宇昕鄒棖
    中國農(nóng)業(yè)科學(xué) 2019年4期
    關(guān)鍵詞:溫帶熱帶區(qū)段

    楊宇昕,鄒棖

    ?

    基于溫帶和熱帶玉米群體全基因組F和XP-EHH的 選擇信號檢測

    楊宇昕,鄒棖

    (中國農(nóng)業(yè)科學(xué)院作物科學(xué)研究所,北京 100081)

    【目的】玉米起源于熱帶地區(qū),經(jīng)過自然和人工選擇,廣泛的種植于溫帶地區(qū)。開花是玉米生長發(fā)育的中心環(huán)節(jié),也是熱帶玉米向溫帶環(huán)境種植的主要適應(yīng)性性狀。鑒定玉米在馴化過程中出現(xiàn)的受選擇基因區(qū)段,并進(jìn)一步挖掘開花候選基因,為玉米的群體改良、開花遺傳機(jī)理解析提供數(shù)據(jù)支撐?!痉椒ā渴紫葐为?dú)分析30份溫帶玉米自交系和21份熱帶玉米自交系的單倍型數(shù)據(jù),通過過濾高缺失和等位基因頻率較低的變異位點(diǎn),得到高質(zhì)量的SNP(single nucleotide polymorphism)標(biāo)記,利用SnpEff軟件對溫帶和熱帶玉米群體的基因組多態(tài)性位點(diǎn)進(jìn)行了功能預(yù)測。其次過濾得到同時存在于溫帶和熱帶玉米的高質(zhì)量SNP標(biāo)記,對溫帶和熱帶玉米的基因型數(shù)據(jù)進(jìn)行主成分分析(principle component analysis,PCA)以確定其群體結(jié)構(gòu),之后利用群體分化指數(shù)(fixation index,F)和群體間擴(kuò)展單倍型純合度(cross population extended haplotype homozygosity,XP-EHH)法分析溫帶和熱帶玉米群體間的選擇信號分布情況,選擇F和XP-EHH值的top 1%為閾值,篩選得到受選擇位點(diǎn)。通過對SNP進(jìn)行功能注釋得到溫?zé)釒в衩兹后w受到選擇的基因。利用agriGO工具對候選馴化基因進(jìn)行功能富集分析。利用相關(guān)的生物信息學(xué)數(shù)據(jù)庫對候選基因進(jìn)行功能注釋,進(jìn)一步鑒定玉米馴化過程中的開花候選基因?!窘Y(jié)果】通過對溫?zé)釒в衩兹后w的高測序深度的SNP進(jìn)行分析,發(fā)現(xiàn)熱帶玉米群體的SNP數(shù)目為14 123 408個,溫帶玉米群體的SNP數(shù)目為8 791 673個,鑒定到的SNP主要分布于基因間區(qū)。2個群體中均存在的SNP標(biāo)記數(shù)目是204 752個。主成分分析表明溫帶和熱帶玉米可以顯著的分為兩個類群。F選擇信號的top 1%是0.3593,共鑒定到557個候選馴化基因,XP-EHH選擇信號法的top 1%是3.2681,共鑒定到1 913個候選基因。鑒定到多個候選基因與玉米的開花調(diào)控密切相關(guān),包括、、如抑制開花基因的表達(dá),導(dǎo)致玉米在長日照環(huán)境下出現(xiàn)晚花表型,是一個重要的開花調(diào)控基因;COL1與開花促進(jìn)因子FT蛋白互作,加速玉米開花以適應(yīng)長日照環(huán)境;的功能注釋揭示該基因是一個光敏色素互作因子,與光周期基因互作?!窘Y(jié)論】熱帶玉米群體具有更高的遺傳多態(tài)性,篩選到一系列參與了熱帶玉米和溫帶玉米的分化候選基因,并且重點(diǎn)挖掘了參與其中的玉米開花調(diào)控相關(guān)基因。

    玉米;選擇信號;群體分化指數(shù);群體間擴(kuò)展單倍型純合度;開花基因

    0 引言

    【研究意義】玉米是非常重要的糧食和飼料作物,考古學(xué)和遺傳學(xué)研究證明現(xiàn)代玉米起源于約9 000年前的墨西哥西南部,其野生祖先是大芻草[1-3]。玉米的起源中心位于低緯度短日照的熱帶環(huán)境,經(jīng)過不斷地人工選擇和馴化,逐漸降低了光周期敏感性,并且被廣泛種植于58°N到40°S的溫帶長日照地區(qū)[4-5]。根據(jù)玉米對地理環(huán)境的適應(yīng)性,其可以劃分為適應(yīng)長日照環(huán)境的溫帶玉米和適應(yīng)短日照環(huán)境的熱帶玉米2種類型[6]。熱帶玉米不僅具有優(yōu)良的農(nóng)藝性狀(例如根系發(fā)達(dá)、抗病蟲性強(qiáng)[7]等特點(diǎn)),而且具有豐富的遺傳變異[8],是進(jìn)行群體改良的優(yōu)良材料。不同的玉米群體在人工選擇和馴化的過程中,逐漸適應(yīng)當(dāng)?shù)氐纳L環(huán)境,導(dǎo)致染色體區(qū)段上控制目標(biāo)性狀基因的優(yōu)勢等位基因頻率逐漸增加,多態(tài)性發(fā)生改變,根據(jù)中性遺傳假說[9],這將會導(dǎo)致與目的基因緊密連鎖的染色體序列多態(tài)性也會隨之發(fā)生改變,這一個過程稱之為搭車效應(yīng)[10]。由于人工選擇而導(dǎo)致的DNA結(jié)構(gòu)變化稱之為即選擇信號(selection signal),其特點(diǎn)是染色體區(qū)段上出現(xiàn)較長距離的單倍型純合、連鎖不平衡值改變以及優(yōu)勢等位基因頻率的增加。因此,利用選擇信號法探究溫?zé)釒в衩兹后w的基因組變化,可以鑒定在熱帶玉米適應(yīng)性改良過程中受到選擇的基因區(qū)段,并且深入的挖掘開花基因,對于揭示玉米的馴化改良的遺傳機(jī)制具有重要意義。【前人研究進(jìn)展】選擇信號檢測方法根據(jù)其分析原理可以分為以下3種類型:(1)基于等位基因頻率譜的方法,代表的計(jì)算方法包括Tajima's D[11]和CLR(composite likelihood ratio)[12]等;(2)基于連鎖不平衡的方法,主要包括EHH(extended haplotype homozygosity)[13]、iHS(integrated haplotype score)[14]和XP-EHH[15]等;(3)基于群體分化進(jìn)行群體遺傳分化的方法,主要包括F法[16-18]。目前,應(yīng)用較為廣泛的選擇信號檢測方法包括群體分化指數(shù)法(F)和群體間擴(kuò)展單倍型純合度法(XP-EHH)。F統(tǒng)計(jì)理論最早由WRIGHT[16]提出,經(jīng)過對其分析方法和計(jì)算理論的不斷改進(jìn),如今使用較多的是由WEIR等[17]提出的無偏估計(jì)的F。F選擇信號法可以掃描全基因組范圍內(nèi)的SNP位點(diǎn),計(jì)算每個SNP的F值,其取值范圍為0—1,0代表群體間所有位點(diǎn)都沒有出現(xiàn)分化,1代表群體間已經(jīng)完全分化。XP-EHH選擇信號是基于基因的連鎖不平衡原理,群體經(jīng)過人工選擇和改良,將會出現(xiàn)較大范圍的染色體重組,由于連鎖作用的存在,導(dǎo)致突變基因附近的中性位點(diǎn)也會逐代的傳遞,因此,在染色體上形成較長范圍的單倍型純合。XP-EHH的計(jì)算方法是EHH和iHS統(tǒng)計(jì)原理的擴(kuò)展,相比于F方法只能計(jì)算出發(fā)生分化的位點(diǎn),XP-EHH統(tǒng)計(jì)量可以得到選擇作用所在的群體,即XP-EHH統(tǒng)計(jì)值為正數(shù)時,表明選擇發(fā)生在試驗(yàn)群體,反之則發(fā)生在參考群體。隨著越來越多物種的參考基因組構(gòu)建完成,利用F和XP-EHH選擇信號法,已經(jīng)在遺傳進(jìn)化、基因挖掘等生物學(xué)領(lǐng)域取得了很大的進(jìn)展。AXELSSON等[19]對狼和狗采取全基因組重測序得到3.8億個SNP遺傳變異標(biāo)記,通過進(jìn)行F選擇信號檢測,在淀粉消化和脂肪代謝中起到重要作用的10個調(diào)控基因表現(xiàn)出明顯的選擇信號,這些調(diào)控基因可以促進(jìn)狗的祖先在以淀粉為主要食用能量的人類社會中生存,相比于肉食性習(xí)性動物的狼,可以得到更多食物,這種可以消化淀粉的適應(yīng)性演化是狼馴化為狗的關(guān)鍵步驟。LIU等[20]利用溫帶、熱帶和亞熱帶玉米中代表性的260個玉米自交系,結(jié)合F的計(jì)算方法,發(fā)現(xiàn)熱帶玉米相比于溫帶玉米具有更高的遺傳多樣性和更多的等位基因位點(diǎn)。通過研究地方品種和自交系的位點(diǎn)間差異,表明現(xiàn)代玉米低于80%的位點(diǎn)來源于地方品種,為利用玉米地方品種進(jìn)行改良提供了有力的試驗(yàn)支持。HE等[21]利用XP-EHH選擇信號的方法,在溫帶玉米改良進(jìn)程中一共鑒定到超過1 100個候選基因區(qū)段,通過基因富集分析,發(fā)現(xiàn)這些基因主要參與蔗糖合成和油分含量調(diào)控等。【本研究切入點(diǎn)】隨著高通量測序技術(shù)在農(nóng)業(yè)領(lǐng)域的廣泛應(yīng)用,高質(zhì)量的玉米參考基因組圖譜[22-23]和單倍型圖譜[24-25]相繼構(gòu)建完成,極大地促進(jìn)了玉米功能基因組學(xué)的發(fā)展,同時這些公共數(shù)據(jù)也為廣大科研人員進(jìn)行其他的生物學(xué)問題研究提供了寶貴的數(shù)據(jù)資源。熱帶玉米在經(jīng)歷人工馴化和改良過程中逐漸適應(yīng)溫帶環(huán)境,因此,有必要利用選擇信號法揭示在改良過程中發(fā)生選擇的位點(diǎn),挖掘候選馴化基因,進(jìn)一步從基因?qū)用嫔咸骄坑衩椎娜后w改良?!緮M解決的關(guān)鍵問題】本研究選擇具有代表性的30個溫帶玉米自交系和21個熱帶玉米自交系作為材料,基于全基因組單倍型數(shù)據(jù),分析溫?zé)釒в衩谆蚪M多態(tài)性,并且結(jié)合F和XP-EHH 2種方法進(jìn)行群體間選擇信號的檢測,挖掘發(fā)生選擇的馴化基因,為玉米馴化的遺傳機(jī)理解析提供理論依據(jù)。

    1 材料與方法

    1.1 溫、熱帶玉米自交系基因型的獲取

    選取具有代表性的30個溫帶玉米自交系和21個熱帶玉米自交系,其基因型數(shù)據(jù)來自于玉米單倍型圖譜第三版[25],選擇單倍型圖譜中測序深度較高的重測序數(shù)據(jù),以確保其準(zhǔn)確率。原始基因型數(shù)據(jù)以VCF文件存儲??紤]到原始基因型數(shù)據(jù)等位基因頻率偏低可能影響后續(xù)分析[26],因此,利用VCFtools軟件將缺失率較高和最小等位基因頻率低于0.05的SNP剔除,相關(guān)參數(shù)為--maf 0.05。首先進(jìn)行溫帶玉米群體和熱帶玉米群體的SNP篩選,分別得到溫帶和熱帶的基因型數(shù)據(jù),這兩個基因型數(shù)據(jù)只進(jìn)行SNP的功能注釋。使用的選擇信號計(jì)算方法要求溫帶和熱帶玉米群體具有相同的SNP數(shù)目,因此,將溫帶玉米和熱帶玉米一起進(jìn)行SNP的篩選,得到共存于溫帶和熱帶玉米群體的高質(zhì)量SNP標(biāo)記,該基因型用于F和XP-EHH選擇信號的計(jì)算。為了評估溫?zé)釒в衩兹后w全基因組水平上的變異情況,利用SnpEff[27](版本號4.3p)軟件對溫?zé)釒в衩兹后w的變異信息進(jìn)行功能注釋。SNP密度分布利用R軟件包CMplot繪制。玉米參考基因組下載于Ensemble數(shù)據(jù)庫。

    1.2 主成分分析

    主成分分析選擇Tassel(版本號5)軟件。首先將過濾得到的溫?zé)釒в衩兹后w的基因型文件導(dǎo)入到Tassel軟件,按照5個成分進(jìn)行分析,選擇前2個主成分進(jìn)行繪圖展示。

    1.3 FST選擇信號檢測

    采用群體分化指數(shù)(F)法進(jìn)行溫?zé)釒в衩组g選擇信號的檢測,其計(jì)算原理是依據(jù)染色體等位基因頻率變化。群體間選擇信號的檢測可以揭示不同群體馴化過程中經(jīng)歷的自然選擇。按照WEIR等[17]的統(tǒng)計(jì)方法進(jìn)行F值的計(jì)算,考慮到基于單位點(diǎn)SNP掃描的方法容易受到遺傳漂變等因素的影響,因此,為降低假陽性,選擇滑動窗口的計(jì)算方法來增加選擇信號的靈敏度[28]。利用VCFtools軟件計(jì)算滑動窗口內(nèi)的F值,相關(guān)參數(shù)設(shè)置為--fst-window-size 200 kb、--fst-window-step 100 kb。利用R包qqman展示全基因組水平上的F值。為了鑒定F值的受選擇位點(diǎn),選擇F值的top 1%作為顯著閾值線[29],高于閾值線的SNP位點(diǎn)定義為“受選擇位點(diǎn)”,F的計(jì)算公式參考WEIR等[17]報(bào)道。

    1.4 XP-EHH選擇信號計(jì)算

    同時利用基于群體間擴(kuò)展單倍型純合度(XP- EHH)的方法檢測群體間存在分化的基因區(qū)段。XP-EHH選擇信號的計(jì)算使用selscan[30]軟件(版本號v1.1.0),將溫帶玉米群體作為試驗(yàn)群體,熱帶玉米群體作為參考群體。當(dāng)基因組某一區(qū)段的XP-EHH的值為正,代表在試驗(yàn)群體中發(fā)生了選擇,反之則表示參考群體的基因組片段發(fā)生選擇。基于XP-EHH選擇信號得到的統(tǒng)計(jì)值近似符合正態(tài)分布[15],對XP-EHH的值進(jìn)行標(biāo)準(zhǔn)化正態(tài)分布處理,使用selscan的norm參數(shù)對原始的XP-EHH值進(jìn)行標(biāo)準(zhǔn)化處理。將標(biāo)準(zhǔn)正態(tài)化處理后的XP-EHH值從大到小排序,取其top 1%[29]作為顯著的閾值線用以判斷溫、熱帶群體間是否發(fā)生選擇,并且篩選出受選擇基因區(qū)段,XP-EHH的計(jì)算原理參考SABETI等[15]的報(bào)道。

    1.5 候選基因的鑒定

    F法是按照滑動窗口計(jì)算得到的受選擇位點(diǎn),因此,將滑動窗口的起始和終止位點(diǎn)各向上游和下游擴(kuò)增50 kb作為受選擇選擇區(qū)段,使用Bedtools[31](版本號v2.26)軟件將受選擇位點(diǎn)附近的候選馴化基因與玉米參考基因組進(jìn)行基因比對,編碼基因若是落在滑動窗口內(nèi)則定義為候選基因。針對XP-EHH法得到的受選擇位點(diǎn),以顯著SNP位點(diǎn)向上下游各擴(kuò)增50 kb作為受選擇的區(qū)域,同樣利用Bedtools軟件篩選出候選基因。由于篇幅所限,利用maizeGDB(https://www. maizegdb.org/)對選擇信號值較高的基因進(jìn)行基因功能注釋。

    1.6 候選基因的富集分析

    利用在線工具agriGO進(jìn)行候選基因的富集分析[32]。分析方法采用奇異富集分析(singular enrichment analysis,SEA),參考數(shù)據(jù)庫選擇ssp. v5a,顯著性GO條目的檢驗(yàn)使用Fisher精確校驗(yàn),顯著性閾值為0.05。分析內(nèi)容包括細(xì)胞組分、生物過程和分子功能等。

    2 結(jié)果

    2.1 溫、熱帶玉米群體的SNP功能注釋

    利用SnpEff軟件評估溫、熱帶玉米群體在基因組水平上的多態(tài)性情況(圖1),基因組注釋信息來源于B73自交系。結(jié)果表明,熱帶玉米染色體的SNP數(shù)目顯著的高于溫帶群體。熱帶玉米群體一共鑒定到14 123 408個SNP,染色體上平均每145個堿基存在一個變異位點(diǎn);在溫帶玉米群體中共鑒定到個8 791 673個SNP,平均每234個堿基存在一個變異位點(diǎn)。溫帶和熱帶玉米群體一起進(jìn)行SNP篩選,最終得到了204 752個符合過濾標(biāo)準(zhǔn)的SNP位點(diǎn)。此外,還分析了溫帶玉米和熱帶玉米群體中SNP的分布區(qū)域(表1),結(jié)果表明,溫帶和熱帶玉米群體中的SNP變異主要都發(fā)生在基因區(qū)間(Intergentic),此外依次是下游區(qū)間(Downstream)、上游區(qū)間(Upstream)、內(nèi)含子區(qū)間(Intron)、外顯子區(qū)間(Exon)等。

    表1 溫帶和熱帶玉米群體的SNP在基因組區(qū)間的分布比例

    a:溫帶玉米群體SNP分布;b:熱帶玉米群體SNP分布。橫軸代表染色體的物理位置,窗口大小為1 Mb區(qū)間。深綠色代表SNP密度小的區(qū)域,紅色代表SNP密度高的區(qū)域

    2.2 溫、熱帶玉米群體的主成分分析

    溫帶玉米和熱帶玉米由于對光周期敏感性的不同,其表型存在較大差異(圖2),例如株高、葉片數(shù)、開花期等。溫帶材料可以在溫帶環(huán)境正常散粉開花,熱帶材料生殖生長期較長,未能在溫帶環(huán)境散粉開花。此外,熱帶材料的株高也顯著的高于溫帶材料。利用溫帶玉米和熱帶玉米基因組信息,借助主成分分析策略來揭示溫帶和熱帶玉米群體之間的遺傳關(guān)系(圖3),按照PC1和PC2 2個維度可以將使用的玉米自交系分為溫帶和熱帶2個類群。

    2.3 FST值分析結(jié)果

    溫帶玉米光周期敏感性低,可以在長日照的溫帶環(huán)境下正常生長結(jié)實(shí);而熱帶玉米則具有較強(qiáng)的光周期敏感性,主要栽培于熱帶環(huán)境。利用2個群體共有的204 752個SNP標(biāo)記計(jì)算得到每個滑動窗口內(nèi)的F值,全基因組水平上F的top 1%是0.3593(圖4),高于閾值線的受選擇位點(diǎn)一共是1 908個,占總變異位點(diǎn)數(shù)的0.9%。其中第4染色體受到選擇的顯著性位點(diǎn)最多,是414個;第7染色體最少,是61個。第1染色體66 460 001—66 510 000區(qū)間內(nèi)含有最高的F值,其值為0.77。在第1、2、3、4、5、6、8、9、10染色體上都有較強(qiáng)的選擇信號。

    圖2 溫帶(a)和熱帶(b)玉米自交系表型

    圖3 溫、熱帶群體的主成分分析

    2.4 XP-EHH的分析

    XP-EHH選擇信號檢測的原理是根據(jù)原有的突變區(qū)段由于連鎖效應(yīng)較難被打破,而產(chǎn)生新的突變需要較長的選擇周期才能達(dá)到很高的基因頻率。因此,當(dāng)某個單倍型區(qū)段在群體中出現(xiàn)頻率較高,代表該基因區(qū)段經(jīng)歷了選擇,可以利用XP-EHH選擇信號的方法鑒定出來。利用selscan軟件計(jì)算全基因組的上每個SNP位點(diǎn)的XP-EHH值,使用溫帶群體作為試驗(yàn)群體,熱帶群體作為參考群體,當(dāng)XP-EHH值為正值時,代表在溫帶群體發(fā)生了選擇,其值為負(fù)值時,表明馴化選擇發(fā)生在熱帶群體。由于主要關(guān)注熱帶玉米向溫帶玉米的馴化,因此,只分析了XP-EHH值大于0的情況,即溫帶玉米群體中發(fā)生選擇的染色體區(qū)段(圖5),結(jié)果表明,在第2、3、7、8、10染色體存在較多的選擇區(qū)段,其top 1%為3.32,高于閾值線的SNP位點(diǎn)數(shù)目是39 664個,XP-EHH值最高的位點(diǎn)是位于Chr.2:69.324 Mb附近,其值為7.3298。以受選擇SNP位點(diǎn)為核心向上下游各擴(kuò)增50 kb得到100 kb的候選基因區(qū)段,利用Bedtools軟件將其與玉米B73參考基因組進(jìn)行比對,共得到1 913個候選基因,其中包含2個與玉米適應(yīng)性改良關(guān)系密切相關(guān)的基因,分別是[33]和[34]。(Chr.9:115 786 897—115 789 787),其F值為3.58,最新的一項(xiàng)研究表明是玉米開花控通路中的一個重要調(diào)節(jié)因子[33],在長日照環(huán)境下,抑制開花促進(jìn)因子[35-36]的表達(dá),導(dǎo)致玉米晚花,通過CRISPR/Cas9技術(shù)介導(dǎo)的敲除,發(fā)現(xiàn)可以促進(jìn)玉米提早開花。鑒定到在溫帶群體受到選擇,因而加速熱帶玉米群體對溫帶環(huán)境的適應(yīng),這與前人報(bào)道一致[33]。(Chr.9:108 447 974—108 449 794)是玉米光周期調(diào)控通路中的重要基因,該基因編碼蛋白可以和開花促進(jìn)因子FT蛋白互作,加速玉米開花以適應(yīng)長日照環(huán)境,是一個重要的開花調(diào)控基因。

    2.5 候選基因的富集分析和功能注釋

    利用F鑒定得到557個候選基因和XP-EHH得到的1 913個候選基因?;蚋患治霰砻?i>F法得到的候選基因富集在生物學(xué)過程(biological process,P),共19個顯著的GO條目(表2)。這些GO條目主要涉及到肌動蛋白發(fā)育調(diào)控、微絲的發(fā)育調(diào)控和細(xì)胞骨架構(gòu)成等。利用XP-EHH得到的1 913個基因沒有富集得到顯著的GO條目,推測是因?yàn)闊釒в衩遵Z化溫帶玉米涉及到眾多農(nóng)藝性狀的改變,并且許多農(nóng)藝性狀是由一系列微效數(shù)量位點(diǎn)控制,因此未能得到顯著的GO條目。盡管如此,XP-EHH的方法仍然可以得到較多的選擇基因,并且這些基因中鑒定得到了與玉米花期調(diào)控相關(guān)的基因,例如和。因此,通過選擇信號的方法可以為解析玉米馴化的遺傳機(jī)理、鑒定候選馴化基因提供實(shí)驗(yàn)依據(jù)。除了這些已經(jīng)報(bào)道過的調(diào)控基因,還統(tǒng)計(jì)了選擇信號值高于top 1%和top 0.1%的候選基因(表3)。另外,通過maizeGDB網(wǎng)站進(jìn)行基因功能注釋,發(fā)現(xiàn)這些基因均是玉米生長發(fā)育過程中的重要調(diào)控因子,例如是一個光敏色素調(diào)節(jié)因子,參與玉米開花通路調(diào)控;參與下胚軸的發(fā)育調(diào)控。

    表2 FST選擇信號得到的基因進(jìn)行GO富集分析的結(jié)果

    表3 選擇信號鑒定得到的候選基因功能注釋

    a受選擇區(qū)域內(nèi)較高F值的基因;b受選擇區(qū)域內(nèi)較高XP-EHH值的基因;c基于maizeGDB的基因功能注釋;d同時受到F和XP-EHH選擇的基因;*代表受到選擇的關(guān)鍵基因(選擇信號值高于top 1%);**代表受到強(qiáng)烈選擇的基因(選擇信號值高于 top 0.1%)

    aGenes with higherFvalues in selected regions;bGenes with higher XP-EHH values in selected regions;cGene function annotation based on maizeGDB;dCandidate genes identified by bothFand XP-EHH methods; *Represents the selected key genes (the selection signal over top 1%); **Represents the strongly selected genes (the selection signal over top 0.1%)

    閾值線代表選擇信號值的top 1% the cutoff lines represent the top 1% of the XP-EHH selection signal

    3 討論

    3.1 溫?zé)釒в衩谆蚪M變異分析

    玉米起源于熱帶環(huán)境的墨西哥西南部,經(jīng)過對其進(jìn)行適應(yīng)性改良,玉米逐漸的擴(kuò)散到世界各地。根據(jù)光周期敏感性的不同,玉米可以被劃分為溫帶玉米和熱帶玉米[6]。熱帶玉米靠近玉米起源中心,很少經(jīng)歷人工選擇和改良,保留著大量的有利變異,是進(jìn)行群體遺傳改良的重要材料[8]。本研究選擇30份溫帶玉米自交系和21份熱帶玉米作為研究材料進(jìn)行重測序。通過SnpEff軟件對變異文件進(jìn)行基因組的功能的注釋和預(yù)測,發(fā)現(xiàn)熱帶玉米基因組具有更多的SNP(圖1),無論是SNP的總數(shù)還是平均分布密度均顯著高于溫帶玉米。該結(jié)果證明了熱帶玉米具有更高的遺傳多樣性,這與前人研究即熱帶玉米具有更高的遺傳多樣性一致[8]。主成分分析發(fā)現(xiàn)溫帶玉米群存在更大的變異幅度(圖2),推測是因?yàn)闇貛в衩兹后w包含較多的亞群,因此其變異幅度更廣,該結(jié)果和前人的主成分分析一致[6]。這也表明本研究所使用的基因型數(shù)據(jù)以及群體劃分的準(zhǔn)確性,為后續(xù)的選擇信號分析奠定了數(shù)據(jù)基礎(chǔ)。玉米作為一個廣適性的作物,其在基因組水平存在多態(tài)性變化,進(jìn)而導(dǎo)致相關(guān)表型變化,是使其成為有重要影響力作物的一個重要條件。研究揭示了溫?zé)釒в衩谆蚪M水平SNP變異情況,可為后續(xù)的玉米功能基因組學(xué)研究提供理論基礎(chǔ)。

    3.2 FST和XP-EHH的計(jì)算原理

    隨著海量生物學(xué)數(shù)據(jù)的產(chǎn)出,使用高密度的SNP標(biāo)記進(jìn)行群體間選擇信號的分析在很多物種中都已成功進(jìn)行,包括人類[37]、綿羊[38]、水稻[39]等。選擇信號是一種重要的群體遺傳學(xué)研究手段,利用合適的選擇信號可以顯著的揭示玉米適應(yīng)性進(jìn)化的遺傳機(jī)理,并且可以深入挖掘與農(nóng)藝性狀相關(guān)的候選基因[21]。本研究利用F和XP-EHH的方法進(jìn)行選擇信號的鑒定。F可以利用等位基因的頻率變化來篩選基因組受到馴化選擇的基因區(qū)段[16],這是因?yàn)橥粋€物種由于馴化的目的、種植環(huán)境、改良時間不同將會導(dǎo)致基因組上的等位基因頻率發(fā)生變化,雖然它可以鑒定出受到選擇的基因,但是不能確定選擇的方向[40],因此,本研究進(jìn)一步引入XP-EHH選擇信號方法。XP-EHH是一種基于基因組上單倍型純合度的方法,鑒定不同群體間發(fā)生選擇作用的區(qū)段,XP-EHH統(tǒng)計(jì)分析方法引入了參考群體和試驗(yàn)群體的概念,因此可以根據(jù)統(tǒng)計(jì)值判斷基因組發(fā)生選擇的方向,這為深入篩選候選馴化基因提供了更有力的數(shù)據(jù)支撐[41]。

    3.3 選擇信號檢測結(jié)果的評價(jià)

    傳統(tǒng)的選擇信號檢測方法主要利用RFLP、AFLP和SSR等分子標(biāo)記。有研究利用玉米和大芻草的99個SSR標(biāo)記作為選擇信號進(jìn)行分析[2],揭示了現(xiàn)代玉米起源于約9 000年前的墨西哥西南部。然而RFLP、AFLP和SSR等分子標(biāo)記開發(fā)程序繁瑣,并且定位精度不高,如今已經(jīng)很少應(yīng)用于選擇信號的研究之中。隨著新一代測序技術(shù)的發(fā)展,基于高通量測序技術(shù)得到的SNP標(biāo)記由于其精度高、密度大以及定位更加準(zhǔn)確等優(yōu)點(diǎn),逐漸取代了傳統(tǒng)的分子標(biāo)記,成為進(jìn)行選擇信號檢測的主要方法[19-21]。本研究選擇具有代表性的溫帶和熱帶玉米群體的重測序數(shù)據(jù),進(jìn)行全基因組水平的F和XP-EHH的2種選擇信號的鑒定,分別篩選到了557和1 913個候選基因。利用F鑒定得到的基因進(jìn)行富集分析,結(jié)果表明這些基因富集在一些重要的生物學(xué)通路上,主要參與肌動蛋白的發(fā)育調(diào)節(jié)(GO:0008064,GO:0032956)、微絲的發(fā)育調(diào)節(jié)(GO:0030832)、細(xì)胞骨架發(fā)育調(diào)控(GO:0030036)等。本研究利用F得到基因進(jìn)行富集分析,得到的顯著GO條目較少,此外XP-EHH選擇信號法鑒定到的候選基因未能得到顯著的GO條目,不能準(zhǔn)確地反映群體分化的結(jié)果。推測是目前所用的參考數(shù)據(jù)庫信息不夠全面,因此未能特異的鑒定出富有生物學(xué)意義的GO條目。本研究通過對具有較高選擇信號的基因進(jìn)行功能注釋,發(fā)現(xiàn)這些基因是在熱帶玉米改良為溫帶玉米過程中受到了極大的選擇,同時也是玉米生長發(fā)育過程中重要的調(diào)控因子(表3)。因此,結(jié)合相關(guān)的表型性狀和選擇信號法可以鑒定玉米適應(yīng)改良過程中發(fā)生變化的基因。其中通過XP-EHH選擇信號法,鑒定到溫?zé)釒в衩走m應(yīng)性改良相關(guān)的調(diào)控基因[33]和玉米開花調(diào)控相關(guān)的基因[34]。盡管前人曾經(jīng)利用過相關(guān)的選擇信號鑒定了玉米在馴化過程中的選擇基因[21],然而本研究著重挖掘了與開花相關(guān)的調(diào)控基因,開花是玉米生育期的一個重要農(nóng)藝性狀,關(guān)系到植株從營養(yǎng)生長到生殖生長的轉(zhuǎn)變。其中包括已經(jīng)報(bào)道過的和。此外在候選基因集中鑒定到了一些開花調(diào)控候選基因,例如,基因功能注釋表明該基因是一個光敏色素互作因子(表3),與光周期基因互作[42]。因此,本研究鑒定得到的候選基因可以為后續(xù)的開花遺傳機(jī)理解析提供理論基礎(chǔ)。這些研究表明結(jié)合F和XP-EHH選擇信號法是進(jìn)行群體遺傳分化、功能富集分析和馴化基因挖掘的一個重要手段。

    4 結(jié)論

    基于對溫帶和熱帶玉米的基因型數(shù)據(jù)分析,發(fā)現(xiàn)熱帶玉米具有更高的遺傳多樣性,表明其是玉米分子育種的重要資源。鑒定到玉米馴化改良過程中受到選擇的基因,且部分受選擇的基因參與玉米開花發(fā)育調(diào)控途徑。

    [1] PIPERNO D R, RANERE A J, HOLST I, IRIARTE J, DICKAU R. Starch grain and phytolith evidence for early ninth millennium BP maize from the Central Balsas River Valley, Mexico., 2009, 106(13): 5019-5024.

    [2] MATSUOKA Y, VIGOUROUX Y, GOODMAN M M, SANCHEZ J, BUCKLER E, DOEBLEY J. A single domestication for maize shown by multilocus microsatellite genotyping., 2002, 99(9): 6080-6084.

    [3] VAN HEERWAARDEN J, DOEBLEY J, BRIGGS W H, GLAUBITZ J C, GOODMAN M M, Gonzalez J d J S, ROSS-IBARRA J. Genetic signals of origin, spread, and introgression in a large sample of maize landraces., 2011, 108(3): 1088-1092.

    [4] BUCKLER E S, HOLLAND J B, BRADBURY P J, ACHARYA C B, BROWN P J, BROWNE C, ERSOZ E, FLINT-GARCIA S, GARCIA A, GLAUBITZ J C, et at.. The genetic architecture of maize flowering time., 2009, 325(5941): 714-718.

    [5] SWARTS K, GUTAKER R M, BENZ B, BLAKE M, BUKOWSKI R, HOLLAND J, KRUSE-PEEPLES M, LEPAK N, PRIM L, ROMAY M C, et at. Genomic estimation of complex traits reveals ancient maize adaptation to temperate North America., 2017, 357(6350): 512-515.

    [6] LU Y L, YAN J B, GUIMARAES C T, TABA S, HAO Z F, GAO S B, CHEN S J, LI J S, ZHANG S H, VIVEK B S, et at.. Molecular characterization of global maize breeding germplasm based on genome-wide single nucleotide polymorphisms., 2009, 120(1): 93-115.

    [7] TALLURY S, GOODMAN M. Experimental evaluation of the potential of tropical germplasm for temperate maize improvement., 1999, 98(1): 54-61.

    [8] HALLAUER A R, CARENA M J. Adaptation of tropical maize germplasm to temperate environments., 2013, 196(1): 1-11.

    [9] KIMURA M. Evolutionary rate at the molecular level., 1968, 217(5129): 624-626.

    [10] SMITH J M, HAIGH J. The hitch-hiking effect of a favourable gene., 1974, 23(1): 23-35.

    [11] TAJIMA F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism., 1989, 123(3): 585-595.

    [12] NIELSEN R, WILLIAMSON S, KIM Y, HUBISZ M J, CLARK A G, BUSTAMANTE C. Genomic scans for selective sweeps using SNP data., 2005, 15(11): 1566-1575.

    [13] SABETI P C, REICH D E, HIGGINS J M, LEVINE H Z, RICHTER D J, SCHAFFNER S F, GABRIEL S B, PLATKO J V, PATTERSON N J, MCDONALD G J, et at.. Detecting recent positive selection in the human genome from haplotype structure., 2002, 419(6909): 832.

    [14] VOIGHT B F, KUDARAVALLI S, WEN X Q, PRITCHARD J K. A map of recent positive selection in the human genome., 2006, 4(3): e72.

    [15] SABETI P C, VARILLY P, FRY B, LOHMUELLER J, HOSTETTER E, COTSAPAS C, XIE X, BYRNE E H, MCCARROLL S A, et at.. Genome-wide detection and characterization of positive selection in human populations., 2007, 449(7164): 913.

    [16] WRIGHT S. The genetical structure of populations., 1949, 15(1): 323-354.

    [17] WEIR B S, COCKERHAM C C. Estimating F-statistics for the analysis of population structure., 1984, 38(6): 1358-1370.

    [18] GIANOLA D, SIMIANER H, QANBARI S. A two-step method for detecting selection signatures using genetic markers., 2010, 92(2): 141-155.

    [19] AXELSSON E, RATNAKUMAR A, ARENDT M L, MAQBOOL K, WEBSTER M T, Perloski M, Liberg O, ARNEMO J M, HEDHAMMAR A, LINDBLAD-TOH K. The genomic signature of dog domestication reveals adaptation to a starch-rich diet., 2013, 495(7441): 360.

    [20] LIU K J, GOODMAN M, MUSE S, SMITH J S, BUCKLER E, DOEBLEY J. Genetic structure and diversity among maize inbred lines as inferred from DNA microsatellites., 2003, 165(4): 2117-2128.

    [21] HE C, FU J J, ZHANG J, Li Y X, ZHENG J, ZHANG H W, YANG X H, WANG J H, WANG G Y. A gene-oriented haplotype comparison reveals recently selected genomic regions in temperate and tropical maize germplasm., 2017, 12(1): e0169806.

    [22] SCHNABLE P S, WARE D, FULTON R S, STEIN J C, WEI F S, PASTERNAK S, LIANG C Z, ZHANG J W, FULTON L, GRAVES T A, et at.. The B73 maize genome: complexity, diversity, and dynamics., 2009, 326(5956): 1112-1115.

    [23] JIAO Y P, PELUSO P, SHI J H, LIANG T, STITZER M C, WANG B, CAMPBELL M S, STEIN J C, WEI X H, CHIN C S, et at.. Improved maize reference genome with single-molecule technologies., 2017, 546(7659): 524.

    [24] CHIA J M, SONG C, BRADBURY P J, COSTICH D, DE LEON N, DOEBLEY J, ELSHIRE R J, GAUT B, GELLER L, GLAUBITZ J C, et at.. Maize HapMap2 identifies extant variation from a genome in flux., 2012, 44(7): 803.

    [25] BUKOWSKI R, GUO X S, LU Y L, ZOU C, HE B, RONG Z Q, WANG B, XU D W, YANG B C, XIE C X, et at.. Construction of the third-generation Zea mays haplotype map., 2017, 7(4): gix134.

    [26] DANECEK P, AUTON A, ABECASIS G, ALBERS C A, BANKS E, DEPRISTO M A, HANDSAKER R E, LUNTER G, MARTH G T, SHERRY S T, et at.. 1000 Genomes Project Analysis Group. The variant call format and VCFtools., 2011, 27(15): 2156-2158.

    [27] CINGOLANI P, PLATTS A, WANG L L, COON M, NGUYEN T, WANG L A, LAND S J, LU X Y, RUDEN D M. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of drosophila melanogaster strain w1118; iso-2; iso-3., 2012, 6(2): 80-92.

    [28] MA Y, DING X, QANBARI S, WEIGEND S, ZHANG Q, SIMIANER H. Properties of different selection signature statistics and a new strategy for combining them., 2015, 115(5): 426.

    [29] HOAGLIN D C, MOSTELLER F, TUKEY J W.. New York: John Wiley & Sons, 1983.

    [30] SZPIECH Z A, HERNANDEZ R D. selscan: an efficient multithreaded program to perform EHH-based scans for positive selection., 2014, 31(10): 2824-2827.

    [31] QUINLAN A R, HALL I M. BEDTools: a flexible suite of utilities for comparing genomic features., 2010, 26(6): 841-842.

    [32] TIAN T, LIU Y, YAN H Y, You Q, Yi X, Du Z, XU W Y, SU Z. agriGO v2. 0: a GO analysis toolkit for the agricultural community, 2017 update., 2017, 45(W1): W122-W129.

    [33] HUANG C, SUN H Y, XU D Y, LIANG Y M, WANG X F, XU G H, TIAN J G, WANG C L, LI D, WU L H, et at..enhances maize adaptation to higher latitudes., 2018, 115(2): E334-E341.

    [34] KHAN S, ROWE S C, HARMON F G. Coordination of the maize transcriptome by a conserved circadian clock., 2010, 10(1): 126.

    [35] GUO L, WANG X, ZHAO M, HUANG C, LI C, LI D, YANG C , YORK A M, XUE W, XU G, LIANG Y, CHEN Q, DOEBLEY J F, TIAN F. Stepwise cis-regulatory changes incontribute to maize flowering-time adaptation., 2018, 28(18): 3005-3015.

    [36] MENG X, MUSZYNSKI M G, DANILEVSKAYA O N. The-likegene functions as a floral activator and is involved in photoperiod sensitivity in maize., 2011, 23(3): 942-960.

    [37] SHEEHAN M J, NACHMAN M W. Morphological and population genomic evidence that human faces have evolved to signal individual identity., 2014, 5: 4800.

    [38] 曾滔, 趙福平, 王光凱, 吳明明, 魏彩虹, 張莉, 李利, 張紅平, 杜立新. 基于群體分化指數(shù)FST的綿羊全基因組選擇信號檢測. 畜牧獸醫(yī)學(xué)報(bào), 2013, 44(12): 1891-1899.

    ZENG T, ZHAO F P, WANG G K, WU M M, WEI C H, ZHANG L, LI L, ZHANG H P, DU L X. Genome-wide detection of selection signatures in sheep populations with use of population differentiation index FST., 2013, 44(12): 1891-1899. (in Chinese)

    [39] HUANG X, SANG T, ZHAO Q, QI F, ZHAO Y, LI C Y ZHU C R, LU T T, ZHANG Z W, LI M, et at.. Genome-wide association studies of 14 agronomic traits in rice landraces., 2010, 42(11): 961.

    [40] MCVICKER G, GORDON D, DAVIS C, GREEN P. Widespread genomic signatures of natural selection in hominid evolution., 2009, 5(5): e1000471.

    [41] 薛周艤源, 宋顯威, 吳林慧, 王露珍, 崔家安, 孫章健, 張政, 馬云龍. 畜禽選擇信號檢測方法及其統(tǒng)計(jì)學(xué)問題. 畜牧獸醫(yī)學(xué)報(bào), 2018, 49(6): 1099-1107.

    XUE Z Y Y, SONG X W, WU L H, WANG L Z, CUI J A, SUN Z J, ZHANG Z, MA Y L. The identification methods of selection signatures in livestock and its statistical problems., 2018, 49(6): 1099-1107. (in Chinese)

    [42] KUMAR I, SWAMINATHAN K, HUDSON K, HUDSON M E. Evolutionary divergence of phytochrome protein function inPIF3 signaling., 2016, 67(14): 4231-4240.

    (責(zé)任編輯 李莉)

    Genome-Wide Detection of Selection Signal in Temperate and Tropical Maize Populations with Use ofFand XP-EHH

    YANG YuXin, ZOU Cheng

    (Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing 100081)

    【Objective】Maize was first domesticated in tropical areas, but it has been cultivated widely in the temperate regions after natural and artificial selection. Flowering time is not only the key component of the entire growth period, but also a major adaptive trait during the dispersal process from tropical to temperate conditions. Thus, identifying the selected gene regions responsible for the adaptation to temperate zones, and discovering the genes that are involved in flowering time could provide a molecular basis for improving maize and for dissecting its flowering mechanism. 【Method】We analyzed the haplotype data of 30 temperate and 21 tropical maize inbred lines. High quality SNP (single nucleotide polymorphism) markers were obtained after filtering out SNPs with high missing rates and low allele frequencies. These high quality SNPs were annotated by SNPeff. Principle component analysis (PCA) of the genotypic data of temperate and tropical maize was performed to further validate the population structure of these samples. Using high quality SNP markers that were present in tropical and temperate populations, we calculated the selection signal using the fixation index (F) and cross population extended haplotype homozygosity (XP-EHH) methods. The top 1% of values was used as a significant threshold to identify the candidate selected signals. The candidate selected genes that we selected from temperate and tropical maize were identified based on their SNP annotation. The function of these selected genes was characterized furtherly by the GO enrichment analysis using agriGO. To identify the genes for flowering time that were under selection, bioinformatics databases were examined that contained relevant data on maize. 【Result】By analyzing the high depth resequencing data, we found 14123408 and 8791673 SNPs in tropical and temperate populations, respectively. The identified SNPs were mainly distributed in the intergenic regions. There were 204752 high quality SNPs that coexisted in temperate and tropical populations. PCA indicated that temperate and tropical maize can be divided into two groups. The top 1% ofFvalue and XP-EHH were 0.3059, 3.2681, and a total of 557 and 1 913 candidate genes were identified byFand XP-EHH methods, respectively. Many candidate genes were highly related to regulation of flowering time, which included,and.is a vital gene for regulating flowering time, and it negatively regulated the floral activator gene, which cause the late flowering time phenotype under long-day conditions. COL1 positively interacts with the FT protein to promote the transition of flowering time to adapt to the long-day environment. Functional annotations ofrevealed that it was a phytochrome interacting factor, and interacts with photoperiod gene. 【Conclusion】Our study revealed that tropical maize had higher genetic diversity than temperate maize. A series of genes that were under selection during the adaptation to tropical to temperate conditions were predicted, and we further explored the genes that were involved in flowering during this process.

    maize; selection signal; fixation index; cross population extended haplotype homozygosity; flowering time genes

    2018-10-30;

    2018-12-09

    國家重點(diǎn)研發(fā)計(jì)劃(2016YFD0100303)、國家自然科學(xué)基金面上項(xiàng)目(31371638)

    楊宇昕,E-mail:yyx0719@126.com。通信作者鄒棖,E-mail:zoucheng@caas.cn

    10.3864/j.issn.0578-1752.2019.04.001

    猜你喜歡
    溫帶熱帶區(qū)段
    第12期 參考答案
    第31期 參考答案
    Facts of Yellowstone
    中老鐵路雙線區(qū)段送電成功
    熱帶風(fēng)情
    女報(bào)(2020年7期)2020-08-17 07:16:05
    熱帶的鳥兒
    站內(nèi)特殊區(qū)段電碼化設(shè)計(jì)
    站內(nèi)軌道區(qū)段最小長度的探討
    淺析分路不良區(qū)段解鎖的特殊操作
    圓滾滾的熱帶“龍”
    老司机午夜十八禁免费视频| 精品人妻1区二区| 99国产精品一区二区三区| 欧美国产精品va在线观看不卡| 人人澡人人妻人| 久久精品国产99精品国产亚洲性色| 亚洲成a人片在线一区二区| 最新在线观看一区二区三区| 无人区码免费观看不卡| av片东京热男人的天堂| 国产极品粉嫩免费观看在线| 欧美成狂野欧美在线观看| 在线观看免费日韩欧美大片| 免费看十八禁软件| 99在线视频只有这里精品首页| 国产熟女午夜一区二区三区| 亚洲成人久久爱视频| 免费看日本二区| 久久午夜综合久久蜜桃| 日日摸夜夜添夜夜添小说| 欧美日本视频| 亚洲最大成人中文| 女人被狂操c到高潮| 我的亚洲天堂| www.精华液| 亚洲男人的天堂狠狠| 成熟少妇高潮喷水视频| 狂野欧美激情性xxxx| 精品欧美一区二区三区在线| 亚洲成国产人片在线观看| √禁漫天堂资源中文www| 日韩欧美在线二视频| 在线播放国产精品三级| 怎么达到女性高潮| 99国产精品99久久久久| 老熟妇乱子伦视频在线观看| 亚洲欧美精品综合一区二区三区| 久久久久久免费高清国产稀缺| 亚洲 国产 在线| 12—13女人毛片做爰片一| 国产在线观看jvid| 少妇被粗大的猛进出69影院| 啦啦啦 在线观看视频| 亚洲性夜色夜夜综合| 国产精品久久久av美女十八| 亚洲国产看品久久| 十八禁网站免费在线| 国产精品一区二区三区四区久久 | 欧美亚洲日本最大视频资源| 男人舔奶头视频| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲欧美日韩无卡精品| 欧美黑人精品巨大| 久久99热这里只有精品18| 国产主播在线观看一区二区| 美女高潮喷水抽搐中文字幕| 成人亚洲精品一区在线观看| 亚洲av第一区精品v没综合| 亚洲精华国产精华精| 亚洲五月婷婷丁香| 国产真实乱freesex| 欧美成人免费av一区二区三区| 精品久久久久久久久久免费视频| 麻豆一二三区av精品| 一二三四在线观看免费中文在| 嫩草影视91久久| 国产一区二区在线av高清观看| 久9热在线精品视频| 怎么达到女性高潮| 亚洲真实伦在线观看| av在线天堂中文字幕| 欧美丝袜亚洲另类 | 精品午夜福利视频在线观看一区| 激情在线观看视频在线高清| 男女下面进入的视频免费午夜 | 好看av亚洲va欧美ⅴa在| 午夜免费成人在线视频| 亚洲aⅴ乱码一区二区在线播放 | 精品国产超薄肉色丝袜足j| 亚洲,欧美精品.| 嫁个100分男人电影在线观看| 中文字幕最新亚洲高清| 国内少妇人妻偷人精品xxx网站 | 久久久精品欧美日韩精品| 很黄的视频免费| 国产亚洲精品综合一区在线观看 | 日韩视频一区二区在线观看| 婷婷丁香在线五月| 国产视频一区二区在线看| 国产精品免费一区二区三区在线| 午夜免费激情av| 精品久久久久久久久久免费视频| 999精品在线视频| 国产一区二区在线av高清观看| 国产国语露脸激情在线看| 一区二区三区国产精品乱码| 精华霜和精华液先用哪个| 国产成人精品久久二区二区91| 欧美日韩亚洲国产一区二区在线观看| 亚洲最大成人中文| 婷婷六月久久综合丁香| 欧美成狂野欧美在线观看| 欧美最黄视频在线播放免费| 亚洲美女黄片视频| 亚洲国产欧美一区二区综合| 国产人伦9x9x在线观看| 国产精品综合久久久久久久免费| 成人亚洲精品一区在线观看| e午夜精品久久久久久久| 久久精品亚洲精品国产色婷小说| 日韩欧美一区二区三区在线观看| 99精品久久久久人妻精品| 精品高清国产在线一区| 久热这里只有精品99| 在线观看午夜福利视频| 久久午夜综合久久蜜桃| 免费看a级黄色片| 亚洲av第一区精品v没综合| 99热这里只有精品一区 | 精品第一国产精品| 久久欧美精品欧美久久欧美| 欧美zozozo另类| 欧美大码av| 不卡一级毛片| 日韩欧美免费精品| xxx96com| 亚洲精品一区av在线观看| 51午夜福利影视在线观看| 亚洲国产中文字幕在线视频| 中文亚洲av片在线观看爽| 免费高清视频大片| 国产av又大| 国产亚洲欧美精品永久| 动漫黄色视频在线观看| 国产免费av片在线观看野外av| 国产精品久久久久久人妻精品电影| √禁漫天堂资源中文www| 在线观看www视频免费| av有码第一页| 精品久久久久久成人av| 久久久久久九九精品二区国产 | 桃色一区二区三区在线观看| 欧美日韩精品网址| 琪琪午夜伦伦电影理论片6080| 国产日本99.免费观看| 两性午夜刺激爽爽歪歪视频在线观看 | 手机成人av网站| 亚洲成人国产一区在线观看| xxxwww97欧美| 国产视频一区二区在线看| 啦啦啦 在线观看视频| 桃色一区二区三区在线观看| 在线观看www视频免费| videosex国产| 久久伊人香网站| 国产一区二区三区视频了| 19禁男女啪啪无遮挡网站| 97超级碰碰碰精品色视频在线观看| 757午夜福利合集在线观看| 他把我摸到了高潮在线观看| 亚洲欧美激情综合另类| 亚洲 欧美 日韩 在线 免费| 免费观看人在逋| 国产久久久一区二区三区| 黄网站色视频无遮挡免费观看| 久久精品国产亚洲av香蕉五月| 欧美成狂野欧美在线观看| 在线国产一区二区在线| 久久精品夜夜夜夜夜久久蜜豆 | 国产视频内射| 中文亚洲av片在线观看爽| 久久久久久久午夜电影| 日韩欧美在线二视频| 每晚都被弄得嗷嗷叫到高潮| 欧美在线一区亚洲| 免费一级毛片在线播放高清视频| 他把我摸到了高潮在线观看| 亚洲免费av在线视频| 又紧又爽又黄一区二区| 天天添夜夜摸| 久久热在线av| 中文字幕人成人乱码亚洲影| av福利片在线| 亚洲专区国产一区二区| 免费看美女性在线毛片视频| 怎么达到女性高潮| 午夜福利一区二区在线看| 麻豆成人午夜福利视频| 777久久人妻少妇嫩草av网站| АⅤ资源中文在线天堂| 日本三级黄在线观看| 国产高清视频在线播放一区| 久久亚洲真实| 女人爽到高潮嗷嗷叫在线视频| 在线永久观看黄色视频| 女警被强在线播放| 久久香蕉精品热| 国产精品爽爽va在线观看网站 | 欧美日本亚洲视频在线播放| 欧美成狂野欧美在线观看| av电影中文网址| 免费电影在线观看免费观看| 老司机午夜福利在线观看视频| 999久久久精品免费观看国产| 欧美性猛交╳xxx乱大交人| 成人18禁在线播放| 日韩成人在线观看一区二区三区| 久99久视频精品免费| 国产精品国产高清国产av| 国产av不卡久久| 成人三级做爰电影| 国产三级在线视频| 国产真实乱freesex| 757午夜福利合集在线观看| 久久久水蜜桃国产精品网| www国产在线视频色| 久久精品91蜜桃| 久久人妻福利社区极品人妻图片| 亚洲天堂国产精品一区在线| 精品久久久久久,| 制服诱惑二区| 国产精品1区2区在线观看.| 成人免费观看视频高清| 成人18禁高潮啪啪吃奶动态图| 欧美又色又爽又黄视频| 一边摸一边抽搐一进一小说| 欧美 亚洲 国产 日韩一| 男女床上黄色一级片免费看| www.精华液| 日韩精品青青久久久久久| 色精品久久人妻99蜜桃| 一本大道久久a久久精品| 日本精品一区二区三区蜜桃| 日韩欧美在线二视频| 嫩草影视91久久| 亚洲精品久久国产高清桃花| 欧美性长视频在线观看| 亚洲av日韩精品久久久久久密| 成人18禁在线播放| 国产真人三级小视频在线观看| 中文字幕人妻熟女乱码| 精品不卡国产一区二区三区| 大型av网站在线播放| 国产激情欧美一区二区| 国产91精品成人一区二区三区| 久久久久久亚洲精品国产蜜桃av| 日韩精品青青久久久久久| 久久天躁狠狠躁夜夜2o2o| 色综合婷婷激情| 香蕉av资源在线| 国产97色在线日韩免费| 精品电影一区二区在线| 国产色视频综合| 欧美日韩瑟瑟在线播放| 久久中文字幕人妻熟女| 波多野结衣巨乳人妻| 香蕉av资源在线| 免费看日本二区| 国产精品综合久久久久久久免费| 国产亚洲av嫩草精品影院| 成人精品一区二区免费| 久久精品人妻少妇| 久久香蕉国产精品| 久久草成人影院| 成人av一区二区三区在线看| 大香蕉久久成人网| 久热爱精品视频在线9| 免费看美女性在线毛片视频| 精品国产国语对白av| 在线免费观看的www视频| 可以免费在线观看a视频的电影网站| 999精品在线视频| 中亚洲国语对白在线视频| 女人爽到高潮嗷嗷叫在线视频| www.999成人在线观看| 国产亚洲欧美精品永久| 丁香欧美五月| 亚洲欧美激情综合另类| 嫩草影视91久久| 国产亚洲av嫩草精品影院| 欧美日韩精品网址| 亚洲av熟女| 99久久99久久久精品蜜桃| 亚洲国产欧洲综合997久久, | 天堂√8在线中文| 久久伊人香网站| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲成a人片在线一区二区| 黄色 视频免费看| 久久天躁狠狠躁夜夜2o2o| 日本熟妇午夜| 男女之事视频高清在线观看| 黑人操中国人逼视频| 久久精品影院6| 亚洲精品国产一区二区精华液| 神马国产精品三级电影在线观看 | 看免费av毛片| www.999成人在线观看| www.熟女人妻精品国产| 一区福利在线观看| 人妻久久中文字幕网| 午夜免费成人在线视频| 中文字幕精品亚洲无线码一区 | 欧美黑人巨大hd| 香蕉久久夜色| 国产精品自产拍在线观看55亚洲| 99久久99久久久精品蜜桃| 久久久水蜜桃国产精品网| 国产aⅴ精品一区二区三区波| av欧美777| 女人被狂操c到高潮| 自线自在国产av| 欧美乱妇无乱码| 少妇 在线观看| 两个人视频免费观看高清| 精品国产国语对白av| 啦啦啦 在线观看视频| 黑人欧美特级aaaaaa片| 亚洲成人国产一区在线观看| 成人18禁高潮啪啪吃奶动态图| 两个人免费观看高清视频| 夜夜夜夜夜久久久久| 18禁观看日本| av福利片在线| 久久精品人妻少妇| 非洲黑人性xxxx精品又粗又长| 国产一级毛片七仙女欲春2 | 午夜福利在线在线| 精品无人区乱码1区二区| 国产爱豆传媒在线观看 | 亚洲欧洲精品一区二区精品久久久| 中文字幕人成人乱码亚洲影| 精品日产1卡2卡| 村上凉子中文字幕在线| e午夜精品久久久久久久| 97碰自拍视频| 在线看三级毛片| 日韩免费av在线播放| 1024手机看黄色片| 日韩免费av在线播放| 国产av又大| 久久精品91蜜桃| 嫁个100分男人电影在线观看| 黄色成人免费大全| a级毛片在线看网站| 岛国在线观看网站| 99国产精品一区二区蜜桃av| 亚洲最大成人中文| 夜夜爽天天搞| 在线观看66精品国产| 国产精品 国内视频| 亚洲最大成人中文| 99国产综合亚洲精品| 午夜日韩欧美国产| 亚洲狠狠婷婷综合久久图片| 韩国av一区二区三区四区| 国产精品爽爽va在线观看网站 | 欧美黑人欧美精品刺激| 精品久久久久久久久久免费视频| 一二三四社区在线视频社区8| 精品第一国产精品| 一二三四社区在线视频社区8| 久久久久久久久免费视频了| 亚洲第一青青草原| 黄色a级毛片大全视频| 欧美黄色片欧美黄色片| 高潮久久久久久久久久久不卡| 亚洲男人天堂网一区| 制服人妻中文乱码| 欧美黄色片欧美黄色片| 免费在线观看日本一区| 成人免费观看视频高清| 在线观看舔阴道视频| 俄罗斯特黄特色一大片| 亚洲va日本ⅴa欧美va伊人久久| 国产亚洲精品久久久久5区| 无限看片的www在线观看| 黄频高清免费视频| 成人国产综合亚洲| 欧美成人一区二区免费高清观看 | 中文字幕精品免费在线观看视频| 欧美成人性av电影在线观看| 国语自产精品视频在线第100页| 免费搜索国产男女视频| 国产黄片美女视频| 一边摸一边抽搐一进一小说| 欧美成人午夜精品| 最近在线观看免费完整版| 欧美亚洲日本最大视频资源| 免费在线观看黄色视频的| 午夜福利高清视频| 97超级碰碰碰精品色视频在线观看| 亚洲真实伦在线观看| 首页视频小说图片口味搜索| 亚洲国产看品久久| 国产高清videossex| 91麻豆精品激情在线观看国产| 成人18禁高潮啪啪吃奶动态图| 日韩欧美三级三区| 成人欧美大片| 黑丝袜美女国产一区| 黑人操中国人逼视频| 亚洲九九香蕉| 欧美绝顶高潮抽搐喷水| 男女视频在线观看网站免费 | 日韩欧美三级三区| 国产精品电影一区二区三区| 天天添夜夜摸| 在线十欧美十亚洲十日本专区| а√天堂www在线а√下载| 欧美日韩瑟瑟在线播放| 视频在线观看一区二区三区| 欧美亚洲日本最大视频资源| 青草久久国产| 欧美性长视频在线观看| 亚洲av日韩精品久久久久久密| 成人一区二区视频在线观看| 国产精品爽爽va在线观看网站 | 女生性感内裤真人,穿戴方法视频| 真人做人爱边吃奶动态| 嫩草影院精品99| 久久久久免费精品人妻一区二区 | 成人一区二区视频在线观看| 国产精品自产拍在线观看55亚洲| 岛国在线观看网站| 国产成人影院久久av| 精品国产超薄肉色丝袜足j| 久9热在线精品视频| 每晚都被弄得嗷嗷叫到高潮| 国产国语露脸激情在线看| 国产男靠女视频免费网站| 亚洲一区高清亚洲精品| 国内揄拍国产精品人妻在线 | 一级a爱片免费观看的视频| 久久久久国产一级毛片高清牌| 一级a爱片免费观看的视频| 99久久综合精品五月天人人| 久久精品夜夜夜夜夜久久蜜豆 | 精品国产国语对白av| 久久中文看片网| 最近最新中文字幕大全电影3 | 久久久久国产一级毛片高清牌| 国产真人三级小视频在线观看| 成人亚洲精品av一区二区| 草草在线视频免费看| 在线av久久热| 日本五十路高清| 99久久精品国产亚洲精品| av天堂在线播放| 欧美乱妇无乱码| 亚洲专区国产一区二区| 久久九九热精品免费| 99久久久亚洲精品蜜臀av| 亚洲欧美激情综合另类| 一区福利在线观看| 精品国产乱子伦一区二区三区| 亚洲熟妇熟女久久| 看黄色毛片网站| 91av网站免费观看| 成人国产综合亚洲| 亚洲成av人片免费观看| 欧美+亚洲+日韩+国产| 成人国语在线视频| 男人舔女人下体高潮全视频| 最近最新中文字幕大全电影3 | 亚洲男人的天堂狠狠| 一进一出好大好爽视频| 亚洲欧美日韩高清在线视频| 欧美日韩中文字幕国产精品一区二区三区| 亚洲专区字幕在线| 久久精品成人免费网站| 欧美一区二区精品小视频在线| 日本撒尿小便嘘嘘汇集6| 国产成人系列免费观看| 中文在线观看免费www的网站 | 精品熟女少妇八av免费久了| 午夜福利在线在线| 人妻久久中文字幕网| 国产精品免费一区二区三区在线| 热re99久久国产66热| 男女之事视频高清在线观看| 久热爱精品视频在线9| 99在线视频只有这里精品首页| 韩国精品一区二区三区| 特大巨黑吊av在线直播 | 色综合婷婷激情| 两个人视频免费观看高清| 中文亚洲av片在线观看爽| 99久久国产精品久久久| 午夜影院日韩av| 丰满人妻熟妇乱又伦精品不卡| 日韩三级视频一区二区三区| 久久精品国产综合久久久| 夜夜爽天天搞| 丝袜在线中文字幕| 国产亚洲欧美98| 国产高清激情床上av| 精品国产乱子伦一区二区三区| 91成年电影在线观看| 老熟妇乱子伦视频在线观看| 欧美人与性动交α欧美精品济南到| 国产野战对白在线观看| 成人18禁在线播放| 色综合站精品国产| 欧美不卡视频在线免费观看 | 可以在线观看毛片的网站| 免费观看人在逋| 亚洲成人免费电影在线观看| 成熟少妇高潮喷水视频| 国产成年人精品一区二区| 搡老熟女国产l中国老女人| 色播亚洲综合网| 一区二区三区激情视频| 在线看三级毛片| 国产欧美日韩一区二区三| 精品久久久久久久人妻蜜臀av| 亚洲精品中文字幕一二三四区| 午夜两性在线视频| 欧美色欧美亚洲另类二区| 国产蜜桃级精品一区二区三区| 一级毛片女人18水好多| 男女之事视频高清在线观看| 在线观看www视频免费| 午夜老司机福利片| 欧美黑人精品巨大| 热re99久久国产66热| 满18在线观看网站| 999精品在线视频| 国产单亲对白刺激| 国产精品二区激情视频| 欧美 亚洲 国产 日韩一| 一本一本综合久久| 美女免费视频网站| 人妻丰满熟妇av一区二区三区| 一夜夜www| 一区二区三区激情视频| 亚洲成国产人片在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 黄片播放在线免费| 男女之事视频高清在线观看| 亚洲精品中文字幕在线视频| 欧美一级毛片孕妇| 亚洲熟妇熟女久久| 激情在线观看视频在线高清| 黄色片一级片一级黄色片| 午夜a级毛片| 国产激情欧美一区二区| 男女视频在线观看网站免费 | 精品久久久久久久毛片微露脸| 亚洲专区国产一区二区| 97超级碰碰碰精品色视频在线观看| 亚洲激情在线av| 黑人操中国人逼视频| 久久久久九九精品影院| 色在线成人网| 50天的宝宝边吃奶边哭怎么回事| 满18在线观看网站| 日韩免费av在线播放| 国产久久久一区二区三区| 成人亚洲精品一区在线观看| 精品少妇一区二区三区视频日本电影| 亚洲成a人片在线一区二区| 国产精品免费一区二区三区在线| 一本综合久久免费| 欧美色视频一区免费| 午夜福利欧美成人| 亚洲av中文字字幕乱码综合 | 日韩精品青青久久久久久| 天堂影院成人在线观看| 日韩欧美免费精品| 丝袜在线中文字幕| 日韩欧美在线二视频| 日韩视频一区二区在线观看| 精品一区二区三区视频在线观看免费| 少妇的丰满在线观看| 99在线人妻在线中文字幕| 国语自产精品视频在线第100页| 亚洲人成伊人成综合网2020| 亚洲av美国av| 色婷婷久久久亚洲欧美| 欧美日韩福利视频一区二区| 免费av毛片视频| 国产精品亚洲美女久久久| 最近在线观看免费完整版| 国产高清激情床上av| 国产亚洲精品av在线| 午夜精品在线福利| 日本 av在线| 精品国产国语对白av| 两性夫妻黄色片| 一区福利在线观看| 91在线观看av| 欧美色视频一区免费| 日本 av在线| 女人爽到高潮嗷嗷叫在线视频| 男人操女人黄网站| 婷婷六月久久综合丁香| av电影中文网址| 国产熟女午夜一区二区三区| 久热这里只有精品99| 国产精品亚洲美女久久久| 久久香蕉激情| 亚洲国产看品久久| 欧美日本视频| 99在线视频只有这里精品首页| 中文字幕人成人乱码亚洲影| 热re99久久国产66热| 两性夫妻黄色片| 国产精品香港三级国产av潘金莲| 久久国产亚洲av麻豆专区| 99riav亚洲国产免费| 99re在线观看精品视频| 看片在线看免费视频| 亚洲国产欧美一区二区综合| 村上凉子中文字幕在线|