寧夏枸杞花藥轉(zhuǎn)錄組分析及花藥發(fā)育相關(guān)基因的篩選

范文強(qiáng)，劉雪霞，馬小蘭，王路瑤，張 興，岳思君，唐建寧，朱金忠，姚新靈，鄭 蕊*

寧夏枸杞花藥轉(zhuǎn)錄組分析及花藥發(fā)育相關(guān)基因的篩選

范文強(qiáng)1，劉雪霞1，馬小蘭1，王路瑤1，張 興1，岳思君1，唐建寧2，朱金忠3，姚新靈1，鄭 蕊1*

1. 寧夏大學(xué)生命科學(xué)學(xué)院，西部特色生物資源保護(hù)與利用教育部重點(diǎn)實(shí)驗(yàn)室，寧夏優(yōu)勢(shì)特色作物現(xiàn)代分子育種重點(diǎn)實(shí)驗(yàn)室，寧夏 銀川 750021 2. 寧夏枸杞產(chǎn)業(yè)發(fā)展中心，寧夏 銀川 750021 3. 中寧縣杞鑫枸杞苗木專業(yè)合作社，寧夏 中寧 755100

為了探知寧夏枸杞花藥發(fā)育過程及挖掘花藥發(fā)育相關(guān)基因。以單雙核花粉時(shí)期和成熟花粉時(shí)期的花藥為研究材料，通過高通量Illumina HiSeq2500測(cè)序平臺(tái)進(jìn)行轉(zhuǎn)錄組分析。組裝共獲得128 151條Unigene，將Unigene與NCBI的非冗余蛋白質(zhì)數(shù)據(jù)庫(kù)（non-redundant protein databasse，NR）、去冗余的蛋白序列數(shù)據(jù)庫(kù)（Swiss prot protein database，Swiss-Prot）、基因本體數(shù)據(jù)庫(kù)（gene ontology，GO）、直源同源群集數(shù)據(jù)庫(kù)（clusters of orthologous groups，COG）、真核同源群數(shù)據(jù)庫(kù)（euKaryotic orthologous groups，KOG）、基因功能和代謝途徑數(shù)據(jù)庫(kù)（Kyoto encyclopedia of genes and genomes，KEGG）數(shù)據(jù)庫(kù)比對(duì)，共有40 140條得到注釋。在錯(cuò)誤發(fā)現(xiàn)率（false discovery rate，F(xiàn)DR）＜0.05且差異倍數(shù)（fold change，F(xiàn)C）≥2的條件下，單雙核花粉時(shí)期與成熟花粉時(shí)期的差異表達(dá)基因共有46個(gè)。KEGG注釋發(fā)現(xiàn)有4個(gè)差異表達(dá)基因參與甘油酯代謝、亞油酸代謝、苯丙氨酸代謝和苯丙烷類生物合成4個(gè)途徑中。富含亮氨酸重復(fù)蛋白、細(xì)胞周期蛋白、類受體蛋白激酶、E3泛素蛋白連接酶等11個(gè)差異表達(dá)基因參與了花藥發(fā)育過程。隨機(jī)選取其中5個(gè)差異表達(dá)基因進(jìn)行qRT-PCR驗(yàn)證，表達(dá)結(jié)果與轉(zhuǎn)錄組數(shù)據(jù)基本一致，表明其數(shù)據(jù)真實(shí)可靠。為進(jìn)一步解釋寧夏枸杞花藥發(fā)育的相關(guān)分子機(jī)制奠定了基礎(chǔ)，也為后續(xù)研究寧夏枸杞花藥發(fā)育的生物學(xué)特征提供了有益的基因資源。

寧夏枸杞；花藥發(fā)育；轉(zhuǎn)錄組；差異表達(dá)基因；生物學(xué)特征

寧夏枸杞L.屬于茄科（Solanaceae）枸杞屬Linn．落葉灌木，是重要的藥食兼用植物，具有很高的食用、藥用和生態(tài)保護(hù)價(jià)值[1]。枸杞是我國(guó)傳統(tǒng)的名貴中藥材和重要經(jīng)濟(jì)作物，具有調(diào)節(jié)免疫、潤(rùn)肺明目等功效，枸杞葉、果和根中都含有人體所需的蛋白質(zhì)、氨基酸和微量元素[2]。枸杞子是一味常用的補(bǔ)肝益腎中藥，現(xiàn)代醫(yī)學(xué)研究證實(shí)其含有胡蘿卜素、維生素A、維生素C、維生素B1及鈣、磷、鐵、鋅、錳等營(yíng)養(yǎng)成分，對(duì)造血功能有促進(jìn)作用，還具有抗衰老、抗腫瘤、抗脂肪肝及降血糖等作用[3]。

