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

利用RNA-Seq發掘玉米葉片形態建成相關的調控基因

2020-02-27 03:50:40郭書磊魯曉民齊建雙魏良明張新韓小花岳潤清王振華鐵雙貴陳彥惠
中國農業科學 2020年1期

郭書磊,魯曉民,齊建雙,魏良明,張新,韓小花,岳潤清,王振華,鐵雙貴,陳彥惠

利用RNA-Seq發掘玉米葉片形態建成相關的調控基因

郭書磊1,2,魯曉民1,齊建雙1,魏良明1,張新1,韓小花1,岳潤清1,王振華1,鐵雙貴1,陳彥惠2

(1河南省農業科學院糧食作物研究所/河南省玉米生物學重點實驗室,鄭州 450002;2河南農業大學農學院,鄭州 450046)

【】葉片寬度和長度等葉形特性是決定植株形態,進而影響種植密度的重要農藝性狀,通過轉錄組測序技術篩選并挖掘玉米葉片形態建成相關的代謝路徑及調控基因,為深入認識葉片發育的分子機理和鑒定葉寬、葉長候選基因奠定基礎。以極端窄葉自交系NL409和寬葉自交系WB665為材料,利用RNA-Seq技術鑒定7葉期第七片葉近基部的差異表達基因(DEGs),通過生物信息學分析,篩選與葉片發育密切相關的代謝通路,利用qRT-PCR驗證不同激素路徑葉形相關基因的表達結果,并結合啟動子區域的序列差異挖掘葉形功能基因。分析對照(WB665)和樣品(NL409)高通量測序結果,在葉寬形成關鍵部位共篩選出5 199個DEGs,其中,2 264(43.55%)個基因表達上調,2 935(56.45%)個基因下調表達,下調基因明顯多于上調基因;GO功能富集分析表明,差異基因主要富集在細胞膜相關的細胞組分中,涉及代謝過程和細胞響應刺激;KEGG富集分析表明,差異基因主要參與到核糖體、植物激素信號轉導、苯丙烷類代謝、乙醛酸和二羧酸代謝等過程,其中核糖體、植物激素信號轉導、鞘脂類代謝下調表達基因較多的路徑與葉片發育密切相關。核糖體路徑富集到多個PRS(PRESSED FLOWER)基因,分析發現()可能在調控窄葉發育過程中發揮重要作用。鞘脂代謝路徑富集的基因幾乎全部下調表達,引起抑制葉片發育的AP1(APETALA1)類和MAPK(Mitogen-Activated Protein Kinase)類基因上調,以及促進葉片發育的()下調,與窄葉發育受抑制的表型一致。植物激素信號轉導路徑富集到的油菜素內酯(BR)響應基因和赤霉素(GA)代謝基因下調,細胞分裂素(CTK)和大部分生長素(Auxin)響應基因上調,與窄葉中DELLA蛋白基因上調表達,抑制GA并促進CTK基因表達的作用模式一致。通過qRT-PCR對18個葉片發育相關基因進行分析,結果表明,其表達趨勢與轉錄組結果一致,分析發現BR相關的、Auxin相關的、e以及TCP類轉錄因子/等基因與窄葉的形成密切相關。明確了一些與玉米葉片發育密切相關的代謝路徑,還發現植物激素間的動態平衡對葉片發育有著重要影響,尤其是生長素與油菜素內酯、細胞分裂素與赤霉素之間的相互作用對調控葉片形態可能發揮重要作用。

玉米;葉寬;葉長;RNA-Seq;形態建成;調控基因

0 引言

【研究意義】玉米作為增產潛力巨大的糧食作物,對于保障糧食安全,促進農業發展具有重要作用。近年來,中國玉米的種植密度維持在4.8—6.77萬株/hm2,與美國8.5—10.95萬株/hm2的種植密度存在顯著差距,提高玉米的耐密植特性,成為實現中國玉米高產突破的重要措施[1-2]。葉片是玉米進行光合作用和抗逆作用的主要營養器官,是耐密植理想株型的關鍵組成因素,對于植株干物質的積累及抗逆性有非常重要的影響。對于大面積密植的農作物,高密度條件下容易出現避蔭綜合癥(shade avoidance syndrome,SAS),迫使群體內個體將能量重新分配,用于葉片和植株伸長生長,從而獲得更多的陽光,盡管避蔭反應對單個植株的適應與生存是有利的,但是由于避蔭反應會削弱植株葉片和儲藏器官的發育[3],對提高大面積種植的作物產量是不利的。高密度條件下,葉片寬大的群體冠層間通風透氣性差,降低植株的抗病性,加劇病蟲害的發生,導致生物量和籽粒產量降低。玉米耐密植株型的合理葉寬、葉長有利于改良理想株型,減少植株的避蔭反應,使個體間競爭最小化,對提高種植密度具有重要的促進作用。葉寬、葉長等葉形特性是決定植株形態進而影響種植密度的重要農藝性狀,因此,選育適宜高密度種植的合理葉寬、葉長玉米,是提高玉米產量的重要技術手段之一。雖然有大量研究揭示了玉米葉發育及形態建成的遺傳基礎,但是鮮有關于葉寬、葉長QTL克隆及調控網絡的研究結果[4],利用轉錄組(RNA-seq)篩選與葉寬、葉長發育相關的調控基因,并發掘其可能的調控路徑,對葉寬、葉長主效基因研究及耐密株型設計育種具有重要意義。【前人研究進展】葉片的發育起始于莖頂端分生組織(shoot apical meristem,SAM)周邊的初始細胞,主要包括葉原基的形成和葉極性軸向的建立[5]。葉原基是由頂端分生組織中心“干細胞”在生長素誘導下形成成簇分布的外緣葉原初生細胞;葉原基形成后在自身遺傳機制和環境因子作用下,建立基-頂軸(葉的基部-頂端)、腹-背軸(莖的近軸面-遠軸面)、中-邊軸(葉的主脈-邊緣)3個極性軸向,這些過程中的發育異常將直接或間接影響葉片的發育和形態建成。玉米葉片的發育和形態建成是由多基因調控的復雜過程,主要由遺傳因素控制,同時受環境影響,葉寬的遺傳相比其他葉形性狀受環境影響較小,而且其具有易觀測的線性形態,更有利于研究葉片發育過程的形態變化和細胞學特征,是研究葉寬最為理想的模式植物[4]。關于擬南芥、水稻、金魚草等植物葉片發育的研究,取得了一系列重要進展,深入解析了葉片發育的調控機理,而玉米葉寬、葉長的研究多為遺傳分析及QTL定位的報道,只有少數文獻關于葉發育基因、、、、、/、、、/和遺傳機理的分析[6-14],揭示了葉寬、葉長發育形成的部分分子機理。高通量轉錄組測序技術,不僅能夠檢測全基因組所有基因的轉錄信息,還可以獲得特定時期某一生物學過程顯著富集的基因,以及特定組織或細胞中的多種轉錄本和優勢表達基因,對研究特定組織的生物學功能具有重要的指導意義,利用轉錄組發掘玉米葉片發育相關基因,為葉片遺傳機理的解析提供了新途徑[15-16]。【本研究切入點】關于葉寬、葉長等葉形主效基因克隆的研究相對有限,已有研究僅揭示了葉形發育的部分分子機理,葉寬、葉長候選基因及其調控路徑有待進一步解析。【擬解決的關鍵問題】本研究通過RNA-seq測序探討葉寬差異極顯著自交系的轉錄表達特性,以期獲得與葉形發育相關的關鍵基因及轉錄調控網絡,更好地理解玉米葉片發育的分子機理,增加對窄葉形成機制的認識,進而為葉形候選基因的克隆及耐密分子育種提供參考。

