摘 要：以寧夏賀蘭縣地方草石蠶（Stachys sieboldii）品種為試材，探討應(yīng)用第三代測序技術(shù)獲得草石蠶全長轉(zhuǎn)錄本信息，應(yīng)用第二代測序技術(shù)獲得3個不同發(fā)育階段草石蠶葉片和塊莖的轉(zhuǎn)錄組信息，對測序結(jié)果進(jìn)行轉(zhuǎn)錄組水平分析，篩選特有差異基因，并進(jìn)行GO和KEGG富集分析，開展草石蠶基因表達(dá)特性的初步研究。結(jié)果表明，第三代測序后Polymerase read的數(shù)據(jù)量為50.82 G，F(xiàn)LNC序列的reads數(shù)為525 593個轉(zhuǎn)錄本；在KEGG等7大數(shù)據(jù)庫的基因功能注釋中均注釋成功的轉(zhuǎn)錄本數(shù)目為6 857個，至少有1個數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為14 078個；與NR數(shù)據(jù)庫比對注釋后，草石蠶與同為唇形目的芝麻（Sesamum indicum）基因序列相似性最高，相似基因個數(shù)為9 149個；與GO數(shù)據(jù)庫比對注釋后，生物學(xué)過程、細(xì)胞成分與分子功能中注釋到基因個數(shù)最多的分別是新陳代謝過程5 093個、細(xì)胞2 004個、鍵聯(lián)結(jié)合6 645個；與KEGG數(shù)據(jù)庫比對注釋后，在細(xì)胞轉(zhuǎn)化、環(huán)境信息處理、遺傳信息處理、新陳代謝和有機(jī)系統(tǒng)功能中注釋到基因數(shù)最多的分別是運輸和分解代謝409個、信號傳導(dǎo)729個、轉(zhuǎn)化626個、碳水化合物代謝601個、內(nèi)分泌系統(tǒng)297個。

第二代測序后，ZLZ_2_S_3和ZLZ_2_S_1差異基因比較中的上調(diào)基因數(shù)最高，為2 303個；ZLZ_2_L_3和ZLZ_2_S_3差異比較中的下調(diào)基因數(shù)最高，為2 033個。葉片3個時期ZLZ_2_L_1、ZLZ_2_L_2與ZLZ_2_L_3組合之間比較差異基因總數(shù)為13 610個，共有差異基因數(shù)為10 203個；ZLZ_2_L_3獨有的差異基因數(shù)最多，為437個。塊莖3個時期ZLZ_2_S_1、ZLZ_2_S_2與ZLZ_2_S_3的組合之間比較差異基因總數(shù)為13 732個，共有的差異基因數(shù)為11 370個，ZLZ_2_S_3獨有的差異基因數(shù)最多，為412個。ZLZ_2_S_3、ZLZ_2_S_2聚類圖顯示高表達(dá)基因的數(shù)值范圍主要為0～2，顯著高于其他處理，可以聚為一類。3個不同時期葉片的差異基因主要集中在光合作用生物碳固定、內(nèi)質(zhì)網(wǎng)中蛋白質(zhì)加工以及乙醛酸和二羧酸代謝等通路上；3個不同時期塊莖的差異基因主要集中在淀粉和蔗糖代謝、糖解與糖代謝合成等通路上。以上結(jié)論將為今后研究草石蠶的生物學(xué)特性、闡明特有表型差異機(jī)制提供參考，為提升草石蠶的基礎(chǔ)理論研究水平提供技術(shù)支撐。

關(guān)鍵詞：第二代轉(zhuǎn)錄組； 第三代全長轉(zhuǎn)錄組； 草石蠶； 基因表達(dá)特性

中圖分類號：S644.5" " " " 文獻(xiàn)標(biāo)識碼：A" " 文章編號：1002-204X（2023）05-0016-14

doi：10.3969/j.issn.1002-204x.2023.05.004

Preliminary Study on Gene Expression Characteristics of Stachys sieboldii

based on Second-Generation and Third-Generation Whole Transcriptome Sequencing Conditions

Li Cheng1， Pei Hongxia1， Qu Ling2*， Guo Song1，

Wang Zhiqiang1， Yang Wanbang1， Liang Pu3， Zhang Guifang4

（1.Horticulture Institute， Ningxia Academy of Agriculture and Forestry Sciences， Yinchuan， Ningxia 750002；

2.Institute of" Wolfberry" Science， Ningxia Academy of Agriculture and Forestry Sciences， Yinchuan， Ningxia 750002； 3.Ningxia Horticulture Extension Station， Yinchuan， Ningxia 752399； 4.Xingqing District Agriculture and Animal Husbandry Bureau， Yinchuan， Ningxia 750002 ）

Abstract Stachys sieboldii， a local variety from Helan County in Ningxia， was used as the experimental material to explore the application of third-generation sequencing technology to obtain the full-length transcriptome information of Stachys sieboldii. Second-generation sequencing technology was used to obtain the transcriptome information of leaves and tubers at three different developmental stages of Stachys sieboldii. The sequencing results were analyzed at the transcriptome level， specific differentially expressed genes were screened， and GO and KEGG enrichment analyses were conducted to conduct a preliminary study on the gene expression characteristics of Stachys sieboldii. The results showed that the data volume of Polymerase read after third-generation sequencing was 50.82 G， and the number of reads of FLNC sequences was 525 593 transcripts. Among the gene function annotations in the seven major databases， the number of annotated transcripts was 6 857， and the number of transcripts with at least one successful annotation in the database was 14 078. After aligning with the NR database， the gene sequence of Stachys sieboldii showed the highest similarity to that of Sesamum indicum， which also belongs to the Lamiales order， with 9 149 similar genes. After aligning with the GO database， the numbers of genes annotated in biological processes， cellular components， and molecular functions were 5 093， 2 004 and 6 645， respectively. After aligning with the KEGG database， the numbers of genes annotated in cell transformation， environmental information processing， genetic information processing， metabolism， and organic system functions were 409， 729， 626， 601 and 297， respectively.