花主要由花萼、花瓣、雄蕊和雌蕊4部分組成，其中花藥作為雄蕊最重要的組成部分，其發(fā)育過程復(fù)雜且精細(xì)，首先由雄蕊原基經(jīng)過減數(shù)分裂形成小孢子四分體結(jié)構(gòu)，然后花藥繼續(xù)膨大生長(zhǎng)，小孢子由四分體結(jié)構(gòu)發(fā)育成花粉粒，藥室破裂釋放花粉，最終花藥組織衰老退化[4]。研究植物的花藥發(fā)育，對(duì)于提高植物的抗逆能力、培育新品種以及提高產(chǎn)量具有重要的理論意義與實(shí)踐價(jià)值[5]?；ㄋ幇l(fā)育相關(guān)基因的研究在許多植物已有報(bào)道，在番茄中，通過RNAi，使基因編碼的9-順式-環(huán)氧類胡蘿卜素雙加氧酶（9--epoxycarotenoid dioxygenase，NCED）活性降低，導(dǎo)致番茄花藥發(fā)育相關(guān)基因、、和的轉(zhuǎn)錄下調(diào)，導(dǎo)致花藥發(fā)育異常，從而產(chǎn)生了大量異?；ǚ踇6]；百合花藥發(fā)育相關(guān)的基因可以在突變體中適當(dāng)表達(dá)，產(chǎn)生正常的花粉粒[7]；水稻中基因的嚴(yán)格表達(dá)對(duì)水稻花藥發(fā)育至關(guān)重要[8]。關(guān)于枸杞生長(zhǎng)發(fā)育[9]，生理生態(tài)[10]以及基因克隆等方面的研究已有相關(guān)報(bào)導(dǎo)，其中枸杞基因已經(jīng)被克隆并進(jìn)行功能驗(yàn)證，結(jié)果表明該基因參與枸杞花藥發(fā)育過程中淀粉供能的調(diào)控[2]。枸杞花藥的研究尚處于起步階段[11]。由于枸杞基因組和轉(zhuǎn)錄組相關(guān)信息相對(duì)缺乏，成為制約研究枸杞花藥發(fā)育及其調(diào)控機(jī)制的一個(gè)因素。因此，通過轉(zhuǎn)錄組學(xué)技術(shù)挖掘并研究一些與枸杞花藥發(fā)育相關(guān)基因的表達(dá)調(diào)控，對(duì)于解析枸杞花藥發(fā)育的分子機(jī)制具有重要意義。

目前，基于高通量測(cè)序技術(shù)的轉(zhuǎn)錄組分析，已經(jīng)在植物、動(dòng)物、微生物轉(zhuǎn)錄組研究中廣泛應(yīng)用[12]。因此，本實(shí)驗(yàn)利用RNA-Seq技術(shù)，選擇以體外更加容易培養(yǎng)的單雙核花粉時(shí)期和減數(shù)分裂完成以后的成熟花粉時(shí)期枸杞花藥進(jìn)行轉(zhuǎn)錄組測(cè)序，得到的Unigene進(jìn)行NCBI的非冗余蛋白質(zhì)數(shù)據(jù)庫(kù)（non-redundant protein databasse，NR）、基因本體數(shù)據(jù)庫(kù)（gene ontology，GO）、直源同源群集數(shù)據(jù)庫(kù)（clusters of orthologous groups，COG）、真核同源群數(shù)據(jù)庫(kù)（euKaryotic orthologous groups，KOG）、基因功能和代謝途徑數(shù)據(jù)庫(kù)（Kyoto encyclopedia of genes and genomes，KEGG）、去冗余的蛋白序列數(shù)據(jù)庫(kù)（Swiss prot protein database，Swiss-Prot）、蛋白質(zhì)家族數(shù)據(jù)庫(kù)（protein family，Pfam）注釋，篩選出與花藥發(fā)育相關(guān)的差異表達(dá)基因，結(jié)果為進(jìn)一步研究枸杞花藥發(fā)育提供基因資源和研究基礎(chǔ)。

1 材料

“寧杞１號(hào)”是寧夏枸杞研究所鐘鉎元[13]選育的高產(chǎn)、優(yōu)質(zhì)、適應(yīng)性強(qiáng)的枸杞新品種。于2016年6月取自寧夏育新枸杞種業(yè)有限公司（106°6'11.516"E，38°31'1.808"N）標(biāo)準(zhǔn)化管理的枸杞實(shí)驗(yàn)田經(jīng)寧夏大學(xué)生命科學(xué)學(xué)院鄭蕊博士鑒定為寧夏枸杞L.。參照徐青等[14]的方法，采集枸杞單雙核花粉時(shí)期（S1）和成熟花粉時(shí)期（S2）的花蕾，設(shè)置3個(gè)生物學(xué)重復(fù)，S1時(shí)期和S2時(shí)期的3個(gè)樣品分別編號(hào)為A1～A3和B1～B3?；ɡ俚牟杉Y(jié)合形態(tài)觀察、縱橫徑測(cè)量和花藥壓片法，于干冰上剝離不同時(shí)期花蕾中的花藥（圖1），并立即置于液氮中速凍，長(zhǎng)期保存于?80 ℃冰箱備用。