1 材料與方法

1.1 植物材料

2017年夏極端窄葉自交系NL409(樣品)和寬葉自交系WB665(對照)種植于河南省農業科學院現代農業科技示范基地,在田間正常生長至7葉期時(對生的葉耳包裹著新鮮的幼莖即將打開),取第7片葉近基部葉片(近基部1/3葉長處,切取2.0 cm),為了消除個體間差異,將相同基因型的材料3株混合取樣,每個材料取3次生物學重復,樣品經液氮速凍后-80℃保存,用于RNA提取。極端窄葉自交系NL409是歐洲雜交種連續自交選育而來,寬葉自交系WB665是昌7-2改良系,兩者親緣關系較遠。

1.2 RNA樣品的提取及轉錄組測序

委托北京安諾優達基因科技有限公司完成RNA提取、文庫構建和轉錄組測序。首先,采用酚/氯仿法提取Total RNA,通過NanoDrop 2000微量分光光度計、Agilent 2100 Bioanalyzer和Agilent RNA 6000 Nano Kit檢測RNA濃度、純度和完整性;其次,檢測合格樣品利用Oligo(dT)富集mRNA,并通過Fragmentation Buffer使其短片段化,以短片段mRNA合成雙鏈cDNA,純化后進行末端修復、加堿基A、加測序接頭處理;最后進行PCR擴增完成測序文庫制備。采用Illumina NextSeq 550AR平臺進行雙末端測序,并對獲得的原始序列(raw reads)進行過濾,去除含接頭、低質量(Q<30)、含N比例大于5%的序列,得到高質量的Clean Reads用于后續分析。

1.3 基因表達數據分析

利用TopHat(v2.0.12)軟件,將6個樣本過濾后的Clean Reads比對到B73基因組(AGPv3 版本),通過FPKM(Fragments per Kilobase per Millon Mapped Fragments)對基因表達量進行標準化,以FPKM>1為表達標準,隨后根據極端窄葉與寬葉的比較,用DEseq(|log2FC(fold change)|≥1和-value/FDR<0.05)確定差異表達基因(differentially expressed genes,DEGs)[17]。利用Uniprot、Swissprot、COG、NR、GO和KEGG等數據庫對DEGs進行功能注釋,通過GO富集分析顯示DEGs顯著富集的功能分類[18],用KEGG分析DEGs主要參與的代謝途徑和信號通路(kyoto encyclopedia of genes and genomes,KEGG)。

1.4 實時熒光定量PCR(qRT-PCR)分析

選擇與葉片發育相關的18個差異表達基因進行qRT-PCR驗證。將轉錄組測序用到的RNA,通過反轉錄試劑盒(GoScript Reverse Transcription Kit購于promega生物技術有限公司)合成cDNA第一鏈,以其為模板,利用GoTaq qPCR Master Mix Kit(購于promega生物技術有限公司)在Bio-Rad CFX96實時熒光定量儀上進行qRT-PCR。以玉米中的作為內參基因,每個基因的定量分析設置3次生物學重復,3次技術重復,采用相對定量2-ΔΔCt法分析結果[19]。

2 結果

2.1 葉寬表型和RNA-Seq測序質量分析

極端窄葉自交系NL409和寬葉自交系WB665葉寬形態穩定時期(散粉后一周)的葉寬表型(圖1),以及第7葉期的葉寬達到極顯著差異(表1),在這兩個時期,NL409的葉寬、葉長、葉面積均顯著小于WB665。6個RNA-Seq文庫原始測序過濾后分析結果(表2),各樣本生物學重復間的相關系數為0.984—0.995,其中約88%高質量序列被比對到外顯子區域,表明測序結果良好,數據可用于進一步分析。

表1 自交系NL409和WB665第7葉期以及穩定時期的葉寬、葉長、葉面積

1)7表示第7葉期,S表示穩定時期。不同大寫字母表示在<0.01水平差異顯著

1)7 indicate the 7th leaf period, S indicate the stable period. different uppercase letters indicate significant differences at<0.01 levels

圖1 自交系NL409(左)和WB665(右)穩定時期的葉片

2.2 基因表達數據分析

選取第7葉期窄葉自交系NL409(T7,樣品)和寬葉自交系WB665(C7,對照)葉片發育最寬處進行轉錄組差異比較分析,發現5 199個差異表達基因,其中2 264(43.55%)個基因上調表達,2 935(56.45%)個基因下調表達。上(下)調1—3倍(fold change)的基因分別占差異表達上(下)調基因總數的29.37%(665)和21.69%(756);上(下)調3—10倍的基因分別占43.86%(993)和29.24%(1 019);上(下)調10—60倍的基因分別占20.76%(470)和12.82%(447);上(下)調60倍以上的基因分別占6.01%(136)和24.29%(713)。從火山圖中可以直觀地看出對照組與比較樣本組之間的表達水平具有顯著差異,而且差異60倍以上的下調基因明顯多于上調基因(圖2)。

表2 樣品reads分布情況

1)71、72、73分別代表7葉期的第1、2、3個樣品

1)71, 72 and 73 indicate the first, second and third samples of the 7th leaf period respectively

GO功能分析表明,差異表達基因被注釋到24個生物過程、22個細胞組分、19個分子功能類別中,差異基因主要富集在細胞膜相關的細胞組分中。其中生物過程顯著富集在次生代謝、花粉雌蕊相互作用、花粉識別、防衛反應和細胞識別5個過程(圖3);細胞組分分析中差異基因顯著富集在細胞外緣、質膜、細胞膜固有成分、細胞膜基本成分、膜部分、共質體、胞間接合和胞間連絲等10個組分類別(圖4);分子功能顯著富集在單加氧酶活性、氧化還原活性、催化活性、纖維素合成酶活性、ADP結合、碳水化合物結合、鐵離子結合、葡糖基轉移酶活性等19個類別(圖5)。

