基于單分子納米孔測序技術揭示大白菜抽薹期全長轉錄組差異的復雜性

關鍵詞：第三代高通量測序技術；單分子納米孔測序技術；大白菜；抽薹期；全長轉錄本

中圖分類號：S634.1 文獻標識號：A 文章編號：1001-4942（2024）11-0001-12

春白菜未熟抽薹是指短縮莖因對低溫敏感從而通過春化在結球后過度伸長生長，致使葉球部分或全部喪失商品價值的現象，是困擾春白菜生產的主要難題。大白菜抽薹是多基因控制的數量性狀，且存在基因型與環境互作的復雜情況。擬南芥中，薹莖發育受6條開花途徑控制，包括受環境因子影響的春化途徑、光周期途徑、環境溫度途徑及受內部因子調控的自主途徑、赤霉素途徑和年齡途徑，這6條途徑的上游基因相互獨立，但下游最終會作用于兩個關鍵整合因子FT和SOCl。春化途徑通過調控因子VRN1～5、VIL1～4、PRC 2、LHP1等協同作用感受低溫后，以表觀調控的方式使FLC（Flower Locus C）沉默，進而解除對FT的抑制作用，其中長鏈非編碼RNA（IncRNA）也起到了關鍵作用：光周期途徑中，關鍵調控因子CO位于FT上游，其周期性變化受各種光受體和COP1/SPA-E3-泛素連接酶復合體的精細調控；在赤霉素途徑中，赤霉素合成與代謝以及由DELLA蛋白介導的赤霉素信號通路都會影響花芽分化和薹莖伸長。因此，植物的開花調控網絡是由多基因參與、多條途徑控制的復雜網絡。

測序技術的進步和大白菜全基因組參考序列的公布，為利用組學技術從全基因組水平揭示大白菜晚抽薹的復雜分子機制提供了新的研究策略。目前，大白菜抽薹相關組學研究已取得了一些進展，研究者們先后采用二代測序技術、甲基化測序揭示出大白菜在春化作用或者不同光周期處理下的轉錄組、甲基化和microRNA的差異。但前人研究所用的多是單一材料在可控條件下的轉錄組變化，而對于晚抽薹和早抽薹材料在接近自然生長條件下通過春化后的全長轉錄組差異分析鮮見報道。因此，本研究選擇大白菜晚抽薹材料06-247與早抽薹材料He102.采用單分子納米孔測序技術對兩者抽薹期帶葉片的花序軸進行全長轉錄組測序，挖掘兩個材料抽薹期全長轉錄組的結構變異和表達差異，并對兩者全長轉錄本之間存在的可變剪切（alternative splicing，AS）、可變多聚腺苷酸位點（alternative polyadenylation，APA）、長鏈非編碼RNA （long non-coding RNA，IncRNA）進行分析，旨在為揭示大白菜晚抽薹分子機制提供新的線索。

1材料與方法

1.1供試材料與樣品采集

所用試材為大白菜晚抽薹自交系06-247與早抽薹自交系He102，二者在抽薹性方面的差異主要表現為在相同條件下完成春化后抽薹期主花莖長度存在極顯著差異。材料種植于山東省農業科學院蔬菜研究所核心試驗基地。2018年12月31日，將經催芽后的兩材料單粒點播于含壤土：草炭：珍珠巖（1：1：1）的育苗缽（口徑8 cm）中，置于冬暖大棚中育苗并使其自然通過春化。冬暖大棚為雙層棚膜結構設計，內層棚膜外設有保溫被，當夜間氣溫低于-3℃時覆蓋保溫被過夜，自然光照。2019年3月1日將幼苗轉移至紗網棚煉苗，3月10日移栽至紗網棚。小區設置為0.8 mx3.0m，雙行種植，行距40cm，株距30 cm，每小區18～20株，3次重復。3月30日，從每個小區至少10個生長正常的單株上采集帶葉片的主花序軸，在低溫下剪碎、混勻，用于總RNA的提取。

1.2全長轉錄組測序

1.2.1總RNA提取及質控 總RNA提取采用Trizol法，用RNase-free DNase I（Ambion，美國）處理總RNA以消除潛在的基因組DNA污染，用Bioanalyzer 2100（Agilent，美國）評估總RNA質量，以RIN≥7且OD_260nm/OD_280nm值≥1.8為閾值作為判斷標準。自檢合格的RNA分成兩份，一份用于轉錄組測序，一份用于qRT-PCR驗證。