2 方法

2.1 RNA提取及測(cè)序文庫(kù)構(gòu)建

RNA提取參照TaKaRa公司的RNAisoplus產(chǎn)品說明書進(jìn)行，提取RNA樣品濃度由微量紫外檢測(cè)儀NanoDrop測(cè)定，并使用瓊脂糖凝膠電泳檢測(cè)RNA完整性，以確保枸杞轉(zhuǎn)錄組測(cè)序文庫(kù)滿足后續(xù)分析的需要。

A-S1時(shí)期的花蕾和花藥 B-S2時(shí)期的花蕾和花藥

2.2 測(cè)序及原始數(shù)據(jù)的處理

通過高通量illuminaHiSeqTM2500測(cè)序平臺(tái)對(duì)已構(gòu)建好的cDNA文庫(kù)進(jìn)行測(cè)序，并獲得下機(jī)數(shù)據(jù)即原始讀段（raw reads）。利用Trinity軟件對(duì)測(cè)序獲得的Clean data進(jìn)行組裝，獲得枸杞轉(zhuǎn)錄組Unigene數(shù)據(jù)庫(kù)，通過測(cè)序錯(cuò)誤率和GC含量對(duì)測(cè)序數(shù)據(jù)的質(zhì)量進(jìn)行評(píng)估，得到Clean data進(jìn)行分析。

2.3 差異基因的功能注釋

使用BLAST軟件將Unigene與NR、GO、COG、KOG、KEGG、Swiss-Prot、Pfam等數(shù)據(jù)庫(kù)比對(duì)（BLAST參數(shù)-value≤1.0×10?5和HMMER 參數(shù)-value≤1.0×10?10），獲得Unigene注釋信息。

2.4 花藥發(fā)育相關(guān)差異表達(dá)基因的篩選

使用DESeq軟件對(duì)“寧杞１號(hào)”花藥S1和S2 2個(gè)發(fā)育時(shí)期的差異表達(dá)基因進(jìn)行分析。在差異表達(dá)分析過程中采用公認(rèn)有效的Benjamini-Hochberg方法對(duì)原有假設(shè)檢驗(yàn)得到的顯著性值進(jìn)行校正，并最終采用校正后的值，即錯(cuò)誤發(fā)現(xiàn)率（false discovery rate，F(xiàn)DR）＜0.05為差異表達(dá)基因篩選的關(guān)鍵指標(biāo)，以降低對(duì)大量基因的表達(dá)值進(jìn)行獨(dú)立的統(tǒng)計(jì)假設(shè)檢驗(yàn)帶來(lái)的假陽(yáng)性。在篩選過程中，將FDR＜0.05且差異倍數(shù)（fold change，F(xiàn)C）≥2作為篩選標(biāo)準(zhǔn)。其中，F(xiàn)C表示2樣品（組）間表達(dá)量的比值。

2.5 qRT-PCR分析

本研究為了驗(yàn)證轉(zhuǎn)錄組數(shù)據(jù)的可靠性，以枸杞組成型表達(dá)基因?yàn)閮?nèi)參（GenBank登陸號(hào)：HQ415754.1），隨機(jī)選取5個(gè)相對(duì)表達(dá)量較高的差異表達(dá)基因進(jìn)行qRT-PCR驗(yàn)證，設(shè)計(jì)qRT-PCR特異引物（表1），利用反轉(zhuǎn)錄試劑盒將RNA反轉(zhuǎn)錄得到cDNA，采用TB Green Premix Ex Taq II試劑盒，在熒光定量PCR儀上進(jìn)行定量檢測(cè)，樣品設(shè)置3個(gè)生物學(xué)重復(fù)，采用2?ΔΔCt法計(jì)算差異基因的相對(duì)表達(dá)量。

表1 qRT-PCR引物

Table 1 Primers used for real-time quantitative PCR

基因ID基因名稱正向引物反向引物 c68267.graph_c0PRX5’-CACCTGTGGATAGGAAGATCA-3’5’-GGAGCACTGTGGATTCATACT-3’ c68844.graph_c0AGPs5’-AGTGAGTGCGAACAATTAAC-3’5’-GTGGTAGTAAGAGATGCGTT-3’ c77007.graph_c0MBR15’-ACGGAAGCTCTGCACATAG-3’5’-GATGACATGGATATTCGTGAT-3’ c75034.graph_c0RPK5’-GAGTATATCTGATGGATATA-3’5’-ACAGGATTAGGTGGCCTA-3’ c72608.graph_c0LRR5’-GATGAACATTGAGTTAACTA-3’5’-CTTCACAACCATTTCAAT-3’ HQ415754.1LbActin5’-GACCTTCAATGTTCCCGCTATG-3’5’-GCCATCACCAGAGTCCAACAC-3’