橫坐標為差異表達基因表達倍數變化,縱坐標為表達量變化的統計學顯著程度

利用MapMan工具對差異表達基因進行pathway富集分析,差異基因被注釋到125條代謝通路中,主要涉及生物合成、能量代謝、信號轉導及次生代謝。以-value≤0.05為篩選標準,發現這些基因顯著富集到11條代謝路徑(表3),其中核糖體路徑富集差異基因數目最多達到83個,37個基因表達上調,46個基因下調表達;其次,植物激素信號轉導通路富集57個差異基因,24個基因上調,33個基因下調;苯丙烷類代謝路徑18個基因上調,15個基因下調;乙醛酸和二羧酸代謝通路15個基因上調,10個基因下調;谷胱甘肽代謝路徑13個基因上調,10個基因下調;鞘脂類代謝路徑2個基因上調,11個基因下調(表4);油菜素內酯生物合成路徑4個基因上調,3個基因下調(表5);類黃酮生物合成路徑6個基因上調,2個基因下調;苯并惡唑嗪類生物合成和苯乙烯類、二芳庚酸類、姜酚類生物合成通路也富集到較多差異基因,分別有7個和5個基因下調表達。基于不同代謝通路中基因上下調比例,推斷核糖體、植物激素信號轉導、鞘脂類代謝這些基因下調表達較多的路徑可能與葉片發育,尤其是窄葉的形成密切相關。

Size代表圓點的大小,表示富集到該GO條目基因的數量,q-value表示該GO條目的富集程度,顏色越趨近于紅色表示富集程度越高。下同

圖4 差異表達基因顯著富集的細胞組分

表3 差異表達基因顯著富集的代謝通路

植物內源激素廣泛參與到葉片發育過程,對葉片形態建成起到一定調節作用。進一步分析植物激素信號轉導路徑發現,在葉寬形成的關鍵部位檢測到大量參與激素信號響應和轉導的差異基因,涉及生長素、油菜素內酯、細胞分裂素、赤霉素、脫落酸、乙烯、茉莉酸等多種激素(圖6)。分析不同響應路徑(表5和圖6)發現,生長素信號轉導通路富集到18個差異基因,有12個響應基因表達上調;油菜素內酯信號響應路徑富集到6個差異基因,全部下調表達;細胞分裂素響應路徑富集的3個基因,表達均上調;赤霉素和水楊酸響應通路分別富集到1個差異基因,表達均上調;脫落酸路徑富集到8個基因下調表達,4個上調表達基因;乙烯路徑中4個基因上調表達,4個基因下調表達;茉莉酸路徑中1個基因上調表達,7個基因下調表達。

表4 鞘脂類代謝路徑及其相關基因

FC:表達量差異倍數(NL409/WB665)

FC:Differential multiple of expression (NL409/WB665)

圖5 差異表達基因顯著富集的分子功能

2.3 熒光定量PCR驗證RNA-Seq結果

研究表明多種轉錄因子參與葉細胞分化及葉片形態建成,如YABBY、MADS類轉錄因子決定葉原細胞命運向不同方向分化形成不同組織,GRF類調控特定部位細胞的增殖和生長;ARF、SBP、TCP、Myb、HB類調控葉細胞的分化和維管的形成;同時ARF、ARR、GRAS類在葉片發育過程分別參與生長素、細胞分裂素、赤霉素信號的響應和生物合成,影響葉片形態建成。為了進一步驗證RNA-Seq結果的準確性,選取與葉片發育密切相關的18個不同的基因進行了qRT-PCR定量分析。對qRT-PCR與RNA-Seq結果進行相關性分析,發現兩者的相關系數2=0.962,說明RNA-Seq結果良好。

3 討論

3.1 葉片發育密切相關的代謝通路

多條路徑顯著富集到大量與玉米葉片發育相關的基因。擬南芥()和()分別編碼不同亞族的核糖體蛋白(PRS),與rRNA形成蛋白復合物促進葉片近軸組織的建立[20-21]。本文核糖體路徑顯著富集到多個和(電子附表1),并鑒定出玉米()啟動子區域存在顯著差異(電子附圖1),其在窄葉中下調表達的作用模式(圖7)與擬南芥、缺失突變導致葉柄變長、葉片變窄的作用機理相似。鞘脂類(Sphingolipid)是細胞膜骨架及脂蛋白的重要組成成分,鞘脂代謝在細胞生長、不同功能細胞分化、細胞內以及細胞間信號傳遞發揮重要作用,鞘脂合成紊亂含量較低時,能夠激活MAPK(mitogen-activated protein kinase)路徑和轉錄因子AP1(APETALA1),進而調控細胞增殖、分化[22-23]。抑制擬南芥花萼葉和莖葉形成,()能夠解除的抑制作用,促進葉片發育[24]。本文鞘脂代謝路徑富集到的基因大部分下調表達(表4),類和MAPK類基因表達上調(表4),表達下調(表4),與窄葉發育受抑制的表型一致。研究人員分析葉片發育的動態轉錄組學,不僅發現大量激素信號基因,細胞壁合成前體、纖維素合成酶、蛋白酶體、翻譯后修飾基因也被大量檢測到[25]。Bolduc等[9]解析調控網絡,不僅在葉片和分生組織中發現大量激素合成和響應基因,而且富集到脂類代謝、次生代謝、細胞壁合成和調控RNA轉錄相關的基因。葉片發育過程不同部位、不同時期的蛋白質組學分析,也富集到大量的尿苷二磷酸-葡糖基/葡聚糖轉移酶、細胞壁前體、細胞壁受體激酶、纖維素合成酶和RNA轉錄相關的基因[26],本文顯著富集的路徑與前人的研究結果一致。由此可知,上述顯著富集的代謝通路(表3)與葉片發育密切相關,這些代謝路徑中的基因可能參與葉寬、葉長的形態建成。