1.2.2全長cDNA文庫構建和測序 用于測序的3次生物學重復的總RNA樣本等量混勻，采用SMART技術構建全長cDNA文庫（SMARTer@PCR cDNA Synthesis Kit，Cat No：634925，TaKaRa生物科技（美國）有限公司生產），得到的cDNA添加switch oligo后再合成互補鏈。對合成的DNA進行損傷修復和末端修復，再利用磁珠對cDNA進行純化。最后委托北京百邁克生物科技有限公司對cDNA文庫進行全長轉錄組測序。測序平臺為Nanopore PromethION48（Oxford Nano-pore Technologies，簡稱ONT，英國），實驗按照ONT公司提供的標準流程執行。

1.3全長轉錄本鑒定與注釋

1.3.1原始讀段的初步處理 測序產生的原始下機序列稱為原始讀段（raw reads），用Min-KNOW2.2軟件包中的Guppy軟件進行base call-ing，過濾掉Qlt;7、長度lt;500 bp的低質量序列和核糖體RNA序列，得到有效讀段（clean reads），并根據序列兩端是否存在引物得到全長讀段（ full-length reads）。用Minmap2軟件對全長讀段與大白菜參考基因組V3.0進行比對，再使用pin-fish軟件分析得到一致性序列（consensus se-quence）。用一致性序列進行后續的所有分析。

1.3.2全長轉錄本以及新基因的篩選 合并每個樣品的一致性序列，再次通過Minimap2與參考基因組進行比對，去除冗余序列，過濾掉identitylt;0.9、coveragelt;0.85的序列，僅合并5’端外顯子有差異的序列，最終獲得高質量非冗余全長轉錄本序列（ collapsed sequence）。將全長轉錄本與大白菜參考轉錄本進行比對，篩選出已知基因或已知基因的異構體。對于不能比對到參考轉錄本中的全長序列，再采用TransDecoder（v3.0.0）軟件，基于開放閱讀框（open reading frame，ORF）長度、對數似然函數值（log-likelihood score）、氨基酸序列與Pfam數據庫蛋白質結構域序列的比對等信息，從中識別可靠的潛在編碼區序列（coding se-quence，CDS），視為具有編碼潛能的新基因。

1.3.3全長轉錄本的注釋 利用Blast工具將上述所有全長轉錄本在Pfam、KOG[24]、Swiss-Prot、NR、eggNOG、GO（Gene Ontolo-gy）和KEGG（Kyoto Encyclopedia of Genes and Genomes）數據庫進行分析，獲得相應的注釋信息。

1.4全長轉錄本結構變異分析

1.4.1可變剪切篩選 可變剪切（alternative spli-cing，AS）的類型主要包括外顯子互斥（mutually excluslve exons，MEE）、外顯子跳躍（skipping ex-on，SE）、內含子保留（retained intron，RI）、可變5′剪切位點（alternative 5′splice-site，A5）、可變3′剪切位點（alternative 3′splice-site，A3）。通過Asta-lavista軟件獲取每個樣品的可變剪切類型，并對各種AS事件類型進行統計。

1.4.2可變多聚腺苷酸位點分析 按下列方法對每個樣品進行可變多聚腺苷酸位點（alternative polyadenylation，APA）的篩選：比較冗余轉錄本的3′端非翻譯區序列（3′-UTR），根據3′-UTR是否完全一致判斷APA的數目。利用MEME比較冗余轉錄本polyA位點上游50 bp的序列，鑒別相應的polyA信號（PAS）。

1.5長鏈非編碼RNA預測

通過對所得非冗余轉錄本的編碼潛能進行分析，篩選出不編碼蛋白的轉錄本，即長鏈非編碼RNA（IncRNA）。再用編碼潛能計算器（CodingPotential Calculato，簡稱CPC）、編碼一非編碼指數分析工具（Coding-Non-Coding Inde，簡稱CP-CI）、編碼潛能評估工具（Coding Potential As-sessment Tool，簡稱CPAT）、蛋白結構域預測工具Pfam對IncRNA進行分類，取4種方法的交集作為高質量IncRNA。

采用兩種方法對IncRNA的靶基因進行預測。方法①：基于IncRNA調控其鄰近基因的表達，主要根據IncRNA與mRNA的位置關系進行預測，定義染色體中每100kb范圍內存在差異表達的IncRNA與差異表達的mRNA，得到順式作用（cis-acting）靶基因。方法②：lncRNA與mRNA由于堿基互補配對而發生作用，主要利用靶基因預測工具LncTar進行靶基因預測，得到反式調控靶基因。

1.6轉錄本定量分析及與抽薹相關的差異基因篩選