After second-generation sequencing， the highest number of upregulated genes in the comparison between ZLZ_2_S_3 and ZLZ_2_S_1 was 2 303. The highest number of downregulated genes in the comparison between ZLZ_2_L_3 and ZLZ_2_S_3 was 2 033. The total number of differentially expressed genes between the three stages of leaves （ZLZ_2_L_1， ZLZ_2_L_2， and ZLZ_2_L_3） was 13 610， with a total of 10 203 differentially expressed genes. ZLZ_2_L_3 had the highest number of unique differentially expressed genes， with 437. The total number of differentially expressed genes between the three stages of tubers （ZLZ_2_S_1， ZLZ_2_S_2， and ZLZ_2_S_3） was 13 732， with a total of 11 370 differentially expressed genes. ZLZ_2_S_3 had the highest number of unique differentially expressed genes， with 412. Cluster analysis showed that ZLZ_2_S_3 and ZLZ_2_S_2 had a higher expression level of genes mainly distributed between 0 and 2， significantly higher than other treatments， and could be clustered into one category. The differentially expressed genes in the three stages of leaves mainly concentrated in pathways such as photosynthesis， carbon fixation in photosynthetic organisms， protein processing in the endoplasmic reticulum， and glyoxylate and dicarboxylate metabolism. The differentially expressed genes in the three stages of tubers mainly concentrated in pathways such as starch and sucrose metabolism， glycolysis/gluconeogenesis， and pyruvate metabolism. These conclusions will provide references for future studies on the biological characteristics of Stachys sieboldii and elucidate the mechanisms of unique phenotypic differences， thereby providing technical support for the improvement of the basic theoretical research level of Stachys sieboldii.

Key words Second-generation transcriptome; Third-generation whole transcriptome; Stachys sieboldii; Gene expression characteristics

唇形科水蘇屬中的草石蠶是多年生草本植物，也是產(chǎn)品為地下塊莖的特色蔬菜，因地下匍匐枝在成熟時，頂端膨大而形成螺旋狀肉質(zhì)塊莖，別名螺絲菜、甘露子、寶塔菜、地蠶、地溜子等。草石蠶起源于我國，唐代《本草拾遺》始稱之為“草石蠶”。其塊莖營養(yǎng)豐富，有特殊香味，可加工成醬腌漬品等，是馳名中外的“八寶菜”之一；同時也可用于提取水蘇糖，具有良好的潤肺益腎等藥用功能[1-4]。草石蠶在寧夏已有近40年的種植歷史，特別是寧夏賀蘭縣螺絲菜作為“寧”字號道地蔬菜享譽(yù)國內(nèi)外。在全國草石蠶產(chǎn)業(yè)中寧夏屬重要的核心產(chǎn)區(qū)，種植面積多年來穩(wěn)定在2 000 hm2以上，形成“老”字號的外向型特色蔬菜產(chǎn)業(yè)，種植和加工總產(chǎn)值過億元。

采用田間壟式栽培對賀蘭縣地方草石蠶品種的植株和塊莖營養(yǎng)體繁殖后獲得試驗材料，基于第三代全長轉(zhuǎn)錄組技術(shù)獲得高質(zhì)量草石蠶不同發(fā)育階段或組織部位的全長轉(zhuǎn)錄本信息，基于第二代測序技術(shù)獲得不同品種（系）或不同發(fā)育階段草石蠶植株和塊莖的轉(zhuǎn)錄組信息，并與第三代測序結(jié)果進(jìn)行轉(zhuǎn)錄組水平的比對。對篩選出的特有差異基因進(jìn)行GO富集和KEGG富集分析，開展基于第二代轉(zhuǎn)錄組、第三代全長轉(zhuǎn)錄組測序條件下草石蠶基因表達(dá)特性的初步研究，以期為今后探究草石蠶的生物學(xué)特性，闡明草石蠶種質(zhì)的特有表型差異機(jī)制提供參考，為提升寧夏“寧”字號特色蔬菜草石蠶的基礎(chǔ)理論研究水平，建設(shè)好種質(zhì)資源創(chuàng)新平臺提供一定的技術(shù)支撐。

本試驗于2020年3月中旬在賀蘭縣開始田間種植試驗，中后期在寧夏農(nóng)林科學(xué)院園藝研究所試驗室、北京諾禾致源科技股份有限公司試驗室進(jìn)行采樣、轉(zhuǎn)錄組測序及分析，2021年5月結(jié)束。

1 材料與方法

1.1 試驗材料

供試草石蠶為寧夏賀蘭縣農(nóng)家品種，由寧夏農(nóng)林科學(xué)院園藝研究所收集。

1.2 試驗方法

2020年3月18日，本試驗在銀川市賀蘭縣四十里店光明村示范戶的種植基地，采用壟式栽培方式在田間定植草石蠶塊莖后，實現(xiàn)了草石蠶的種植和塊莖繁殖[5]。本試驗分為3個時期。第一個時期：花期，匍匐莖初生（8月6日采樣）；第二個時期：地下莖伸長，塊莖初生時期，9月2日采樣；第三個時期：塊莖生長最旺盛時期，10月1日采樣，完成對草石蠶葉片和塊莖等試驗樣品的3個批次田間采集，取樣部位包括葉片組織和地下塊莖共2個部位，每個組織取樣重復(fù)3次，共取樣18份。每個樣取500 mg新鮮組織。18個標(biāo)號包括第一個時期的L-1-1、S-1-1、L-1-2、S-1-2、L-1-3、S-1-3，第二個時期的L-2-1、S-2-1、L-2-2、S-2-2、L-2-3、S-2-3，第三個時期的L-3-1、S-3-1、L-3-2、S-3-2、L-3-3、S-3-3。其中：在L-1-1標(biāo)號中，“L”代表葉片，第一個“-1”代表在第一個時期采樣，第二個“-1”代表第1株樣本。在S-3-3標(biāo)號中，“S”代表地下塊莖，第一個“-3”代表在第三個時期采樣，第二個“-3”代表第三株樣本，以此類推。在分別取得不同處理的樣本后，用凍存管存放并標(biāo)記，分組分樣快速置于-80 ℃液氮中冷凍保存；隨后，在寧夏農(nóng)林科學(xué)院園藝研究所試驗室整理裝入自封袋，采用干冰保存法運輸?shù)奖本┲Z禾致源科技股份有限公司分子試驗室，經(jīng)過所有樣品總RNA提取和4個步驟檢驗合格后，開展第二代轉(zhuǎn)錄組、第三代全長轉(zhuǎn)錄組相關(guān)測序。