2.6 統(tǒng)計(jì)學(xué)方法

本研究使用統(tǒng)計(jì)學(xué)軟件SPSS22和Excel2016對(duì)數(shù)據(jù)進(jìn)行分析。

3 結(jié)果與分析

3.1 枸杞花藥轉(zhuǎn)錄組測(cè)序和clean data組裝

樣品經(jīng)過轉(zhuǎn)錄組測(cè)序以后，共獲得44.72 Gb clean data，Q30堿基百分比在92.07%以上（表2），表明轉(zhuǎn)錄組數(shù)據(jù)豐富，質(zhì)量較高。利用Trinity軟件對(duì)clean data數(shù)據(jù)進(jìn)行組裝，共得到128151條Unigene，序列總長(zhǎng)度為93 244 314 bp，平均長(zhǎng)度為727.61 bp，N50為1170 bp，其中長(zhǎng)度區(qū)間位于200～300 bp的Unigene數(shù)量最多，為42 785條，占33.39%（表3），表明組裝完整性較高，可用于后續(xù)分析。

3.2 Unigene功能注釋

使用BLAST軟件將Unigene序列與NR、GO、COG、KOG、KEGG、Swiss-Prot、Pfam等7個(gè)數(shù)據(jù)庫(kù)比對(duì)，最終獲得40 140個(gè)有注釋信息的Unigene，占Unigene總數(shù)的31.32%（表4）。

3.3 差異表達(dá)基因分析

采用DESeq方法篩選，設(shè)置篩選差異條件FDR＜0.05，且FC≥2，作為篩選標(biāo)準(zhǔn)，共獲得差異表達(dá)基因46個(gè)，且相對(duì)于S1時(shí)期，這些基因均在S2時(shí)期下調(diào)表達(dá)（圖1）。

表2 轉(zhuǎn)錄組測(cè)序產(chǎn)量

Table 2 Statistics of transcriptome sequencing output

樣品讀取數(shù)堿基數(shù)GC/%Q30/% A129 290 3758 634 418 32643.1792.44 A223 477 5066 921 070 68243.0292.87 A321 597 6166 404 336 10043.5192.14 B122 641 0256 720 401 74042.9992.07 B228 977 3308 608 669 92443.1092.46 B325 044 7647 431 459 82442.8692.36

表3 組裝結(jié)果

Table 3 Statistics of assembly results

序列長(zhǎng)度/nt基因（占比/%） 200～30042 785（33.39） 300～50034 095（26.61） 500～100027 280（21.29） 1000～200014 326（11.18） ≥20009665（7.54） 總數(shù)/條128 151 總長(zhǎng)度/bp93 244 314 N50長(zhǎng)度/bp1170 平均長(zhǎng)度/bp727.61

表4 Unigene注釋信息

Table 4 Unigene annotation information

數(shù)據(jù)庫(kù)Unigene數(shù)量百分率/% NR39 35230.71 GO17 50113.66 COG9 184 7.17 KOG19 04014.86 KEGG11 092 8.66 Pfam22 73917.74 Swiss-prot21 53016.80 總計(jì)40 14031.32

3.3.1 差異表達(dá)基因的GO功能分類 對(duì)差異表達(dá)基因進(jìn)行GO分析，結(jié)果表明有8個(gè)差異表達(dá)基因分別被注釋到生物學(xué)過程（biological process，BP）細(xì)胞組分（cellular component，CC）和分子功能（molecular function，MF）中（圖2）。

其中，分子功能中催化活性（catalytic activity）和蛋白結(jié)合（binding）所占比例最高，分別為7個(gè)（87.5%）和6個(gè)（75%）；細(xì)胞組分中最多的是細(xì)胞膜部分（membrane），共有1個(gè)，占12.5%；生物學(xué)過程中則是代謝過程（metabolic process）、細(xì)胞過程（cellular process）和單組織過程（single-organism process）比重較大，分別為4個(gè)（50%）、3個(gè)（37.5%）和2個(gè)（25%），如表5所示。

3.3.2 差異表達(dá)基因COG/KOG分類 對(duì)差異表達(dá)基因進(jìn)行COG分類（圖3），結(jié)果顯示，有17個(gè)差異表達(dá)基因分布于8類基因家族，注釋最多的是僅一般功能預(yù)測(cè)（general function prediction only，4）；其次是翻譯（translation，3），復(fù)制、重組和修復(fù)（reproduction、reorganization and repair，3），信號(hào)轉(zhuǎn)導(dǎo)機(jī)制（signal transduction mechanisms，3）；注釋最少的是氨基酸的運(yùn)輸和代謝（transport and metabolism of amino acids，1），脂質(zhì)的運(yùn)輸和代謝（transport and metabolism of lipids，1），碳水化合物運(yùn)輸與代謝（carbohydrate transport and metabolism，1）以及次生產(chǎn)物合成運(yùn)輸及代謝（secondary metabolites biosynthesis transport and catabolism，1）。