Tryptophan metabolism:色氨酸代謝;Auxin:生長素;Ubiquitin mediated proteolysis:泛素介導的蛋白水解;Cell enlargement:細胞增大;Plant growth:植物生長;Zeatin biosynthesis:玉米素的生物合成;Cytokinine:細胞分裂素;Cell division:細胞分裂;Shoot initiation:根的起始;Diterpenoid biosynthesis:二萜類化合物的生物合成;Gibberellin:赤霉素;Stem growth:莖的生長;Induced germination:誘導萌發;Carotenoid biosynthesis:類胡蘿卜素的生物合成;Abscisic acid:脫落酸;Stomatal closure:氣孔閉合;Seed dormancy:種子休眠;Cysteine and methionine metabolism:半胱氨酸和蛋氨酸代謝;Ethylene:乙烯;Endoplasmic reticulum(ER):內質網;Fruit ripening:果實成熟;Senescence:衰老;Brassinosteroid biosynthesis:油菜素內酯生物合成;Brassinosteroid:油菜素甾醇;Proteasomal degradation:蛋白酶體降解;Cell elongation:細胞伸長;Alpha-Linolenic acid metabolism:α-亞麻酸代謝;Jasmonic acid:茉莉酸;Monoterpenoid biosynthesis類單萜生物合成;Indole alkaloid biosynthesis:吲哚生物堿生物合成;Stress response:應激反應;Phenylalanine metabolism:苯丙氨酸代謝;Salicylic acid:水楊酸;Disease resistance:抗病性。紅色表示富集到的差異基因表達上調,綠色表示基因下調Red indicates that DEGs are up-regulated, and green indicates down-regulation of DEGs

表5 植物激素信號轉導路徑相關基因

FC:表達量差異倍數(NL409/WB665);BR:油菜素內酯Brassinosteroid;CTK:細胞分裂素Cytokinine;GAs:赤霉素信號Gibberellin signaling;GAmrg:赤霉素代謝相關基因Gibberellin metabolism related genes

圖7 RNA-Seq與qRT-PCR結果比較

本研究發現油菜素內酯信號轉導路徑顯著富集的基因全部下調(表5),不僅檢測到前人已研究表達量降低的BZR類和BRI類基因[9],還發現BR路徑中BSK類、TCH類基因也下調表達。生長素信號轉導路徑富集到的基因大部分上調表達(表5),其中ARF(GRMZM2G078274、GRMZM2G030710和GRMZM5G874163)、AUX-IAA(GRMZM2G104176、GRMZM2G077356、GRMZM2G079957、GRMZM2G074427)和Auxin transporter(GRMZM2G127949)在調控葉片發育過程中發揮著重要作用[9]。此外,高表達時,抑制GA信號及葉基部細胞的生長[25],本文DELLA家族蛋白(GRMZM2G013016)在近基部窄葉中高表達(表5,圖6),且赤霉素代謝相關的基因幾乎全部下調表達(表5),這與窄葉發育受抑制的結果一致。由此推斷油菜素內酯、生長素、赤霉素在葉片發育過程發揮重要作用。

3.2 不同激素途徑發掘的葉形重要功能基因

植物器官發育過程中極性建立是器官形態建成的核心。葉片的發育起始于葉原基細胞的分裂增殖,伴隨著一系列特定基因的精確調控,葉原基建立多維空間的極性軸向,促使葉原基細胞朝特定的方向分化增殖,并最終影響葉片的形態和大小。

TCP類轉錄因子通過AUX、BR、GA、CTK等多種激素調控葉片、花發育以及植株形態建成等生長過程。葉片發育過程可調控影響擬南芥葉片及植株形態,而BRI是發揮功能所必需的;此外,多個負向調控葉片邊緣區細胞分化增殖,其中、分別負調控、和、影響生長素分布及含量[9,27-28,31]。過量表達TCP II型基因()使西紅柿和擬南芥的葉片變小,或表現卷曲葉[27-28],擬南芥的同源基因/(AC205574.3_FG006)在窄葉中上調表達的作用模式(圖7,電子附表1)與前人研究結果一致。因此,推測/可能影響多種激素調控玉米葉片發育。

激素BR可影響玉米葉片形態建成。編碼一種細胞色素P450蛋白,作用于葉片縱向極性細胞的伸長生長,其功能缺失影響擬南芥葉柄變短,葉片變大變寬[29],擬南芥BR合成相關的P450家族和突變后影響葉細胞變短變少[30-31]。(VIVIPAROUS1)、(VP1/ABI3-LIKE)、(Related to ABI3/VP1)、(Related to ABI3/ VP1-Like)分別編碼不同的B3結構DNA結合蛋白,擬南芥缺失突變后影響葉片變小且赤霉素含量降低[32],作為BR受體負向調控擬南芥蓮座葉的生長[33],與本文中下調和表達上調(電子附表2),以及赤霉素路徑基因下調表達的結果一致(表5)。通過E-box(CACCTG/ CACGTG/CACATG)結合位點激活BR響應和合成基因的表達,維持BR平衡對生長的調節作用[34-35],本文中BR響應基因的啟動子區域存在大量該元件(電子附圖4—9),說明植物激素BR與葉片發育密切相關。此外,(GRMZM2G143235)和(GRMZM2G344521)在葉寬形成的關鍵部位上調表達(圖7),且位于已定位的葉寬區間內[36],同時啟動子區域的序列差異(電子附圖2),表明和可能在葉長和葉寬發育過程中發揮著重要作用,由此推測BR通過某種途徑影響玉米葉長、葉寬的形成。

生長素合成、分布異常影響葉片形態。玉米與水稻()和()同源,、、和編碼一種生長素合成酶,調控葉原基的形成,其功能缺失突變影響葉細胞的分裂和維管束的形成,導致水稻葉片變窄[37-40],與本文在窄葉中下調表達(圖7)的作用模式一致。編碼一種AGO7型蛋白,是合成ta-siARF所必需的,通過ta-siARF調控的轉錄,影響玉米葉片腹背面生長素分布,導致突變體葉片變窄卷曲[41]。此外,位于葉寬研究的熱點區域[4],而且其啟動子區域序列存在顯著差異(電子附圖3)。由此可知,(GRMZM2G480386)、e(GRMZM5G892991)可能通過影響生長素調控窄葉的形成。

玉米葉片發育過程大量與生長素相關的基因參與其形成。玉米的ChIP分析富集大量生長素AUX-IAA和ARF基因,轉錄組分析發現生長素基因(、)在玉米葉細胞朝不同功能區生長的過渡區上調表達[9],與本研究中、(圖7)表達結果一致。擬南芥受葉原基近軸端特異表達tasiR-ARF的抑制,在近軸端不表達,在遠軸端()與、相互作用的同時,影響、和表達上調,促進葉原細胞的分化及遠軸組織的形成[42-43];miR165/166抑制HD-ZIP lll類基因、在葉原基遠軸端的表達,促進近軸組織的形成[44]。、和(圖7)位于已定位的葉寬位點內[4,36]。此外,、、(圖7)分別受miR167、miR164、miR166調控并且與葉長、葉寬位點顯著相關[36],由此可知(GRMZM2G017709)、(GRMZM5G874163)、(GRMZM2G175827)、(GRMZM2G078274)、(GRMZM2G055585)、(GRMZM2G029692)在葉片發育過程可能發揮特定作用。