1.3 測定項目及方法

該試驗中，提取各樣本的總RNA后，分別檢驗，檢驗合格后，應(yīng)用PacBio Sequel測序系統(tǒng)（第三代測序）實現(xiàn)單分子實時測序，通過對草石蠶3個不同發(fā)育階段（花期或匍匐莖初生、塊莖初生時期和塊莖生長最旺盛時期）和2個不同組織部位（葉片、塊莖），共6個獨立樣本混樣后的總樣本開展第三代測序，以期獲得草石蠶全長轉(zhuǎn)錄本信息，對數(shù)據(jù)量的要求是30 G/樣以上。

應(yīng)用Illumina Hiseq測序系統(tǒng)（第二代測序）獲得3個不同發(fā)育階段（3個時期）和2個不同組織部位（葉片、塊莖），共18個草石蠶獨立樣本的轉(zhuǎn)錄組信息，對數(shù)據(jù)量的要求是6 G/樣。通過對試驗結(jié)果的建庫上機(jī)分析后，對第二、三代測序結(jié)果進(jìn)行轉(zhuǎn)錄組水平的分析和比對，初步進(jìn)行基因差異表達(dá)分析，篩選特有差異基因，對這些差異基因進(jìn)行GO富集和KEGG富集分析，開展了草石蠶基因表達(dá)特性的初步研究[6-11]。

2 結(jié)果與分析

2.1 草石蠶第三代全長轉(zhuǎn)錄本的測序結(jié)果與分析

2.1.1 consensus統(tǒng)計

PacBio Sequel測序平臺為環(huán)狀測序，測序過程中單分子產(chǎn)出的高質(zhì)量測序reads稱為Polymerase read（酶聚合序列）。在本試驗中，草石蠶第三代測序后Polymerase read的數(shù)據(jù)量為50.82 G，Polymerase read的數(shù)目為814 910個。為了保證數(shù)據(jù)分析的質(zhì)量及可靠性，需要對原始數(shù)據(jù)再過濾，去除接頭和長度小于50 bp的原始下機(jī)數(shù)據(jù)，得到Subreads序列。本試驗中有效插入片段Subreads的數(shù)據(jù)量大小為49.38 G；隨后開展CCS統(tǒng)計，CCS序列是通過每個ZMW（零模波導(dǎo)孔）孔中的（Subreads）子序列經(jīng)過自身的比對校正后得到的一致性序列，無需進(jìn)行參考序列的比對。CCS序列要求每個插入序列中至少包含1個完整（full-pass）的Subreads，在本試驗中，通過聚類之后得到的CCS序列的reads數(shù)為692 084個轉(zhuǎn)錄本。由SMRT Analysis軟件所定義的兩端同時含有3'引物、5'引物，以及3'引物前含有polyA尾的序列統(tǒng)稱為全長序列（Full-Length （FL） Read）；反之，則為非全長序列（non-full-length read）。全長序列中的非嵌合序列稱為全長非嵌合序列（Full-Length-Non-Chimeric Read， FLNC），聚類之后得到的FLNC序列的reads數(shù)為525 593個轉(zhuǎn)錄本，F(xiàn)LNC數(shù)目占CCS數(shù)目的比例為75.94%。最后，全長轉(zhuǎn)錄組使用hierarchical n*log（n）算法將同一轉(zhuǎn)錄本的FLNC序列進(jìn)行聚類去冗余，得到Consensus序列，再利用arrow軟件對得到的Consensus序列進(jìn)行校正，最終獲得Polished consensus序列進(jìn)行后續(xù)分析，統(tǒng)計結(jié)果見圖1。從圖1中可以看出，聚類之后得到的Polished consensus reads數(shù)為29 889個轉(zhuǎn)錄本，其中：最短的Polished consensus長度為54 bp，最長的Polished consensus長度為10 680 bp，Polished consensus的平均長度為2 459 bp。

2.1.2 ZLZ_3_LS的基因功能注釋

2.1.2.1 ZLZ_3_LS大數(shù)據(jù)庫中的基因功能注釋 為獲得全面的基因功能信息，對使用CD-HIT軟件去冗余之后的序列進(jìn)行基因功能注釋，所使用的7大數(shù)據(jù)庫包括NR、NT、Pfam、KOG、SwissProt、KEGG、GO。這7種數(shù)據(jù)庫注釋結(jié)果如圖2所示。

NR數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為13 954個，SwissProt數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為12 382個，KEGG數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為13 354個，KOG數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為9 001個，GO數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為10 498個，NT數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為12 367個，Pfam數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為10 498個，在所有數(shù)據(jù)庫中均注釋成功（all Databases）的轉(zhuǎn)錄本數(shù)目為6 857個，至少有1個數(shù)據(jù)庫注釋成功（at least one Database）的轉(zhuǎn)錄本數(shù)目為14 078個。

此外，在7大數(shù)據(jù)庫注釋結(jié)果中選出5個數(shù)據(jù)庫繪制的維恩圖如圖3所示。

僅由NR數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為98個，僅由KEGG數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為3個，僅由KOG數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為1個，僅由GO數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為32個，僅由NT數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為68個，5個數(shù)據(jù)庫中都注釋成功（at least one Database）的轉(zhuǎn)錄本數(shù)目為6 942個。