綠色和紅色的點(diǎn)代表有顯著性表達(dá)差異的基因，綠色代表基因表達(dá)量下調(diào)，紅色代表基因表達(dá)量上調(diào)，黑色的點(diǎn)代表無(wú)顯著性表達(dá)差異的基因

KOG分類發(fā)現(xiàn)共有10個(gè)基因被注釋到5個(gè)功能中（圖4），其中，含轉(zhuǎn)錄本數(shù)目最多的是僅一般功能的預(yù)測(cè)（general function prediction only，5）；其次是翻譯后修飾（post-translational modifications，3）、蛋白質(zhì)轉(zhuǎn)換（protein turnover，3）和伴侶蛋白（chaperones，3)；最少的是細(xì)胞周期調(diào)控（cell cycle control，1）、細(xì)胞分裂（celldivision，1）、染色體分裂（chromosome partitioning，1）、氨基酸的轉(zhuǎn)運(yùn)和代謝（amino acid transport and metabolism，1）以及信號(hào)轉(zhuǎn)導(dǎo)機(jī)制（signal transduction mechanisms，1）。

圖2 差異表達(dá)基因GO二級(jí)節(jié)點(diǎn)注釋

表5 差異表達(dá)基因GO功能分類

Table 5 GO functional classification of differentially expressed genes

GO功能分類功能注釋DEGs數(shù)量百分率/% CC細(xì)胞膜112.5 MF催化活性蛋白結(jié)合7687.575.0 BP代謝過程細(xì)胞過程單組織過程43250.037.525.0

3.3.3 差異表達(dá)基因KEGG注釋 對(duì)差異表達(dá)基因進(jìn)行KEGG注釋，發(fā)現(xiàn)有3個(gè)差異表達(dá)基因分別參與甘油酯代謝、亞油酸代謝、苯丙氨酸代謝和苯丙烷類生物合成4條代謝通路（表6）。

橫坐標(biāo)為COG各分類內(nèi)容，縱坐標(biāo)為基因數(shù)目

橫坐標(biāo)為KOG各分類內(nèi)容，縱坐標(biāo)為基因數(shù)目

表6 差異表達(dá)基因KEGG功能注釋

Table 6 KEGG functional annotation ofdifferentially expressed genes

編號(hào)KEGG功能注釋基因數(shù)量基因ID ko00561甘油酯代謝1c77096.graph_c1 ko00591亞油酸代謝1c31079.graph_c0 ko00940苯丙烷類生物合成1c68267.graph_c0 ko00360苯丙氨酸代謝1c68267.graph_c0

3.4 花藥發(fā)育相關(guān)基因篩選

通過對(duì)枸杞轉(zhuǎn)錄組測(cè)序結(jié)果進(jìn)行功能注釋、功能分類及代謝途徑分析發(fā)現(xiàn)，共篩選到11個(gè)與花藥發(fā)育相關(guān)的差異表達(dá)基因（表7）。

3.5 差異表達(dá)基因的qRT-PCR驗(yàn)證

隨機(jī)選取5個(gè)差異表達(dá)基因進(jìn)行qRT-PCR分析。5個(gè)差異表達(dá)基因即過氧化物酶（peroxidase，PRX）、阿拉伯半乳糖蛋白（arabinogalactan protein，AGP）、E3泛素蛋白連接酶（E3 ubiquitin-protein ligase MBR1）、類受體蛋白激酶（receptor protein kinase-like，RPK）和富含亮氨酸重復(fù)蛋白（leucine-rich repeat protein，LRR），在S2時(shí)期的表達(dá)量相對(duì)于S1時(shí)期的表達(dá)量均下調(diào)，且5個(gè)差異表達(dá)基因表達(dá)水平在S2時(shí)期極顯著低于S1時(shí)期，與轉(zhuǎn)錄組測(cè)序結(jié)果基本一致（圖5），表明轉(zhuǎn)錄組數(shù)據(jù)結(jié)果可靠。

表7 花藥發(fā)育相關(guān)差異表達(dá)基因

Table 7 Differentially expressed genes related to anther development

基因ID基因描述變化倍數(shù)基因表達(dá) c68267.graph_c0過氧化物酶?4.590 996 613下調(diào) c77096.graph_c1二?；视王；D(zhuǎn)移酶?2.919 934 769下調(diào) c75034.graph_c0過氧化物酶?3.479 978 950下調(diào) c72608.graph_c0富含亮氨酸重復(fù)蛋白?3.232 330 297下調(diào) c73778.graph_c0硫氧還蛋白?2.932 284 773下調(diào) c77007.graph_c0E3泛素蛋白連接酶?3.571 161 469下調(diào) c68844.graph_c0阿拉伯半乳糖蛋白?3.730 754 896下調(diào) c77579.graph_c0限制性內(nèi)切酶?3.097 051 914下調(diào) c77438.graph_c0內(nèi)切酶?2.956 411 029下調(diào) c69446.graph_c0細(xì)胞周期蛋白?3.158 118 976下調(diào) c31079.graph_c0脂氧合酶?3.124 745 110下調(diào)