植物激素間的的相互作用可調節葉片發育。BR信號路徑中的BIN2激酶能夠磷酸化ARF2,導致其結合Auxin響應元件的功能喪失,影響擬南芥葉片等器官增大[45-46],同時BR能夠誘導BZR1靶基因的表達[47]。與相互作用促進擬南芥下胚軸變長[48]。Auxin能夠誘導擬南芥表達,同時通過BR路徑促進細胞分裂生長[49]。此外,Auxin還可促進表達,加強BR信號轉導,而擬南芥通過影響生長素促進細胞膨大生長[27],ARF能夠特異結合Auxin響應元件(TGTCTC)[50],本文BR響應基因啟動子富含該元件(電子附圖4—電子附圖9),而且BR響應基因全部下調,Auxin響應基因上調(表5,圖7),這暗示本文中(GRMZM2G030710)、(GRMZM2G078274)(GRMZM2G065635)(GRMZM2G15877)(表5,電子附表1)可能介導Auxin和BR的平衡來調節植物器官以及葉片的發育。

HB類轉錄因子包含多個KNOX、bHLH和HD-ZIP lll基因。()和屬于KNOX家族,在頂端分生組織特定表達,調控細胞分化向膨大生長的過渡,影響葉片多維空間的發育,尤其在基部-末端軸向的形成過程中,與、、相互作用調控葉片發育[9]。KNOX蛋白能夠促進細胞分裂素(CTK)合成基因上調,同時抑制GA3和GA20氧化酶合成,促進SAM細胞分裂形成葉原基[9,51-52];提高GA含量可以解除DELLA對葉片發育的抑制作用[52]。與本文中KNOX家族基因(和)(圖7)、CTK路徑基因(表5)、DELLA基因上調(表5),GA氧化酶基因下調(表5)影響葉片變窄的結果一致。GA可通過DELLA蛋白負向調控CTK響應基因進而影響CTK變化[53]。擬南芥GA信號的抑制因子SPY(SPINDLY)能夠正向調控CTK促進植株及葉片的發育[54]。由此可知(GRMZM2G017709)和(GRMZM2G433591)可能通過影響細胞分裂素與赤霉素的平衡參與玉米窄葉的形成。

生長調節因子GRF的互作因子AN3/GIF1作為葉片發育過程的調節開關,主要影響葉細胞橫向分裂,分析其調控機理,在葉片基部分裂區富集到大量上調表達的GRF基因,在葉部伸長區僅發現和高量表達,其他GRF基因在葉部伸長區表達量逐漸降低,過量表達和使玉米葉片變小[13-14],與本文中檢測到(GRMZM2G041223)和(GRMZM2G018414)在窄葉葉片發育的過渡區域上調表達(圖7)的趨勢一致,同時和位于葉寬研究的熱點區域[4]。(GRMZM2G143235)編碼一種TGA結合域的堿性亮氨酸拉鏈蛋白,主要在SAM和葉片與葉鞘連接處表達,作用于葉基部葉耳、葉舌的形成,影響葉片的發育[55]。ZmARF2(GRMZM2G082836)是ADP核糖基化轉錄因子,能夠調控擬南芥葉面積和株高[56]。由此推斷、、、(圖7)在玉米葉片發育以及葉寬形成過程中可能發揮重要作用。

4 結論

明確了一些與玉米葉片發育密切相關的代謝路徑,發現植物激素間的動態平衡對葉片發育有著重要影響,尤其是Auxin與BR、CKT與GA之間的相互作用對調控葉片形態可能發揮重要作用。

[1] 明博, 謝瑞芝, 侯鵬, 李璐璐, 王克如, 李少昆. 2005—2016年中國玉米種植密度變化分析. 中國農業科學, 2017, 50(11): 1960-1972.

MING B, XIE R Z, HOU P, LI L L, WANG K R, LI S K. Changes of maize planting density in China., 2017, 50(11): 1960-1972. (in Chinese)

[2] 郭書磊, 陳娜娜, 齊建雙, 岳潤清, 韓小花, 燕樹鋒, 盧彩霞, 傅曉雷, 郭新海, 鐵雙貴. 不同密度下玉米倒伏相關性狀與產量的研究. 玉米科學, 2018, 26(5): 71-77.

GUO S L, CHEN N N, QI J S, YUE R Q, HAN X H, YAN S F, LU C X, FU X L, GUO X H, TIE S G. Study on the relationship between yield and lodging traits of maize under different planting densities., 2018, 26(5): 71-77. (in Chinese)

[3] YANG C W, XIE F M, JIANG Y P, LI Z, HUANG X, LI L. Phytochrome A negatively regulates the shade avoidance response by increasing auxin/indole acidic acid protein stability., 2018, 44(1): 29-41.

[4] 郭書磊, 張君, 齊建雙,岳潤清, 韓小花, 燕樹鋒, 盧彩霞, 傅曉雷, 陳娜娜, 庫麗霞, 鐵雙貴. 玉米葉形相關性狀的Meta-QTL及候選基因分析. 植物學報, 2018, 53(4): 487-501.

GUO S L, ZHANG J, QI J S,YUE R Q, HAN X H, YAN S F, LU C X, FU X L, CHEN N N, KU L X, TIE S G. Analysis of meta-quantitative trait loci and their candidate genes related to leaf shape in maize., 2018, 53(4): 487-501. (in Chinese)

[5] 崔曉峰, 黃海. 葉發育的遺傳調控機理研究進展. 植物生理學報, 2011, 47(7): 631-640.

CUI X F, HUANG H. Recent progresses from studies of leaf development., 2011, 47(7): 631-640. (in Chinese)

[6] MORENO M A, HARPER L C, KRUEGER R W, DELLAPORTA S L, FREELING M. liguleless1 encodes a nuclear-localized protein required for induction of ligules and auricles during maize leaf organogenesis., 1997, 11(5): 616-628.

[7] WALSH J, WATERS C A, FREELING M. The maize geneliguleless2 encodes a basic leucine zipper protein involved in the establishment of the leaf blade-sheath boundary., 1998, 12(2): 208-218.

[8] ZHANG X L, MADI S, BORSUK L, NETTLETON D, ELSHIRE R J, BUCKNER B, BUCKNER D J, BECK J, TIMMERMANS M, SCHNABLE P S, SCANLON M J. Laser microdissection of narrow sheath mutant maize uncovers novel gene expression in the shoot apical meristem., 2007, 3(6): 1040-1052.