2.1.2.2 ZLZ_3_LS數(shù)據(jù)庫的注釋和分類 NR是一個非冗余的蛋白質(zhì)數(shù)據(jù)庫，由NCBI創(chuàng)建并維護(hù)，其特點在于內(nèi)容比較全面，同時注釋結(jié)果中包含有物種信息，可作物種分類用。通過與NR數(shù)據(jù)庫進(jìn)行比對注釋后，可以獲取草石蠶基因序列與近緣物種基因序列的相似性，以及草石蠶基因的功能信息。統(tǒng)計注釋到的物種及基因數(shù)目統(tǒng)計結(jié)果如圖4所示。

草石蠶與同為唇形目的近緣物種芝麻（Sesamum indicum）的基因序列相似性最高，比對后相似基因個數(shù)為9 149個；與唇形目的斑點猴面花（Erythranthe guttata）的基因序列相似性較高，比對后相似基因個數(shù)為2 501個；與女貞病毒（Ligustrum virus）、旋蒴苣苔（Dorcoceras hygrometricum）和蠶豆（Broad bean）的基因序列相似的基因個數(shù)分別為370、222、111個，特別是注釋出女貞病毒（Ligustrum virus）。草石蠶因病毒病發(fā)病造成種性發(fā)生退化的研究目前基本上是空白，通過第三代全長轉(zhuǎn)錄組測序方法發(fā)現(xiàn)的這一情況可能會為今后的研究提供新的思路。其他的10余個物種的基因序列相似的基因個數(shù)都在100個以內(nèi)，相對較少。

GO是基因功能描述的分類系統(tǒng)，GO分為3大類：①細(xì)胞組分（Cellular Component），用于描述亞細(xì)胞結(jié)構(gòu)、位置和大分子復(fù)合物，如核仁、端粒和識別起始的復(fù)合物；②分子功能（Molecular Function），用于描述基因、基因產(chǎn)物個體的功能，如與碳水化合物結(jié)合或ATP水解酶活性等；③生物過程（Biological Process），用來描述基因編碼的產(chǎn)物所參與的生物過程，如有絲分裂或嘌呤代謝等。對基因進(jìn)行GO注釋之后，將注釋成功的基因按照GO 3大類的下一層級進(jìn)行分類，結(jié)果見圖5。

GO term的3種基本分類中的生物學(xué)過程中注釋到metabolic process（新陳代謝過程）的基因個數(shù)為5 093個，注釋到cellular process（細(xì)胞過程）的基因個數(shù)為4 757個，注釋到single-organism process（單有機(jī)體過程）的基因個數(shù)為2 854個，注釋到localization（定位）的基因個數(shù)為1 435個，注釋到biological regulation（生物調(diào)節(jié)）的基因個數(shù)為1 341個，注釋到regulation of biological process（生物過程調(diào)節(jié)）的基因個數(shù)為1 285個；細(xì)胞成分中注釋到cell（細(xì)胞）的基因個數(shù)為2 004個，注釋到cell part（細(xì)胞區(qū)域）的基因個數(shù)為2 004個，注釋到organelle（細(xì)胞器）的基因個數(shù)為1 425個，注釋到membrane（細(xì)胞膜）的基因個數(shù)為1 150個，注釋到macromolecular complex（大分子絡(luò)合物）的基因個數(shù)為1 112個，注釋到membrane part（薄膜部分）的基因個數(shù)為1 081個；分子功能中主要注釋到binding（鍵聯(lián)結(jié)合）的基因個數(shù)為6 645個，注釋到catalytic activity（催化活性）的基因個數(shù)為5 108個。

KEGG系統(tǒng)是分析基因產(chǎn)物和化合物在細(xì)胞中的代謝途徑及這些基因產(chǎn)物的功能的數(shù)據(jù)庫。其整合了基因組、化學(xué)分子和生化系統(tǒng)等方面的數(shù)據(jù)，KO（KEGG ORTHOLOG）系統(tǒng)將各個KEGG注釋系統(tǒng)聯(lián)系在一起，KEGG已建立了一套完整KO注釋的系統(tǒng)，可完成新測序物種的基因組或轉(zhuǎn)錄組的功能注釋?；蜃鱇O注釋后，可根據(jù)它們參與的KEGG代謝通路進(jìn)行分類，結(jié)果如圖6所示。

在Cellular Processes（細(xì)胞轉(zhuǎn)化）功能中，主要注釋到Transport and catabolism（運輸和分解代謝）的基因數(shù)為409個，注釋到Cell growth and death（細(xì)胞生長和死亡）的基因數(shù)為203個；在Environmental Information Processing（環(huán)境信息處理）功能中，主要注釋到Signal transduction（信號傳導(dǎo)）的基因數(shù)為729個；在Genetic Information Processing（遺傳信息處理）功能中，主要注釋到Translation（轉(zhuǎn)化）的基因數(shù)為626個，注釋到Folding， sorting and degradation（折疊、分揀和降解）的基因數(shù)為447個，注釋到Transcription（轉(zhuǎn)錄）的基因數(shù)為283個；在Human Diseases（人類疾病）功能中，主要注釋到Infectious diseases： Viral（病毒性傳染?。┑幕驍?shù)為364個，注釋到Cancers： Overview（癌癥綜述）的基因數(shù)為299個，注釋到Infectious diseases： Bacterial（細(xì)菌性感染）的基因數(shù)為257個；在Metabolism（新陳代謝）功能中，主要注釋到 Carbohydrate metabolism（碳水化合物代謝）的基因數(shù)為601個，注釋到Global and overview maps（全局和總覽圖）的基因數(shù)為414個，注釋到Energy metabolis（能量代謝）的基因數(shù)為387個，注釋到Amino acid metabolism（氨基酸代謝）的基因數(shù)為381個，注釋到Lipid metabolism（脂類代謝）的基因數(shù)為249個；在Organismal Systems（有機(jī)系統(tǒng)）功能中，主要注釋到Endocrine system（內(nèi)分泌系統(tǒng)）的基因數(shù)為297個，注釋到Nervous system（神經(jīng)系統(tǒng)）的基因數(shù)為239個，注釋到Immune system（免疫系統(tǒng)）的基因數(shù)為236個，注釋到Environmental adaptation（環(huán)境適應(yīng)）的基因數(shù)為219個。