***代表差異極顯著，P＜0.001

3.6 相關(guān)性分析

利用統(tǒng)計(jì)學(xué)軟件對(duì)S1和S2 2個(gè)時(shí)期的轉(zhuǎn)錄組數(shù)據(jù)和qRT-PCR結(jié)果進(jìn)行相關(guān)性分析，結(jié)果表明S1、S2 2個(gè)時(shí)期差異表達(dá)基因的qRT-PCR結(jié)果與轉(zhuǎn)錄組數(shù)據(jù)相關(guān)系數(shù)分別為0.25、0.41，呈正相關(guān)，說明轉(zhuǎn)錄組數(shù)據(jù)可靠。

4 討論

本研究對(duì)“寧杞１號(hào)”花藥發(fā)育過程中S1時(shí)期和S2時(shí)期的花藥進(jìn)行轉(zhuǎn)錄組測(cè)序，測(cè)序結(jié)果與NR、GO、COG、KOG、KEGG、Swiss-Prot、Pfam 7個(gè)數(shù)據(jù)庫(kù)比對(duì)，有注釋信息的Unigene占Unigene總數(shù)的31.32%，未被注釋到的Unigene占68.68%，可能原因是這些Unigene為非編碼序列，無(wú)法進(jìn)行同源序列比對(duì)；另外，由于枸杞基因組信息缺乏，這些沒有被注釋到的基因可能是一些未知的基因。對(duì)于篩選的差異基因數(shù)目較少，推測(cè)原因可能是大多數(shù)基因在S1時(shí)期和S2時(shí)期差異表達(dá)不明顯。

花藥發(fā)育轉(zhuǎn)錄組分析顯示，在S1和S2 2個(gè)時(shí)期，有46個(gè)基因差異表達(dá)，且在S2時(shí)期全部下調(diào)表達(dá)。46個(gè)差異表達(dá)基因中，有33個(gè)基因在NR數(shù)據(jù)庫(kù)中得到注釋。將33個(gè)基因在PubMud數(shù)據(jù)庫(kù)檢索發(fā)現(xiàn)，PRX、二酰基甘油?；D(zhuǎn)移酶、RPK、富含亮氨酸重復(fù)基因、硫氧還蛋白、E3泛素蛋白連接酶、AGP、限制性內(nèi)切酶等11個(gè)基因（表7）與擬南芥[15-16]等植物花藥發(fā)育過程有關(guān)。GO和KEGG功能分析表明，有4個(gè)差異表達(dá)基因被注釋到甘油酯代謝、亞油酸代謝、苯丙氨酸代謝和苯丙烷類生物合成4個(gè)代謝途徑中，其中編碼PRX的基因既在苯丙氨酸代謝中表達(dá)，也在苯丙烷生物合成中表達(dá)，說明該基因的功能與這2個(gè)代謝途徑比較密切。

PRX是苯丙烷類生物合成中的一類關(guān)鍵酶，苯丙烷類生物合成作為一個(gè)重要的代謝通路，是苯丙素代謝網(wǎng)絡(luò)的核心。苯丙素生物合成會(huì)影響植物直立生長(zhǎng)，植物體內(nèi)長(zhǎng)途水分運(yùn)輸[16]以及花器官的發(fā)育[17-18]。徐青等[14]的研究結(jié)果顯示，S1時(shí)期，絨氈層細(xì)胞退化，降解產(chǎn)物流入藥室，匯聚形成中央大液泡；S2時(shí)期，絨氈層細(xì)胞完全消失，表明絨氈層降解產(chǎn)物是花藥發(fā)育必須的營(yíng)養(yǎng)物質(zhì)；Jacobowitz等[15]的研究發(fā)現(xiàn)，基因是擬南芥花藥發(fā)育過程中維持絨氈層和小孢子細(xì)胞壁完整性必不可少的，PRX基因突變會(huì)導(dǎo)致絨氈層的肥大而導(dǎo)致雄性不育；同時(shí)，Hecht等[16]的研究發(fā)現(xiàn)，受體樣蛋白激酶2基因?qū)儆贚RR-RLK家族，是花藥中層和絨氈層分化所必需的，受體樣蛋白激酶2基因突變的花藥中不存在中間層，且絨氈層顯示肥大，表明受體樣蛋白激酶2會(huì)影響絨氈層的發(fā)育來(lái)控制花粉的成熟。本研究結(jié)果顯示，基因和受體樣蛋白激酶2基因在花藥發(fā)育S2時(shí)期的表達(dá)量下調(diào)，由此推測(cè)，在S1時(shí)期花藥發(fā)育需要絨氈層的降解產(chǎn)物提供營(yíng)養(yǎng)，此時(shí)基因和受體樣蛋白激酶2基因的表達(dá)量較高，而S2時(shí)期花粉已完全成熟，不再需要絨氈層的降解來(lái)提供營(yíng)養(yǎng)，所以這2個(gè)基因的表達(dá)量較S1時(shí)期下調(diào)。這與前人的研究結(jié)果一致，因此推測(cè)和受體樣蛋白激酶2基因是S1時(shí)期降解絨氈層為花藥發(fā)育提供營(yíng)養(yǎng)物質(zhì)必需的基因。此外，甘油酯代謝過程受到二?；视王；D(zhuǎn)移酶基因（diacylglycerol acyltransferase，DGAT）的表達(dá)調(diào)控，Banilas等[19]的研究結(jié)果顯示，基因在成熟或衰老的橄欖組織中高度表達(dá)，同時(shí)基因轉(zhuǎn)錄物也在花藥和子房發(fā)育的后期積累，表明基因表達(dá)主要影響著花藥成熟時(shí)期的發(fā)育；這與本課題組的研究結(jié)果不同，推測(cè)可能原因是S2時(shí)期花藥的發(fā)育受到基因和其他基因共同調(diào)控。