1.6.1轉錄本定量分析及差異表達轉錄本鑒定 轉錄組測序可以模擬成一個隨機抽樣的過程，為了讓片段數目能真實地反映轉錄本表達水平，首先對樣品中Mapped Reads的數目進行歸一化。采用CPM（counts per million） 作為衡量轉錄本或基因表達水平的指標，CPM計算公式如下：

CPM=N₁/N₂×1000000。

其中，CPM表示CPM；N₁表示reads mapped to tran-script，即比對到某一轉錄本上的reads數目；N₂表示total reads aligned in sample，即比對到參考轉錄組的片段總數。

對于沒有生物學重復的樣本，使用EBSeq進行差異分析，參數設置FC≥2且FDRlt;0.01。FC（Fold Change）表示差異倍數，即兩（組）樣品間表達量的比值。FDR是False Discovery Rate的簡稱，是通過對差異顯著性P值（P-value）進行校正得到的。

1.6.2與抽薹有關的差異候選基因篩選 在差異表達基因集的GO注釋中，用“vernalization”、“photoperiod”、“gibberellins”和“gibberellic acid”為關鍵詞進行基因條目的篩選，合并所有條目并去除重復項，獲得與抽薹直接相關的差異候選基因。

1.6.3差異候選基因qRT-PCR驗證 總RNA的逆轉錄采用天根的Quant cDNA第一鏈合成試劑盒（KR103）進行。qRT-PCR在Applied Biosys-tems QuantStudio 3型qPCR儀上進行，采用天根RealUniversal彩色熒光定量預混試劑（SYBR Green，FP201）并按說明書操作，所用引物見表1。反應總體積20μL，包含如下成分：2×RealUniver-sal PreMi×10μL，Primer F（10μmol·L^-1）0.6μL，Primer R（10μmol.L^-1） 0.6μL，0.1～100 ng cDNA模板，添加ddH2O至20μL。以葡萄糖-6-磷酸脫氫酶基因（G6PD）作為內參基因，按照2^-ΔΔCt）法計算特異基因相對表達量。統計分析用SAS Version 9.0軟件進行，依據Duncan’s multiple range test法計算0.05水平差異顯著性。

2結果與分析

2.1測序數據質控分析

He102與06-247分別得到6962774條和5659251條原始讀段，總堿基數分別為82.02Gb和62.64 Gb；最大長度分別是11902bp和10407bp，居中長度（N50）分別為1329bp和1211bp，平均長度分別是1178bp和1106bp（表2）。He102與06-247的原始讀段長度大部分介于1kb～10kb，其中分布reads數最多的長度均為1kb（圖1A、B）。原始讀段的質量（Q）值介于1～18之間，讀段的質量值分布最多的是11（圖1C、D），平均質量值為11（表2）。上述數據表明本次測序產生了高質量的數據。

2.2全長轉錄本的鑒定

He102與06-247測序產生的高質量有效讀段（clean reads）分別是6 896 814條和5586307條，兩側都有引物序列的全長讀段（full-length reads）分別有6051557條和4844987條，去冗余后分別得到64228條和55123條一致性序列（consensus isoforms）；過濾掉identitylt;0.9、cover-agelt;0.85的序列，合并僅5′端外顯子有差異的reads，兩者分別得到37309個和30415個高質量非冗余全長轉錄本，其居中長度分別為1628bp和1537bp，平均長度分別為1413bp和1337bp，最大長度分別為7074bp和5642bp（表3）。

2.3全長轉錄本序列的比對及新基因的數據庫注釋

整合各樣品高質量非冗余轉錄本，累計獲得46485個unigene，包括44797個已知基因和1688個新基因。經過嚴格的序列比對，發現22971個已知基因的序列與參考轉錄組同源基因的序列之間存在多態性差異，它們是已知基因的異構體（isoforms）。全部1688個新基因在GO、KEGG、COG、KOG、Pfam、Swiss-Prot、eggNOG或NR 8個數據庫中獲得注釋。

2.4可變剪切分析

從He102與06-247中分別鑒定出3571個和2065個AS事件，分別獨有2421個和993個AS事件，涉及的基因數目分別是1832個和819個，He102獨有的AS基因數目是06-247的兩倍多。這些AS事件包括外顯子互斥、內含子保留、外顯子跳躍、可變5′剪切位點和可變3′剪切位點共5種類型（表4）。其中，內含子保留是最常見的可變剪切類型，分別占He102與06-247可變剪切的60.52%和57.58%；外顯子互斥則是比較罕見的可變剪切類型，在He102與06-247中分別鑒定出10個和4個，占比均不足3%0（圖2）。2.5可變多聚腺苷酸位點分析