2.1.3 ZLZ_3_LS轉(zhuǎn)錄因子分析

轉(zhuǎn)錄因子（TF）是一群能與基因5'端上游特定序列專一性結(jié)合，從而保證目的基因以特定的強(qiáng)度在特定的時間與空間表達(dá)的蛋白質(zhì)分子。使用iTAK軟件進(jìn)行植物轉(zhuǎn)錄因子預(yù)測。將注釋到轉(zhuǎn)錄本數(shù)目最多的前29個轉(zhuǎn)錄因子家族進(jìn)行柱形圖分析（圖7），包括C3H、WRKY、bHLH、SNF2、bZIP、C2H2、Others、MYB-related、NAC、AP2/ERF-ERF、FAR1、GRAS、B3-ARF、TRAF、HB-HD-ZIP、HSF、PHD、SET、Jumonji、AUX/IAA、MYB、IWS1、GARP-G2-like、Trihelix、SBP、B3、TCP、mTERF、TUB等不同的轉(zhuǎn)錄因子家族，其TF數(shù)目分別為51、46、44、37、36、35、33、32、28、28、27、26、26、25、24、23、23、21、17、16、16、15、14、14、14、12、11、11、11個。

2.2 草石蠶第二代轉(zhuǎn)錄組測序的結(jié)果與分析

2.2.1 參考序列比對分析

應(yīng)用Illumina Hiseq測序系統(tǒng)（第二代測序）獲得3個不同發(fā)育階段（3個時期）和2個不同組織部位（葉、匍匐莖或塊莖），共6個處理（ZLZ_2_L_1、ZLZ_2_S_1、ZLZ_2_L_2、ZLZ_2_S_2、ZLZ_2_L_3、ZLZ_2_S_3），每個處理3次重復(fù)。通過對18個草石蠶樣本的轉(zhuǎn)錄組信息建庫，并上機(jī)分析后，利用CD-HIT軟件對校正后的consensus序列去冗余，得到的轉(zhuǎn)錄本作為基因的參考序列（ref），然后將Illumina測序得到的每個樣品的clean reads比對到ref上。在此過程中采用RSEM軟件比對，比對統(tǒng)計結(jié)果如表1所示。Total reads：測序序列經(jīng)過質(zhì)控后的數(shù)量統(tǒng)計，Total mapped：能比對到參考序列上測序序列的數(shù)量統(tǒng)計。ZLZ_2_L_1、ZLZ_2_S_1、ZLZ_2_L_2、ZLZ_2_S_2、ZLZ_2_L_3、ZLZ_2_S_3能比對到參考序列上測序序列的數(shù)量比例分別是92.27%、93.71%、94.89%、93.42%、90.26%和88.71%。

2.2.2 差異基因篩選分析

火山圖可直觀展示每個比較組合的差異基因分布情況，差異比較分析中的組合安排如表2所示?；鹕綀D中橫坐標(biāo)表示基因在處理組和對照組中的表達(dá)倍數(shù)變化（log2（fold change）），縱坐標(biāo)表示基因在處理組和對照組中表達(dá)差異的顯著性水平（-log10（padj）或-log10（pvalue））。有顯著性差異表達(dá)的基因用紅色點（上調(diào)）和綠色點（下調(diào)）表示，無顯著性差異表達(dá)的基因用藍(lán)色點表示?；鹕綀D結(jié)果見圖8。

組合1（ZLZ_2_L_1和ZLZ_2_L_2）差異比較中，上調(diào)基因數(shù)為856個，下調(diào)基因數(shù)為239個；組合2（ZLZ_2_L_2和ZLZ_2_L_3）差異比較中，上調(diào)基因數(shù)為405個，下調(diào)基因數(shù)為1 224個；組合3（ZLZ_2_L_3和ZLZ_2_L_1）差異比較中，上調(diào)基因數(shù)為661個，下調(diào)基因數(shù)為640個；組合4（ZLZ_2_S_1和ZLZ_2_S_2）差異比較中，上調(diào)基因數(shù)為571個，下調(diào)基因數(shù)為594個；組合5（ZLZ_2_S_2和ZLZ_2_S_3）差異比較中，上調(diào)基因數(shù)為161個，下調(diào)基因數(shù)為355個；組合6（ZLZ_2_S_3和ZLZ_2_S_1）差異比較中，上調(diào)基因數(shù)為2 303個，下調(diào)基因數(shù)為617個；組合7（ZLZ_2_L_1和ZLZ_2_S_1）差異比較中，上調(diào)基因數(shù)為833個，下調(diào)基因數(shù)為378個；組合8（ZLZ_2_L_2和ZLZ_2_S_2）差異比較中，上調(diào)基因數(shù)為665個，下調(diào)基因數(shù)為1 750個；組合9（ZLZ_2_L_3和ZLZ_2_S_3）差異比較中，上調(diào)基因數(shù)為913個，下調(diào)基因數(shù)為2" 033個。綜合看，組合6（ZLZ_2_S_3和ZLZ_2_S_1）差異比較中的上調(diào)基因數(shù)最高，為2 303個；組合9（ZLZ_2_L_3和ZLZ_2_S_3）差異比較中的下調(diào)基因數(shù)最高，為2" 033個。

2.2.3 差異基因維恩圖分析

韋恩圖可展示不同比較組合間差異基因的重疊情況，通過韋恩圖，可篩選比較組合共有或獨有的差異基因。草石蠶3個不同發(fā)育階段（3個時期）和2個不同組織部位（葉片、塊莖），共6個處理（ZLZ_2_L_1、ZLZ_2_S_1、ZLZ_2_L_2、ZLZ_2_S_2、ZLZ_2_L_3、ZLZ_