隨機(jī)選取5個(gè)差異表達(dá)基因進(jìn)行qRT-PCR驗(yàn)證和相關(guān)性分析，qRT-PCR結(jié)果與轉(zhuǎn)錄組數(shù)據(jù)一致。相關(guān)性分析結(jié)果表明，S1、S2 2個(gè)時(shí)期5個(gè)差異表達(dá)基因的qRT-PCR結(jié)果與轉(zhuǎn)錄組數(shù)據(jù)呈正相關(guān)，說明轉(zhuǎn)錄組數(shù)據(jù)真實(shí)可靠。本研究篩選到的差異表達(dá)基因，在花藥發(fā)育過程中發(fā)揮重要的作用，這些差異表達(dá)基因的發(fā)現(xiàn)為后續(xù)研究枸杞花藥發(fā)育提供基因資源和奠定研究基礎(chǔ)。

利益沖突 所有作者均聲明不存在利益沖突

[1] Zhao J H, Li H, Yin Y,. Fruit ripening inandis associated with distinct gene expression patterns [J].2020, 10(8): 1550-1567.

[2] 張興, 周麗, 喬亞蕊, 等. 枸杞Lb14-3-3c基因克隆及轉(zhuǎn)化馬鈴薯的研究 [J]. 植物遺傳資源學(xué)報(bào), 2019, 20(6): 1523-1534.

[3] 陳廷卓. 枸杞子的化學(xué)成分與藥理作用研究進(jìn)展 [J]. 中醫(yī)臨床研究, 2010, 2(22): 24.

[4] 胡若琳, 袁超, 牛義, 等. 植物MYB轉(zhuǎn)錄因子在花藥發(fā)育中的調(diào)控作用 [J]. 生物工程學(xué)報(bào), 2020, 36(11): 2277-2286.

[5] 張?jiān)趯? 趙海, 胡夢(mèng)輝, 等. 組學(xué)在花藥發(fā)育研究中的應(yīng)用進(jìn)展I: 轉(zhuǎn)錄組學(xué) [J]. 生物技術(shù)進(jìn)展, 2019, 9(5): 433-439.

[6] Wang X H, Xu N, Dong K,. SlNCED1 affects pollen maturation in tomato by regulating the expression of anther-specific genes [J]., 2021, 95(2): 191-205.

[7] Yuan G, Wu Z, Liu X,. Characterization and functional analysis of LoUDT1, a bHLH transcription factor related to anther development in the lily oriental hybrid Siberia (spp.) [J]., 2021, 166: 1087-1095.

[8] Ortolan F, Fonini L S, Pastori T,. Tightly controlled expression of OsbHLH35 is critical for anther development in rice [J]., 2021, 302: 110716.

[9] 鄭國(guó)琦. 寧夏枸杞果實(shí)結(jié)構(gòu)、發(fā)育與糖分積累關(guān)系研究 [D]. 西安: 西北大學(xué), 2011.

[10] 李捷, 崔永濤, 柏延文, 等. 兩種枸杞對(duì)干旱脅迫的生理響應(yīng)及抗旱性評(píng)價(jià) [J]. 甘肅農(nóng)業(yè)大學(xué)學(xué)報(bào), 2019, 54(5): 79-87.

[11] 徐蕊. 枸杞LbMYB103基因在植物花器官發(fā)育中的功能初探 [D]. 銀川: 寧夏大學(xué), 2018.