在He102與06-247中分別有19636個和18 996個基因具有APA，APA數為1～29個不等，基因數與APA數統計如圖3所示。其中3 286個基因僅在He102中有APA現象，2646個基因僅在06- 247中有APA現象。雖然二者有16 350個相同基因發生了APA數的變化，但仍然有很多基因在兩者中存在APA數和polyA位置的差異。

2.6 IncRNA的鑒定與分析

IncRNA是不編碼蛋白的全長轉錄本，采用CPC、CNCI、CPAT和Pfam 4種方法分別鑒定出1358、554、1209條和453條IncRNA序列，4種方法的交集有453條高質量IncRNA（圖4A），其中He102獨有81條，06-247獨有55條。80%以上的IncRNA屬于long intergemc non-coding RNA（lincRNA），來自基因間區（圖4B）。靶基因預測結果顯示，451個IncRNA具有順式調控靶基因的潛能，97個具有反式調控潛能。

2.7與開花相關的差異基因篩選

隨機抽取不同數量的測序數據單獨進行轉錄本定量分析，得到He102與06-247轉錄本的飽和度模擬檢驗結果，如圖5所示，最終發現有29 062個基因符合定量比較的要求，這些基因的差異比較火山圖和MA圖如圖6A、B所示。在設定閾值條件下篩選到差異表達轉錄本5197個，這些差異表達轉錄本的聚類分析結果如圖6C所示，其中He102_vs_06-247上調的有2453個，下調的有2744個。

差異表達的轉錄本中有GO注釋的4024個，以“vernalization”、“photoperiod”、“gibberellins”和“gibberellic acid”為關鍵詞分別在GO注釋中篩選，累計獲得41個差異轉錄本，它們在兩材料中的差異關系以及NR注釋見表5。其中有DNA甲基化修飾酶、赤霉素受體GIDIA、赤霉素調控蛋白6和7、開花調控因子FRI-like蛋白、FLC異構體、SOCl-like基因AGIA2和AGL72、FT以及BFT等多種與開花相關的基因。

隨機選擇19個差異表達基因（12個下調、7個上調）進行qRT-PCR驗證，檢測結果如表6所示。總計有16個基因（占84.21%）的調控趨勢與測序結果一致，其中10個基因的差異達到顯著水平（Plt;0.05）。雖然有3個基因的驗證結果與測序結果不一致，但在兩個材料之間的差異均未達到顯著性水平（Pgt;0.05）。表明本研究篩選到的差異基因可信度較高，可以為后續研究提供可靠的參考。

3討論

晚抽薹是春白菜育種最重要的目標之一。晚抽薹性是多基因控制的數量性狀，且受溫度、光照等環境因素影響大。近年來，隨著測序技術的進步，基于二代測序技術的比較轉錄組學研究已經用于晚抽薹和早抽薹大白菜材料之間的差異基因篩選，這為更全面闡釋大白菜晚抽薹分子機理提供了重要參考。但由于二代測序的讀段相對較短，通常僅有200～500 bp，而大白菜基因組在長期演化過程中經歷了基因組的三倍化過程，有很多基因存在多拷貝現象，二代測序技術的短讀段無疑不能區分來自同源基因的同源片段，這使得比較轉錄組學產生的差異基因信息存在先天的系統誤差。納米孔測序技術是第三代測序技術，它的最大優點是可以使超長堿基構成的DNA鏈在一個單通道中就能夠被解碼和識別，而不需要將長鏈分割成小的短鏈，因而可以對全長轉錄組進行測序，對于區別同源基因、識別新基因、篩選可變剪切等基因結構變異有更大的優勢。但國內外尚未見采用納米孔測序技術分析晚抽薹和早抽薹大白菜材料的帶葉片花序軸全長轉錄組的相關報道。因此，本研究基于納米孔測序技術，對大白菜晚抽薹材料06-247與早抽薹材料He102春化完成后的帶葉片的花序軸進行了全長轉錄組分析，所得結果與二代測序的結果可以互相印證，互為補充，將對闡釋大白菜晚抽薹分子機理提供重要參考。

本研究通過納米孔測序技術對單個文庫的數據產出分別是82.02 Gb（He102）和62.64 Gb（06-247），是前人二代測序單個文庫數據產出（6 Gb或8.81Gb）的7.1～13.7倍，成為轉錄本結構分析、新轉錄本／基因篩選以及定量比較的基礎。另外，本研究對新轉錄本的篩選質控嚴格，僅保留identity≥0.9、coverage≥0.85且5′端外顯子有差異的序列，以確保獲得高質量全長新轉錄本序列，最終從He102與06-247中分別獲得37309個和30415個高質量全長轉錄本。將兩材料的全長轉錄本序列合并，累計獲得44797個已知基因和1688個新基因，這一結果與Dai等采用二代測序技術獲得的數據（分別是44799個和2280個）相當。進一步進行嚴格的序列比對，發現在44797個已知基因中有22971個是參考轉錄組同源基因的異構體，這些大量存在的異構體基因（即等位基因的多樣性）從分子水平為種質資源多樣性提供了有力證據，而關于大白菜抽薹期同源基因異構體的挖掘尚未見報道。