2_S_3）間的差異基因維恩圖分析見圖9所示。通過分析，草石蠶葉片3個時期ZLZ_2_L_1、ZLZ_2_L_2與ZLZ_2_L_3的組合之間比較差異基因總數(shù)為13 610個，3個組合間共有的差異基因數(shù)為10 203個，ZLZ_2_L_1獨有的差異基因數(shù)為344個，ZLZ_2_L_2獨有的差異基因數(shù)為115個，ZLZ_2_L_3獨有的差異基因數(shù)為437個?？傮w上ZLZ_2_L_3獨有的差異基因數(shù)最多，其次是ZLZ_2_L_1，最少的是ZLZ_2_L_2。草石蠶塊莖3個時期ZLZ_2_S_1、ZLZ_2_S_2與ZLZ_2_S_3的組合之間比較差異基因總數(shù)為13 732個，3個組合間共有的差異基因數(shù)為11 370個，ZLZ_2_S_1獨有的差異基因數(shù)為130個，ZLZ_2_S_2獨有的差異基因數(shù)為90個，ZLZ_2_S_3獨有的差異基因數(shù)為412個?？傮w上ZLZ_2_S_3獨有的差異基因數(shù)最多，其次是ZLZ_2_S_1，最少的是ZLZ_2_S_2。

2.2.4 差異基因聚類分析

采用層次聚類可以對基因的FPKM值（每千個堿基的轉(zhuǎn)錄每百萬映射讀取的片段）進(jìn)行聚類分析，對行（row）進(jìn)行均一化處理，將log10（FPKM+1）值進(jìn)行歸一化轉(zhuǎn)換并聚類，紅色表示高表達(dá)基因，藍(lán)色表示低表達(dá)基因，顏色從紅到藍(lán)，表示log10（FPKM+1）從大到小。熱圖中表達(dá)模式相近的基因或樣本被聚集在一起，每個方格中的顏色反映的不是基因表達(dá)值，而是表達(dá)數(shù)據(jù)的行進(jìn)行均一化處理后得到的數(shù)值，所以熱圖中的顏色需要橫向比較（同一基因在不同樣本中的表達(dá)情況）。草石蠶3個不同發(fā)育階段（3個時期）和2個不同組織部位（葉片、塊莖），共6個處理（ZLZ_2_L_1、ZLZ_2_S_1、ZLZ_2_L_2、ZLZ_2_S_2、ZLZ_2_L_3、ZLZ_2_S_3）間的聚類構(gòu)成整體FPKM層次聚類圖如圖10所示。

因為ZLZ_2_S_3的時期正好是草石蠶塊莖的采收期，聚類圖顯示高表達(dá)基因的數(shù)值范圍主要為0～2，顯著高于其他5個處理，同時說明ZLZ_2_S_3高表達(dá)基因的表達(dá)值顯著高于其他5個處理；ZLZ_2_S_2的時期是草石蠶塊莖的膨大期，高表達(dá)基因的數(shù)值范圍也主要為0～2，顯著高于其他4個處理，ZLZ_2_S_3和ZLZ_2_S_2可以聚為一類，ZLZ_2_L_1和ZLZ_2_L_3聚為一類，ZLZ_2_S_1和ZLZ_2_L_2聚為一類。

2.2.5 差異基因KEGG富集散點圖分析

散點圖是KEGG富集分析結(jié)果的圖形化展示方式。KEGG富集程度通過Rich factor、Qvalue和富集到此通路上的基因個數(shù)來衡量。其中：Rich factor指該pathway中富集到的差異基因個數(shù)（Sample number）與注釋基因個數(shù)（Background number）的比值。Rich factor越大，表示富集的程度越大。Qvalue是經(jīng)過多重假設(shè)檢驗校正之后的Pvalue，Qvalue的取值范圍為0～1，取值越接近于零，表示富集越顯著。點的大小表示此pathway中差異表達(dá)基因個數(shù)多少，而點的顏色對應(yīng)于不同的Qvalue范圍。不同處理間具體的差異基因KEGG富集散點圖如圖11所示。