[9] BOLDUC N, YILMAZ A, MEJIA-GUERRA M K, MOROHASHI K, O'CONNOR D, GROTEWOLD E, HAKE S. Unraveling the KNOTTED1 regulatory network in maize meristems., 2012, 26(15): 1685-1690.

[10] TIMMERMANS M C, SCHULTES N P, JANKOVSKY J P, NELSON T. Leafbladeless1 is required for dorsoventrality of lateral organs in maize., 1998, 125(15): 2813-2823.

[11] DOUGLAS R N, WILEY D, SARKAR A, SPRINGER N, TIMMERMANS MARJA C P, SCANLON M J. Ragged seedling2 Encodes an ARGONAUTE7-like protein required for mediolateral expansion, but not dorsiventrality, of maize leaves., 2010, 22(5): 1441-1451.

[12] Candela H, Johnston R, Gerhold A, FOSTER T, HAKE S. The milkweed pod1 gene encodes a KANADI protein that is required for abaxial/adaxial patterning in maize leaves., 2008, 20(8): 2073-2087.

[13] NELISSEN H, EECKHOUT D, DEMUYNCK K, PERSIAU G, WALTON A, BEL M V, VERVOORT M, CANDAELE J, BLOCK J D, AESAERT S, LIJSEBETTENS M V, GOORMACHTIG S, VANDEPOELE K, LEENE J V, MUSZYNSKI M, GEVAERT K, INZE D, JAEGER G D. Dynamic changes in ANGUSTIFOLIA3 complex composition reveal a growth regulatory mechanism in the maize leaf., 2015, 27(6): 1605-1619.

[14] WU L, ZHANG D, XUE M, QIAN J J, HE Y, WANG S C. Overexpression of the maize GRF10, an endogenous truncated growth-regulating factor protein, leads to reduction in leaf size and plant height., 2014, 56(11): 1053-1063.

[15] 吳慶飛, 秦磊, 董雷, 丁澤紅, 李平華, 杜柏娟. 玉米光合突變體hcf136(high chlorophyll fluorescence 136)的轉錄組分析. 作物學報, 2018, 44(4): 493-504.

WU Q F, QIN L, DONG L, DING Z H, LI P H, DU B J. Transcriptome analysis on a maize photosynthetic mutant hcf136 (high chlorophyll fluorescence 136)., 2018, 44(4): 493-504. (in Chinese)

[16] 邱化榮, 周茜茜, 何曉文, 張宗營, 張世忠, 陳學森, 吳樹敬. 基于轉錄組分析蘋果水楊酸特異響應基因MdWRKY40的啟動子鑒定. 中國農業科學, 2017, 50(20): 3970-3990.

QIU H R, ZHOU Q Q, HE X W, ZHANG Z Y, ZHANG S Z, CHEN X S, WU S J. Identification of MdWRKY40 promoter specific response to salicylic acid by transcriptome sequencing., 2017, 50(20): 3970-3990. (in Chinese)

[17] TRAPNELL C, WILLAMS B A, PERTEA G, MORTAZAVI A, KWAN G, BAREN M J, SALZBERG S L, WOLD B JHTER L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation., 2010, 28(5): 511-515.

[18] WANG L K, FENG Z X, WANG X, WANG X W, ZHANG X G. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data., 2009, 26(1): 136-138.

[19] LIVAK K J, SCHMITTGEN T D. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCTmethod., 2001, 25(4): 402-408.

[20] SCHIPPERS J H M, MUELLER-ROEBER B. Ribosomal composition and control of leaf development., 2010, 179(4): 307-315.

[21] YAO Y, LING Q H, WANG H, HUANG H. Ribosomal proteins promote leaf adaxial identity., 2008, 135(7): 1325-1334.

[22] SPIEGEL S, MERRILL Jr A H. Sphingolipid metabolism and cell growth regulation., 1996, 10(12): 1388-1397.

[23] MERRILL Jr A H, SANDHOFF K. Sphingolipids: Metabolism and cell signaling//.Amsterdam, Elsevier Science Press, 2002, 373-407.

[24] MONNIAUX M, MCKIM S M, CARTOLANO M, THEVENON E, PARCY F, MILTANTIS T, HAY A. Conservation vs divergence in LEAFY and APETALA1 functions betweenand., 2017, 216(2): 549-561.

[25] LI P H, PONNALA L, GANDOTRA N, WANG L, SI Y Q, TAUSTA S L, KEBROM T H, PROVART N, PATEL R, MYERS C R, REIDEL E J, TURGEON R, LIU P, SUN Q, NELSON T, Brutnell T P. The developmental dynamics of the maize leaf transcriptome., 2010, 42(12): 1060-1067.

[26] FACETTE M R, SHEN Z X, Bj?rnsdóttir F R, BRIGGS S P, SMITH L G. Parallel proteomic and phosphoproteomic analyses of successive stages of maize leaf development., 2013: 2798-2812.

[27] NICOLAS M, CUBAS P. The role of TCP transcription factors in shaping flower structure, leaf morphology, and plant architecture//:,. Cambridge, MA: Academic Press, 2016: 249-267.

[28] MANASSERO N G U, VIOLA I L, WELCHEN E, GONZALEZ D H. TCP transcription factors: architectures of plant form., 2013, 4(2): 111-127.

[29] KIM G T, TSUKAYA H, UCHIMIYA H. The ROTUNDIFOLIA3 gene ofencodes a new member of the cytochrome P-450 family that is required for the regulated polar elongation of leaf cells., 1998, 12(15): 2381-2391.

[30] FUJIOKA S, LI J, CHOI Y H, SETO H, TAKATSUTO S, NOGUCHI T, WATANABE T, Kuriyama H, YOKOT T, CHORY J, SAKURAI A. Thedeetiolated2 mutant is blocked early in brassinosteroid biosynthesis., 1997, 9(11): 1951-1962.

[31] GUO Z, FUJIOKA S, BLANCAFLOR E B, MIAO S, GOU X P, LI J. TCP1 modulates brassinosteroid biosynthesis by regulating the expression of the key biosynthetic gene DWARF4 in., 2010, 22(4): 1161-1173.

[32] SUZUK M, WANG H H Y, MCCARTY D R. Repression of the LEAFY COTYLEDON 1/B3 regulatory network in plant embryo development by VP1/ABSCISIC ACID INSENSITIVE 3-LIKE B3 genes., 2007, 143(2): 902-911.

[33] HU Y X, WANG Y H, LIU X F, LI J Y.RAV1 is down-regulated by brassinosteroid and may act as a negative regulator during plant development., 2004, 14(1): 8.

[34] YUAN W Y, LUO X, LI Z C, YANG W N, WANG Y Z, LIU R, DU J M, HE Y H. A cis cold memory element and a trans epigenome reader mediate Polycomb silencing of FLC by vernalization in., 2016, 48(12): 1527-1534.