在比較轉錄組學研究中，為了得到可靠的結果，常設生物學重復。如Dai等在晚抽薹材料JWW和早抽薹材料XBJ不同春化處理時期的葉片的比較轉錄組學研究中，每個處理設3次重復，發現JWW在春化處理5個階段均存在表達差異的基因有2 245個，而XBJ中有3278個，但并沒有比較兩個材料之間的差異表達基因。Wei等以花原基為試材，同樣設3次重復，將晚抽薹DH系“Y410-1”與早抽薹菜芯DH系“CX14-1”以及晚抽薹DH系“SY2004”與早抽薹菜芯DH系“CX14-1”進行比較，分別獲得12932個和12711個差異表達的基因。上述如此龐大的差異基因對后續選擇候選基因以及深入闡釋大白菜耐抽薹分子機制造成一定困難。考慮到一個物種特定時期表達的基因數目是有限的，隨著測序量的增加，檢測到的轉錄本數目會趨于飽和，而且全長轉錄組測序的成本也很高，因此本研究舍棄了生物學重復。但是表達量越高，轉錄本越容易被檢測定量，對于表達量低的轉錄本，則需要更大的數據量才能被準確定量。基于此，本研究采用大的數據產出來彌補不設重復的缺陷，并隨機抽取轉錄本進行了飽和度檢測，僅僅選擇飽和度檢驗合格的轉錄本進行定量比較，以使得定量分析結果更加可靠。經檢驗，最終有29062個基因符合定量比較的要求，在設定的閾值條件下篩選到5197個差異表達基因。這一數據比Dai等的報道結果略多但遠遠低于Wei等的結果。由于兩個大白菜材料的背景差異比較大，這些差異基因未必都與晚抽薹性有關，鑒于春化、光周期、赤霉素等是影響大白菜抽薹性比較重要的因素，我們在差異基因中進一步篩選GO注釋中含有關鍵詞“vernalization”、“photoperiod”、“gibberellins”或“gibberellic acid”的基因，累計得到41個與抽薹直接相關的差異基因，其中有與DNA表觀修飾有關的DNA（cytosine-5）-methyltransferase DRM2，有與赤霉素信號有關的赤霉素受體GIDIA和赤霉素調控蛋白，有調控開花的轉錄因子FRI-like蛋白、FLC異構體、SOCl-like基因AGIA2和AGL72、FT以及BFT等。眾所周知，FLC是春化途徑的重要調控因子，FT、CONSTANS等是光周期途徑的關鍵調控因子，DNA甲基化在春化介導的染色質結構調控方面發揮重要調控作用。因此，這些候選基因可為深入解析大白菜耐抽薹分子機制提供功能基因。qRT-PCR結果表明，將近85%的差異表達基因的表達趨勢與驗證結果相一致，這進一步增加了候選基因的可信度。

另外，本研究還發現了大量He102與06-247獨有的AS、APA及IncRNA，而已有研究證明開花調控因子FLC的表達受IncRNA的調控，因此不排除與春化相關的基因可能也存在AS和APA方面的差異。本研究獲得的大白菜不同耐抽薹性自交系抽薹期轉錄組在結構方面的差異將為揭示大白菜晚抽薹的分子機制提供新的參考。

4結論

大白菜的抽薹性是多基因控制的數量性狀，本研究采用單分子納米孔測序技術對晚抽薹材料06-247與早抽薹材料He102抽薹期花莖組織進行全長轉錄組測序及分析，獲得了1688個新基因，對新的轉錄本進行了數據庫注釋。篩選得到5197個差異表達基因，并在此基礎上進一步以關鍵詞篩選的方法挖掘出了41個與抽薹密切相關的差異基因，部分差異基因經qRT-PCR驗證準確可靠，這為下一步選擇功能基因研究兩個材料抽薹性差異的分子機制提供了信息。同時，我們還發現兩個材料的轉錄本在AS、APA、IncRNA等方面存在廣泛的變異，這些變異可能對大白菜抽薹性具有一定影響。綜上，全長轉錄組測序揭示出大白菜抽薹期轉錄組存在多方面的復雜差異，這些結果可為進一步闡釋大白菜晚抽薹的分子機制提供新的參考。