在組合1中，ZLZ_2_L_1和ZLZ_2_L_2處理之間，在Protein processing in endoplasmic reticulum（內(nèi)質(zhì)網(wǎng)中蛋白質(zhì)加工）通路中的差異表達(dá)基因個數(shù)為43個，在Carbon fixation in photosynthetic organisms（光合作用生物的碳固定）通路中的差異表達(dá)基因個數(shù)為19個，在Glyoxylate and dicarboxylate metabolism（乙醛酸和二羧酸代謝）通路中的差異表達(dá)基因個數(shù)為16個，在Phenylpropanoid biosynthesis（苯丙素的生物合成）通路中的差異表達(dá)基因個數(shù)為16個，在Photosynthesis（光合作用）通路中的差異表達(dá)基因個數(shù)為11個，均呈現(xiàn)顯著富集；在組合2中，ZLZ_2_L_2和ZLZ_2_L_3處理之間，在Photosynthesis-Antenna proteins（光合作用-觸角蛋白）通路中的差異表達(dá)基因個數(shù)為11個，在Glycerolipid metabolism（新陳代謝）通路中的差異表達(dá)基因個數(shù)為17個，在Tryptophan metabolism（色氨酸代謝）通路中的差異表達(dá)基因個數(shù)為9個，在Lysine degradation（賴氨酸退化）通路中的差異表達(dá)基因個數(shù)為10個，均呈現(xiàn)顯著富集；在組合3中，ZLZ_2_L_3和ZLZ_2_L_1處理之間，在Carbon fixation in photosynthetic organisms（光合作用生物的碳固定）通路中的差異表達(dá)基因個數(shù)為25個，在Glyoxylate and dicarboxylate metabolism（乙醛酸和二羧酸代謝）通路中的差異表達(dá)基因個數(shù)為21個，在Protein processing in endoplasmic reticulum（內(nèi)質(zhì)網(wǎng)中蛋白質(zhì)加工）通路中的差異表達(dá)基因個數(shù)為35個，在Glycine， serine and threonine metabolism（甘氨酸、絲氨酸和蘇氨酸代謝）通路中的差異表達(dá)基因個數(shù)為18個，在Phenylpropanoid biosynthesis（苯丙素的生物合成）通路中的差異表達(dá)基因個數(shù)為20個，均呈現(xiàn)顯著富集；在組合4中，ZLZ_2_S_1和ZLZ_2_S_2處理之間，在Tyrosine metabolism（絡(luò)氨酸代謝）通路中的差異表達(dá)基因個數(shù)為17個，在Isoquinoline alkaloid biosynthesis（異喹啉生物堿生物合成）通路中的差異表達(dá)基因個數(shù)為13個，在Protein processing in endoplasmic reticulum（內(nèi)質(zhì)網(wǎng)中蛋白質(zhì)加工）通路中的差異表達(dá)基因個數(shù)為35個，在Phenylpropanoid biosynthesis（苯丙素的生物合成）通路中的差異表達(dá)基因個數(shù)為20個，均呈現(xiàn)顯著富集；在組合5中，ZLZ_2_S_2和ZLZ_2_S_3處理之間，在Starch and sucrose metabolism（淀粉和蔗糖代謝）通路中的差異表達(dá)基因個數(shù)為16個，Photosynthesis-Antenna proteins（光合作用-觸角蛋白）通路中的差異表達(dá)基因個數(shù)為5個，都呈現(xiàn)顯著富集；在組合6中，ZLZ_2_S_3和ZLZ_2_S_1處理之間，在Tyrosine metabolism（絡(luò)氨酸代謝）通路中的差異表達(dá)基因個數(shù)為28個，在Isoquinoline alkaloid biosynthesis（異喹啉生物堿生物合成）通路中的差異表達(dá)基因個數(shù)為17個，在Starch and sucrose metabolism（淀粉和蔗糖代謝）通路中的差異表達(dá)基因個數(shù)為41個，在Fatty acid degradation（脂肪酸降解）通路中的差異表達(dá)基因個數(shù)為17個，在Glycolysis/Gluconeogenesis（糖解與糖代謝合成）通路中的差異表達(dá)基因個數(shù)為36個，均呈現(xiàn)顯著富集；在組合7中，ZLZ_2_L_1和ZLZ_2_S_1處理之間，在Photosynthesis（光合作用）通路中的差異表達(dá)基因個數(shù)為23個，在Carbon fixation in photosynthetic organisms（光合作用生物的碳固定）通路中的差異表達(dá)基因個數(shù)為28個，在Glyoxylate and dicarboxylate metabolism（乙醛酸和二羧酸代謝）通路中的差異表達(dá)基因個數(shù)為22個，在Glycolysis/Gluconeogenesis（糖解與糖代謝合成）通路中的差異表達(dá)基因個數(shù)為27個，在Starch and sucrose metabolism（淀粉和蔗糖代謝）通路中的差異表達(dá)基因個數(shù)為28個，均呈現(xiàn)顯著富集；在組合9中，ZLZ_2_L_3和ZLZ_2_S_3處理之間，在Ribosome（核糖體）通路中的差異表達(dá)基因個數(shù)為129個，在Ribosome biogenesis in eukaryotes（核糖體生物發(fā)生）通路中的差異表達(dá)基因個數(shù)為41個，在Carbon fixation in photosynthetic organisms（光合作用生物的碳固定）通路中的差異表達(dá)基因個數(shù)為38個，在Phenylpropanoid biosynthesis（苯丙素的生物合成）通路中的差異表達(dá)基因個數(shù)為42個，在Starch and sucrose metabolism（淀粉和蔗糖代謝）通路中的差異表達(dá)基因個數(shù)為48個，均呈現(xiàn)顯著富集。從差異基因KEGG富集散點圖可以得出，草石蠶的3個時期、2個部位之間的差異基因富集情況各有特點。總體來說，3個不同時期的葉片部位差異基因主要集中在光合作用生物的碳固定、內(nèi)質(zhì)網(wǎng)中蛋白質(zhì)加工，以及乙醛酸和二羧酸代謝等通路上；3個不同時期的塊莖部位差異基因主要集中在淀粉和蔗糖代謝、糖解與糖代謝合成、絡(luò)氨酸代謝和異喹啉生物堿生物合成等通路上；葉片部位與塊莖部位之間的差異基因主要集中在光合作用生物的碳固定、淀粉和蔗糖代謝等通路上。

3 結(jié)論與討論

應(yīng)用PacBio Sequel測序系統(tǒng)實現(xiàn)單分子實時測序，通過對草石蠶3個不同發(fā)育階段（花期或匍匐莖初生、塊莖初生時期和塊莖生長最旺盛時期）和2個不同組織部位（葉片、塊莖）開展第三代測序，獲得的草石蠶全長轉(zhuǎn)錄本信息；應(yīng)用Illumina Hiseq測序系統(tǒng)獲得3個不同發(fā)育階段（3個時期）和2個不同組織部位（葉片、塊莖）的轉(zhuǎn)錄組信息，對第二代、第三代測序結(jié)果進(jìn)行轉(zhuǎn)錄組水平的分析和比對，初步進(jìn)行了基因差異表達(dá)分析，篩選特有差異基因，對這些差異基因進(jìn)行GO富集和KEGG富集分析，開展了草石蠶基因表達(dá)特性的初步研究。