[35] JE B I, PIAO H L, PARK S J,PARK S H, KIM C M, XUAN Y H, PARK S H, HUANG J, CHOI Y D, AN G, WONG H L, FUJIOKA S, KIM M C, SHIMAMOTO K, HAN C. RAV-Like1 maintains brassinosteroid homeostasis via the coordinated activation of BRI1 and biosynthetic genes in rice., 2010, 22(6): 1777-1791.

[36] TIAN F, BRADBURY P J, BROWN P J, HUNG H Y, SUN Q, FLINT-GARCIN S, ROCHEFORD T R, MCMULLEN M D, HOLLAND J B, BUCKLER E S. Genome-wide association study of leaf architecture in the maize nested association mapping population., 2011, 43(2): 159-165.

[37] Ishiwata A, Ozawa M, Nagasaki H, KATO M, NODA Y, YAMAGUCHI T, NOSAKA M, SHIMIZU-SATO S,NAGASAKI A,MAEKAWA M,HIRANO H Y,SATO Y. Two WUSCHEL-related homeobox genes, narrow leaf2 and narrow leaf3, control leaf width in rice., 2013, 54(5): 779-792.

[38] KUBO F C, YASUI Y, KUMAMARU T, SATO Y, HIRANO H Y. Genetic analysis of rice mutants responsible for narrow leaf phenotype and reduced vein number., 2016, 91(4): 235-240.

[39] JIANG D, FANG J J, LOU L, ZHAO J F, YUAN S J, YIN L, SUN W, PENG L X, GUO B T, LI X Y. Characterization of a null allelic mutant of the rice NAL1 gene reveals its role in regulating cell division., 2015, 10(2): e0118169.

[40] FUJINO K, MATSUDA Y, OZAWA K, NISHIMURA T, KOSHIBA T, FRAAIJE M W, SEKIGUCHI H. NARROW LEAF 7 controls leaf shape mediated by auxin in rice., 2008, 279(5): 499-507.

[41] DOUGLAS R N, WILEY D, SARKAR A, SPRINGER N, TIMMERMANSM C P, SCANLON M J.ragged seedling2 encodes an ARGONAUTE7-like protein required for mediolateral expansion, but not dorsiventrality, of maize leaves., 2010: 1441-1451.

[42] DOTTO M C, PETSCH K A, AUKERMAN M J, BEATTY M, HAMMELL M, TIMMERMANS M C P. Genome-wide analysis of leafbladeless1-regulated and phased small RNAs underscores the importance of the TAS3 ta-siRNA pathway to maize development., 2014, 10(12): e1004826.

[43] XIE Y, STRAUB D, EGUEN T, BRANDT R, STAHL M, Jaime F. MARTINEZ-GAECIA J F, WENKEL S. Meta-analysis ofKANADI1 direct target genes identifies a basic growth-promoting module acting upstream of hormonal signaling pathways., 2015, 169(2): 1240-1253.

[44] EMERY J F, FLOYD S K, AIVAREZ J, ESHED Y, HAWKER N P, IZHAKI A, BAUM S F, BOWMAN J L. Radial patterning ofshoots by class III HD-ZIP and KANADI genes., 2003, 13(20): 1768-1774.

[45] OKUSHIMA Y, MITINA I, QUACH H L, THEOLOGIS A. AUXIN RESPONSE FACTOR 2 (ARF2): a pleiotropic developmental regulator., 2005, 43(1): 29-46.

[46] SCHRUFF M C, SPIELMAN M, TIWARI S, ADAMS S, FENBY N, SCOTT R J. The AUXIN RESPONSE FACTOR 2 gene oflinks auxin signalling, cell division, and the size of seeds and other organs., 2006, 133(2): 251-261.

[47] TIAN H Y, LV B S, DING T T, BAI M Y, DING Z J. Auxin-BR interaction regulates plant growth and development., 2018, 8: 2256.

[48] OH E, ZHU J Y, BAI M Y, ARENHART R A, SUN Y, WANG Z Y. Cell elongation is regulated through a central circuit of interacting transcription factors in thehypocotyl., 2014, 3: e03031.

[49] CHUNG Y, MAHARJAN P M, LEE O, FUJIOKA S, JANG S, KIM B, TAKATSUTO S, TSUJIMOTO M, KIM H, CHO S, PARK T, CHO H, HWANG I, CHOE S. Auxin stimulates DWARF4 expression and brassinosteroid biosynthesis in., 2011, 66(4): 564-578.

[50] SAKAMOTO T, MORINAKA Y, INUKAI Y, KITANO H, FUJIOKA S. Auxin signal transcription factor regulates expression of the brassinosteroid receptor gene in rice., 2013, 73(4): 676-688.

[51] JASINSKI S, PIAZZA P, CRAFT J, HAY A, WOOLLEY L, RIEU I, PHILLIPS A, HEDDEN P, TSIANTIS M. KNOX action inis mediated by coordinate regulation of cytokinin and gibberellin activities., 2005, 15(17): 1560-1565.

[52] NELISSEN H, RYMEN B, JIKUMARU Y, DEMUYNCK K, LIJSEBETTENS M V, KAMIYA Y, INZE D, BEEMSTER G T S. A local maximum in gibberellin levels regulates maize leaf growth by spatial control of cell division., 2012, 22(13): 1183-1187.

[53] FONOUNI-FARDE C, KISIALA A, BRAULT M, EMERY R J N, ANOUCK D, FLORIAN F. DELLA1-mediated gibberellin signaling regulates cytokinin-dependent symbiotic nodulation., 2017, 175(4): 1795-1806.

[54] GREENBOIM-WAINBERG Y, MAYMON I, BOROCHOV R, ALVAREZ J, OLSZEWSKI N, ORI N, ESHED Y, WEISS D. Cross talk between gibberellin and cytokinin: theGA response inhibitor SPINDLY plays a positive role in cytokinin signaling., 2005, 17(1): 92-102.

[55] WALSH J, WATERS C A, FREELING M. The maize geneliguleless2 encodes a basic leucine zipper protein involved in the establishment of the leaf blade-sheath boundary., 1998, 12(2): 208-218.

[56] WANG Q L, XUE X J, LI Y L, DONG Y B, ZHANG L, ZHOU Q, DENG F, MA Z Y, QIAO D H, HU C H, REN Y L. A maize ADP-ribosylation factor ZmArf2 increases organ and seed size by promoting cell expansion in., 2016, 156(1): 97-107.

Explore Regulatory Genes Related to Maize Leaf Morphogenesis Using RNA-Seq