[12] Grabherr M G, Haas B J, Yassour M,. Full-length transcriptome assembly from RNA-Seq data without a reference genome [J]., 2011, 29(7): 644-652.

[13] 葉國(guó)盛. 枸杞特優(yōu)品種寧1號(hào) [J]. 農(nóng)技服務(wù), 1996, (9): 24.

[14] 徐青, 秦墾, 馮愛玲, 等. 枸杞不育系與可育系花藥發(fā)育的細(xì)胞學(xué)研究 [J]. 寧夏大學(xué)學(xué)報(bào): 自然科學(xué)版, 2009, 30(3): 263-267.

[15] Jacobowitz J R, Doyle W C, Weng J K. PRX9 and PRX40 are extensin peroxidases essential for maintaining tapetum and microspore cell wall integrity duringanther development [J]., 2019, 31(4): 848-861.

[16] Hecht V, Vielle-Calzada J P, Hartog M V,. TheSOMATIC EMBRYOGENESIS RECEPTOR KINASE 1 gene is expressed in developing ovules and embryos and enhances embryogenic competence in culture [J]., 2001, 127(3): 803-816.

[17] Li Y, Kim J I, Pysh L,. Four isoforms of4-coumarate: CoA ligase have overlapping yet distinct roles in phenylpropanoid metabolism [J]., 2015, 169(4): 2409-2421.

[18] Kumar A, Ellis B E. The phenylalanine ammonia-lyase gene family in raspberry. Structure, expression, and evolution [J]., 2001, 127(1): 230-239.

[19] Banilas G, Karampelias M, Makariti I,. The olive DGAT2 gene is developmentally regulated and shares overlapping but distinct expression patterns with DGAT1 [J]., 2011, 62(2): 521-532.

Anther transcriptome and screening of genes related to anther development in

FAN Wen-qiang1, LIU Xue-xia1, MA Xiao-lan1, WANG Lu-yao1, ZHANG Xing1, YUE Si-jun1, TANG Jian-ning2, ZHU Jin-zhong3, YAO Xin-ling1, ZHENG Rui1

1. Key Laboratory of Protection and Utilization of Special Biological Resources in Western China, Ministry of Education, Key Laboratory of Modern Molecular Breeding for Dominant and Special Crops in Ningxia, College of Life Science, Ningxia University, Yinchuan 750021, China 2.Wolfberry Industry Development Center, Yinchuan 750021, China 3. Qixin Wolfberry Seedling Professional Cooperatives of Zhongning County, Zhongning 755100, China

To investigate the anther development process of, and explore the genes related.The anther transcriptome analysis was performed by the high-throughput Illumina HiSeq2500 sequencing platform using anthers of single microspore and bi-cellular pollen stage and mature pollen stage as research materials.A total of 128,151 unigene were obtained after assembly, and 40,140 unigenes were annotated by searching databases of NR, Swiss-Prot, GO, COG, KOG and KEGG. Under the condition of FDR (false discovery rate) ＜0.05 and FC (fold change) ≥2, a total of 46 differentially expressed genes were found between single microspore and bi-cellular pollen stage and mature pollen stage. KEGG annotations identified four differentially expressed genes involved in four pathways: glycerol ester metabolism, linoleic acid metabolism, phenylalanine metabolism and phenylpropanoid biosynthesis.There were eleven differentially expressed genes were involved in the anther development process, such as leucine-rich repeat protein (LRR), cyclin (Cyclin), receptor-like protein kinase (RPK), E3 ubiquitin protein ligase (MBR1), and so on. Five of the different expression genes among the eleven genes were randomly selected for qRT-PCR verification, and the expression results were basically consistent with the transcriptome data, indicating that the data were true and reliable.This study lays the foundation for further explaining the molecular mechanisms related to anther development in., and also provides useful genetic resources for subsequent studies on the biological characteristics of anther development in.

L.; anther development; transcriptome; differentially expressed genes; biological characteristics

R282.12

A

0253 - 2670(2022)03 - 0827 - 08

10.7501/j.issn.0253-2670.2022.03.023

2021-10-06

寧夏自然科學(xué)基金資助項(xiàng)目（2020AAC03097）；寧夏枸杞產(chǎn)業(yè)發(fā)展中心項(xiàng)目“寧夏枸杞核心種質(zhì)資源差異性分析及功能基因挖掘”；寧夏回族自治區(qū)重點(diǎn)研發(fā)計(jì)劃（重大）重點(diǎn)項(xiàng)目（2019BBF02022）

范文強(qiáng)，男，碩士，研究方向?yàn)橹参锷飳W(xué)。E-mail: 1534204028@qq.com

鄭 蕊，女，博士，研究方向?yàn)橹参锷锛夹g(shù)。E-mail: xlzheng@126.com

[責(zé)任編輯 時(shí)圣明]