結(jié)果表明，草石蠶第三代測序后Polymerase read的數(shù)據(jù)量為50.82 G，得到的FLNC序列的reads數(shù)為525 593個轉(zhuǎn)錄本；在KEGG、GO等7大數(shù)據(jù)庫的基因功能注釋中，在所有數(shù)據(jù)庫中均注釋成功的轉(zhuǎn)錄本數(shù)目為6 857個；至少有1個數(shù)據(jù)庫注釋成功的轉(zhuǎn)錄本數(shù)目為14 078個；與NR數(shù)據(jù)庫進(jìn)行比對注釋后，草石蠶與同為唇形目的近緣物種芝麻（Sesamum indicum Linn.）的基因序列相似性最高，比對后相似基因個數(shù)為9 149個；與GO數(shù)據(jù)庫比對注釋后，生物學(xué)過程、細(xì)胞成分與分子功能中注釋到基因個數(shù)最多的分別是metabolic process（新陳代謝過程）5 093個、cell（細(xì)胞）2 004個、binding（鍵聯(lián)結(jié)合）6 645個。與KEGG數(shù)據(jù)庫比對注釋后，在細(xì)胞轉(zhuǎn)化、環(huán)境信息處理、遺傳信息處理、新陳代謝和有機(jī)系統(tǒng)功能中注釋到基因數(shù)最多的分別是Transport and catabolism（運輸和分解代謝）409個、Signal transduction（信號傳導(dǎo)）729個、Translation（轉(zhuǎn)化）626個、Carbohydrate metabolism（碳水化合物代謝）601個、Endocrine system（內(nèi)分泌系統(tǒng)）297個。

第二代測序后，ZLZ_2_S_3（塊莖生長最旺盛時期的塊莖組織）和ZLZ_2_S_1（匍匐莖初生時期的塊莖組織）差異基因比較中的上調(diào)基因數(shù)最高，為2 303個；ZLZ_2_L_3（塊莖生長最旺盛時期的葉片組織）和ZLZ_2_S_3差異比較中的下調(diào)基因數(shù)最高，為2 033個。葉片3個時期ZLZ_2_L_1、ZLZ_2_L_2與ZLZ_2_L_3的組合之間比較差異基因總數(shù)為13 610個，3個組合間共有的差異基因數(shù)為10 203個，ZLZ_2_L_3獨有的差異基因數(shù)最多，為437個；塊莖3個時期ZLZ_2_S_1、ZLZ_2_S_2與ZLZ_2_S_3的組合之間比較差異基因總數(shù)為13 732個，3個組合間共有的差異基因數(shù)為11 370個，ZLZ_2_S_3獨有的差異基因數(shù)最多，為412個。ZLZ_2_S_3、ZLZ_2_S_2聚類圖顯示高表達(dá)基因的數(shù)值范圍主要為0～2，顯著高于其他處理的數(shù)值，可以聚為一類。3個不同時期葉片的差異基因主要集中在光合作用生物的碳固定、內(nèi)質(zhì)網(wǎng)中蛋白質(zhì)加工，以及乙醛酸和二羧酸代謝等通路上；3個不同時期塊莖的差異基因主要集中在淀粉和蔗糖代謝、糖解與糖代謝合成、絡(luò)氨酸代謝和異喹啉生物堿生物合成等通路上。以上結(jié)論將為繼續(xù)開展基于第二代轉(zhuǎn)錄組、第三代全長轉(zhuǎn)錄組測序條件下草石蠶基因表達(dá)特性的研究，對今后探究草石蠶的生物學(xué)特性，闡明草石蠶種質(zhì)的一些特有表型差異機(jī)制具有重要意義；可為提升寧夏“寧”字號特色蔬菜草石蠶的基礎(chǔ)理論研究水平，建設(shè)好種質(zhì)資源創(chuàng)新平臺提供技術(shù)支撐。

參考文獻(xiàn)：

[1] 張德純. 草石蠶[J]. 中國蔬菜，2015（6）：10.

[2] 鄭美娟，張洪，馬明霞，等. 草石蠶研究文獻(xiàn)綜述[J]. 廣東化工，2015，42（13）：169-170.

[3] 張振賢，喻景權(quán)，于賢昌，等. 蔬菜栽培學(xué)[M]. 北京：中國農(nóng)業(yè)大學(xué)出版社，2003.

[4] 耿廣東，楊祎，張素勤. 草石蠶塊莖膨大相關(guān)SSH文庫構(gòu)建及分析[J]. 分子植物育種，2020，18（17）：5650-5657.

[5] 李程，郭松，王志強(qiáng)，等. 草石蠶新型壟式復(fù)合栽培技術(shù)應(yīng)用及對產(chǎn)業(yè)可持續(xù)發(fā)展的初步探討[J]. 寧夏農(nóng)林科技，2019，60（3）：6-9.

[6] 王芬，木仁，張雯，等. 基于第三代納米孔測序技術(shù)的茶樹全長轉(zhuǎn)錄組分析[J]. 福建農(nóng)業(yè)學(xué)報，2021，36（10）：1177-1186.

[7] 周潤澤，王翊豪，崔志瑩，等. 基于PacBio的三代測序的三七基因組測序分析[J]. 時珍國醫(yī)國藥，2021，32（3）：719-723.

[8] 顏秀娟，何鑫，趙云霞，等. 銀川羊角椒雄性不育兩用系轉(zhuǎn)錄組與代謝組關(guān)聯(lián)分析[J]. 分子植物育種，2021，19（13）：4252-4258.

[9] 王瑞嫻，李川. 全長轉(zhuǎn)錄組測序技術(shù)在非模式植物轉(zhuǎn)錄組學(xué)研究中的應(yīng)用[J]. 分子植物育種，2019，17（2）：502-508.

[10] 朱菲. 鐵皮石斛三代基因組及其在功能基因與性狀相關(guān)基因挖掘中的應(yīng)用[D]. 南京：南京師范大學(xué)，2021.

[11] 鐘偉民，張興坦，趙茜，等. 三代測序PacBio在轉(zhuǎn)錄組研究中的應(yīng)用[J]. 福建農(nóng)林大學(xué)學(xué)報（自然科學(xué)版），2018，47（5）：524-529.