GUO ShuLei1,2, LU XiaoMin1, QI JianShuang1, WEI LiangMing1, ZHANG Xin1, HAN XiaoHua1, YUE RunQing1, WANG ZhenHua1, TIE ShuangGui1, CHEN YanHui2

(1Cereal Crops Institute, Henan Academy of Agricultural Sciences/Henan Provincial Key Lab of Maize Biology, Zhengzhou 450002;2College of Agronomy, Henan Agricultural University, Zhengzhou 450046)

【】Leaf shape characteristics are one of important agronomic traits that determine plant morphology and affect planting density. However, the molecular mechanism related to leaf shape remain unknown in maize. Here, transcriptome sequencing technology was used to screen and explore genome-wide analysis of regulatory genes and metabolic pathways involved in leaf morphogenesis. This study will lay the foundation for further understanding the regulator mechanism of leaf development in plant and identifying candidate genes of leaf shapes, such as leaf width and leaf length.【】Extreme narrow-leaf inbred line NL409 and wide-leaf line WB665 were selected as the experimental materials. By RNA-Seq technology, the differentially expressed genes (DEGs) of the seventh leaf base between these two lines were identified during the 7th leaf stage. Furthermore, metabolic pathways closely related to leaf development were also analyzed using a series of bioinformatics analysis. qRT-PCR was used to validate the expression level of DEGs in different hormone pathways, and the further promoter analysis were performed to explore leaf-shape functional genes.【】By analyzing the high-throughput sequencing in WB665 and NL409, a total of 5 199 DEGs were obtained at the primary section of leaf width formation. Of which, 2 264 (43.55%) genes were up-regulated, whereas down-regulated genes were significantly more than up-regulated genes with 2 935 (56.45%) decreased genes. GO enrichment analysis showed that these DEGs were mainly enriched in cell membrane-associated function terms of cellular components, including metabolic process and cell stimulus response. KEGG enrichment analysis showed that these DEGs were mainly involved in ribosome, plant hormone signal transduction, sphingolipid metabolism pathways, phenylpropanoid biosynthesis, glyoxylate and dicarboxylate metabolism and other processes. among which ribosome, plant hormone signal transduction, sphingolipid metabolism pathways with more down-regulated genes were closely related to leaf development. One of PRS (PRESSED FLOWER) family genes, which were enriched in the ribosomal pathway in this study,PRS13 (PFL2) was identified to participate in regulating the development of narrow leaves. The expression pattern of genes enriched in sphingolipid metabolism pathway and its related MAP kinase, AP1-like, and LFY-like were consistent with the result of the inhibited development of narrow leaves. Notably, all of BR (Brassinosteroid) response genes and most of GA (Gibberellin) metabolic genes were down-regulated in plant hormone signal transduction pathway, while the expression level of all the CTK (Cytokinine) response genes and Auxin genes are mostly increased. The action of up-regulated expression of DELLA protein gene affecting the GA and CTK pathways was consistent with the phenotypic result of narrow leaves. Eighteen genes were validated by qRT-PCR. The result showed that the expression trend was consistent with the transcriptome data. Moreover, the BR-related, auxin-related,and TCP-like transcription factorswere identified to be closely associated with the formation of narrow leaves.【】Summarily, this study unveils several metabolic pathways closely related to leaf development in maize, and find the dynamic balance between plant hormones plays an important role in leaf development, especially the interaction between Auxin and BR as well as CTK and GA.

maize (); leaf width; leaf length; RNA-Seq; morphogenesis; regulatory gene

10.3864/j.issn.0578-1752.2020.01.001

2019-05-20;

2019-07-11

中國博士后科學基金(2017M612404)、河南省科技攻關計劃項目(182102110122)、河南省農業科學院優秀青年科技基金

郭書磊,E-mail:guosl1309@163.com。通信作者王振華,E-mail:wzh201@126.com。通信作者鐵雙貴,E-mail:tieshuangg@126.com。通信作者陳彥惠,E-mail:chy9890@163.com

(責任編輯 李莉)

主站蜘蛛池模板: 亚洲妓女综合网995久久| 久久国产精品嫖妓| 亚洲综合第一区| 91尤物国产尤物福利在线| 亚洲视频免| 亚洲一区二区三区在线视频| 欧美日韩专区| 99精品在线看| 亚洲男人的天堂在线| 永久免费av网站可以直接看的 | 蜜桃视频一区二区三区| 日韩国产综合精选| 国产精品成人观看视频国产| 亚洲无线一二三四区男男| 国产精品第页| 欧美中文字幕第一页线路一| 亚洲精品无码日韩国产不卡| 精品欧美日韩国产日漫一区不卡| 99ri精品视频在线观看播放| 国产成人毛片| 999精品免费视频| 99福利视频导航| 少妇高潮惨叫久久久久久| 久久99蜜桃精品久久久久小说| 五月综合色婷婷| 青青青国产视频手机| 亚洲国产精品一区二区第一页免| 国产成人8x视频一区二区| 成人免费一级片| 久久天天躁夜夜躁狠狠| 国产在线专区| 四虎国产精品永久一区| 一级毛片在线免费视频| 亚洲精选高清无码| 国产精品永久不卡免费视频| 日韩精品毛片人妻AV不卡| 亚欧乱色视频网站大全| 内射人妻无套中出无码| 日韩午夜福利在线观看| 992tv国产人成在线观看| 久久久噜噜噜| 波多野结衣爽到高潮漏水大喷| 日本午夜精品一本在线观看 | 亚洲美女久久| 蜜桃视频一区二区| 亚洲国产成人精品一二区| 精品国产女同疯狂摩擦2| 欧美成人影院亚洲综合图| 国产99在线| 92午夜福利影院一区二区三区| 99免费视频观看| 最新亚洲人成网站在线观看| 国产精品lululu在线观看 | 天天综合色天天综合网| 午夜视频免费试看| 国产伦精品一区二区三区视频优播| 久久久久久久97| 免费高清a毛片| 国产精品一区不卡| 免费一级成人毛片| 精品国产三级在线观看| 日本午夜三级| 91无码视频在线观看| 日韩精品一区二区三区中文无码| 国产迷奸在线看| 亚洲成人黄色在线| 无码免费的亚洲视频| 特级做a爰片毛片免费69| 91无码视频在线观看| 99热免费在线| 久久亚洲黄色视频| 亚洲国产综合精品一区| 国外欧美一区另类中文字幕| 重口调教一区二区视频| 免费黄色国产视频| 二级特黄绝大片免费视频大片| 国内精品小视频在线| 国产麻豆aⅴ精品无码| 色哟哟国产精品一区二区| 成人午夜亚洲影视在线观看| 超碰精品无码一区二区| 日韩在线1|