銀杏小孢子母細胞響應秋水仙堿低頻多倍化的細胞學和轉錄組分析

柯 夢 鄧厚銀 郭文斌 李艷星 米躍騏 李 云 孫宇涵，*

（1北京林業大學生物科學與技術學院/國家林業和草原局刺槐工程技術研究中心/林木育種與生態修復國家工程研究中心/林木、花卉遺傳育種教育部重點實驗室，北京 100083；2襄垣縣森防檢疫站，山西 長治 046200）

多倍化是指生物體中多套染色體在同一細胞核內穩定遺傳，產生擁有3 套或以上完整染色體組個體的過程，是植物進化的重要因素［1］。自然界中至少一半的被子植物都經歷過多倍化事件［2］，而在裸子植物中，銀杏的多倍化研究也一直是重點關注內容［3］。通過形成2n 配子參與受精的有性多倍化是自然界多倍體產生的主要途徑［4］，但自然發生多倍體的頻率并不高。人工誘導能極大提高多倍體得率，目前已有銀灰楊（Populus canescens）［5］、杜仲（Eucommia ulmoides）［6］等多個樹種通過誘導2n 配子成功獲得多倍體植株。人工誘導多倍體的方法包括生物學誘導、物理誘導和化學誘導3 類，其中化學誘導類的秋水仙堿處理一直被認為是多倍體育種中最常用、最經濟的手段。

銀杏（Ginkgo biloba L.）為銀杏科（Ginkgoaceae）落葉喬木，是現存裸子植物門中最古老的孑遺植物，我國重要的經濟林樹種，具有廣泛的栽培價值。然而由于市場需求旺盛，目前已獲得的銀杏新品種資源仍然供不應求［7］。研究表明，三倍體在生長速度、果實產量、抗逆性等方面發揮巨大作用［8］，因此培育銀杏三倍體對于提高其綜合利用價值、滿足市場需求至關重要。秋水仙堿處理能夠誘導銀杏小孢子母細胞產生2n 雄配子，但2n 雄配子得率低，最高只能達到7%［9］，不足以用于銀杏多倍體培育。然而，即使已有研究初步推測此現象可能與銀杏小孢子母細胞中的蛋白質有關［10］，但秋水仙堿誘導銀杏2n雄配子低得率的具體原因和分子調控模式仍有待明確。

鑒于此，本研究利用免疫熒光染色技術，結合考馬斯亮藍G-250 法、茚三酮顯色法、四氮唑藍法、愈創木酚法和硫代巴比妥酸法等細胞生理生化檢測技術，從細胞骨架和生理生化角度出發，連續跟蹤觀察和分析秋水仙堿誘導的銀杏小孢子母細胞在各時間點的不同變化，最終經過篩選并確定可能影響銀杏小孢子母細胞受秋水仙堿加倍誘導過程中的關鍵時期。然后通過轉錄組測序技術挖掘響應秋水仙堿誘導的相關基因，以期為解決秋水仙堿誘導銀杏2n 雄配子低得率這一問題提供數據參考，同時為從分子水平探索秋水仙堿誘導生殖細胞加倍過程中的調控機制，從而進一步提高林木有性多倍化個體得率奠定理論基礎。

1 材料與方法

1.1 試驗材料

以北京林業大學鷲峰林場（40.06787N，116.08134E）的一株30 年生銀杏為試驗材料，初春時取帶雄球花的枝條于溫室中水培（20～30 ℃），觀察其生殖細胞發育狀況。當其小孢子母細胞即將進入減數分裂的粗線期時［10］，對雄球花進行0.6%秋水仙堿（colchicine）避光瓶浸法處理2 d（處理組，簡稱為CC）和蒸餾水處理2 d（對照組，簡稱為CK）［11］。

1.2 生理生化變化檢測

收集秋水仙堿和蒸餾水處理0、12、24、36、48、60、72、84、96 h 后的小孢子囊，-80 ℃冷凍保存后用于生化指標測定，每項指標重復3次。其中，采用考馬斯亮藍G-250 法測定可溶性蛋白含量［12］，茚三酮顯色法測定游離脯氨酸（proline，PRO）含量［13］，四氮唑藍法測定超氧化物歧化酶（superoxide dismutase，SOD）活性［13］，硫代巴比妥酸法測定丙二醛（malondialdehyde，MDA）含量［13］，愈創木酚法測定過氧化物酶（peroxidase，POD）活性［13］。

1.3 碳蠟切片

定期隨機摘取處理組和對照組發育良好的3～5個銀杏花芽，用福爾馬林-乙酸-乙醇（formalin-acetoalcohol，FAA）固定液［V（乙醇）∶V（甲醛）∶V（冰醋酸）=18∶1∶1］固定60 min 后，使用聚乙烯二醇（polyethlene glycol，PEG）包埋銀杏樣本，切片觀察［14］。

1.4 免疫熒光觀察

選取秋水仙堿和蒸餾水處理0、12、24、48 和96 h后的銀杏小孢子囊，固定劑處理60 min，用雙抗免疫熒光［一抗：anti-α-tubulin（T-9026，Sigma 公司，美國），二抗：anti-mouse IgG FITC conjugated（F-0257，Sigma公司，美國）］染色，在激光共聚焦掃描顯微鏡下觀察小孢子微管蛋白的變化情況［14］。

1.5 RNA-seq測序及差異表達基因分析

根據免疫熒光觀察和生化檢測結果，選擇不同時間點（0、24、48、96 h）對照組和處理組的材料作為測序樣本，3次生物學重復。

首先利用Sequel 測序儀第三代（PacBio 公司，美國）對樣本全長cDNA 進行測序，將獲得的全長轉錄本序列與已發表的銀杏全基因組測序結果［15］比對，并將三代轉錄組注釋結果和原基因組結果進行合并，作為進行二代轉錄組檢測差異表達基因分析的參考序列比對注釋，然后利用二代高通量測序技術Illumina HiSeq?（Illumina公司，美國）對隨機抽取的處理組和對照組各3 份重復樣本材料進行文庫構建和測序分析，最后經過濾得到clean reads 比對到前期合并后所得參考序列，經過差異表達分析篩選出樣品間差異表達基因，并利用基因本體富集分析（Gene Ontology，GO）和京都基因和基因組百科全書富集分析（Kyoto Encyclopedia of Genes and Genomes，KEGG）對篩選的差異基因進行進一步的功能富集，從而確定其在秋水仙堿加倍誘導下所行使的主要生物學功能。

1.6 數據處理

采用GraphPad Prism 9.0 進行數據統計和生理生化指標變化分析，并應用SPSS 25.0軟件的T檢驗進行統計學檢驗，統計的顯著性水平為P<0.05。

2 結果與分析

2.1 銀杏小孢子母細胞生理生化變化差異

以蒸餾水處理為對照，選擇小孢子母細胞發育的0、0.5、1、3、6、12、24、36、48、60、72、84 和96 h 進行生理生化檢測。結果發現，對照組和處理組的生化指標隨時間增加均有不同程度的變化。與對照組相比，除12 h 以外，48 h 前處理組可溶性蛋白含量始終高于對照組，而48 h后，處理組可溶性蛋白含量呈下降趨勢且低于對照組（圖1-A）。相較于可溶性蛋白，對照組和處理組中游離脯氨酸含量變化趨勢均保持一致，且處理組始終高于對照組，但在36～48 h 對照組含量持續下降時，處理組游離脯氨酸含量開始出現緩慢上升，并在48 h下降（圖1-B）。處理組的超氧化物歧化酶活性在0～12 h 呈下降趨勢，對照組則緩慢上升，延至12 h后，兩種處理方式的超氧化物歧化酶活性均在極速上升（12～24 h）后趨于平緩，且二者間差異較小（圖1-C）。另外，對照組和處理組的過氧化物酶活性在0～48 h 均表現為下降趨勢，但處理間差異較小，48 h 后，處理組的過氧化物酶活性開始緩慢上升且始終高于變化不明顯的對照組（圖1-D）。秋水仙堿誘導的丙二醛含量在0～96 h 均比對照組高，且在6～12 h間，處理組的丙二醛含量整體呈增加趨勢，并在36～48 h后呈穩定的上升趨勢（圖1-E）。上述結果表明，隨時間增加，經秋水仙堿誘導的銀杏小孢子母細胞生理生化指標表現出不同程度的差異，尤其在處理后12、24、48和96 h差異較明顯。

圖1 銀杏小孢子母細胞內不同生理生化指標隨秋水仙堿處理時間的變化情況Fig.1 Changes of different physiological and biochemical indexes in G. biloba microspore mother cells treated with colchicine at different time

2.2 秋水仙堿對銀杏小孢子母細胞微管骨架的影響

2.2.1 銀杏小孢子母細胞微管骨架結構對秋水仙堿的響應 用熒光顯微鏡分別觀察處理組和對照組0、12、24、48 和96 h 的小孢子母細胞。結果發現，與對照組相比，經秋水仙堿處理的小孢子母細胞微管結構在0～12 h變化不明顯（圖2-A、B、F、G），但在24 h，處理組微管結構開始出現減弱現象（圖2-H），48 h時，指示微管骨架的熒光強度減到最弱，甚至消失（圖2-I），而對照組的微管結構在24～48 h 無明顯變化（圖2-C、D）。值得注意的是，至96 h時，對照組中小孢子母細胞正常進行減數分裂，處理組的綠色熒光也在逐漸顯像，但依舊呈彌散狀分布在細胞核內（圖2-E、G）。上述結果表明，48 h 可能是秋水仙堿誘導銀杏小孢子母細胞多倍化形成的關鍵時間點。

圖2 秋水仙堿處理銀杏小孢子母細胞細胞骨架顯微結構Fig.2 Cytoskeletal microstructure of G. biloba microspore mother cells treated with colchicine after different time

2.2.2 銀杏小孢子母細胞微管蛋白分布和微管厚度差異 通過對秋水仙堿處理48 h的銀杏小孢子母細胞微管進行高倍顯微三維掃描發現，對照組銀杏小孢子母細胞內的微管結構呈聚合狀（圖3-A），而處理組小孢子母細胞的微管結構則呈彌散狀（圖3-B）。分析其微管骨架厚度發現，48 h對照組微管結構完整，呈聚合狀，由藍到紅，微管蛋白層次由少到多（圖3-C），處理組則呈彌散狀，微管骨架厚度較薄（圖3-D）。上述結果表明，秋水仙堿處理銀杏小孢子母細胞48 h時，其微管蛋白的分布和厚度均出現了明顯變化，因此可以進一步推測48 h可能是秋水仙堿誘導銀杏小孢子母細胞多倍化的關鍵節點。

圖3 48 h銀杏小孢子母細胞微管高倍顯微掃描結果Fig.3 High magnification microscopy of microtubules after colchicine treatment of G. biloba microspore mother cells for 48 h

2.3 轉錄組測序挖掘與銀杏小孢子母細胞多倍化相關的差異基因

選取秋水仙堿處理組和對照組0、24、48 和96 h 的小孢子母細胞進行轉錄組分析，相關性分析表明樣品間相關性較好（圖4）。轉錄組測序的Q30（質量值≥30）堿基百分比至少為88%，表明測序結果可靠，可用于后續分析。另外，在GO、KEGG 等數據庫中分別注釋到的基因數目為16 807、11 063，注釋率達56.36%。

圖4 樣品間皮爾遜相關性檢驗Fig.4 Pearson correlation among different examples

以矯正后的P值錯誤發現率（false discovery rate，FDR）≤0.05，差異倍數（fold-change，FC）≥2 為篩選條件，鑒定到4 690 個差異表達基因（圖5）。經統計分析發現：與0 h 相比，蒸餾水處理銀杏小孢子母細胞24、48 和96 h 后，差異表達基因的數量變化不明顯；但相比之下，秋水仙堿誘導的差異表達基因個數在24、48和96 h 分別增加了22.1%、83.9%、60.6%，且隨著秋水仙堿處理時間延長，差異表達的基因數目也隨之增多，96 h 篩選到的差異表達基因數目（9 201）比24 h（5 144）多了78.9%（圖5-A）。分析同一時間發現，處理24 h 時，共篩選到了193 個差異表達基因，其中上調基因135 個，下調基因58 個；至48 h 時，鑒定到上調基因684 個，下調基因741 個，共1 425 個差異表達基因；到96 h時，挖掘到的差異基因明顯增多，共3 400個，其中2 165 個上調基因，1 235 個下調基因（圖5-B）。通過Venn 圖分析發現，在24、48 和96 h 存在23 個重疊差異表達基因（圖5-C）。

圖5 差異表達基因統計圖Fig.5 Statistics of number of differentially expressed gene

2.4 差異表達基因功能注釋

將所篩選的差異表達基因注釋到GO 數據庫中，并對不同時間的前20 個GO 富集結果進行統計分析。結果顯示，隨著秋水仙堿處理時間延長，差異表達基因富集的數量明顯增多（圖6）。具體來看，秋水仙堿處理24 h 時，富集到的差異基因主要參與DNA 模板轉錄、細胞內大分子的生物合成、RNA 代謝、RNA 生物合成、含堿基化合物代謝、含堿基化合物生物合成等生物過程（圖6-A）。48 h 時，雖然有少數基因富集到“催化活性”這一分子功能條目上，但是大多還是參與核苷酸結合相關過程、蛋白修飾、ATP結合等生物過程。值得注意的是，48 h富集到與細胞或亞細胞成分運動、微管運動、微管結合等生物過程相關基因，且大分子生物合成過程在前20 個條目中沒有富集（圖6-B）。而96 h時，大部分差異基因富集到細胞內成分、細胞器、細胞質等細胞成分一類上，卻沒有與微管相關的條目，并且細胞大分子生物合成、細胞生物合成、大分子生物合成等生物過程又重新在96 h富集（圖6-C）。

圖6 差異表達基因GO富集前20條目Fig.6 Top 20 terms for GO enrichment analysis of differentially expressed genes

對篩選的差異基因進行KEGG 富集分析，并統計前20個代謝途徑。結果表明，僅“植物-病原菌相互作用”在3 個時間都有富集；植物晝夜節律、酮化物生物合成、類黃酮生物合成、戊糖和葡萄糖醛酸相互轉化在24、48 h 的前20 條通路中存在；另外，24 和96 h 中都注釋到蛋白酶體、GTP 結合蛋白質通路和植物激素信號轉導途徑相關基因；萜類化合物生物合成、谷胱甘肽代謝、轉運蛋白通路相關基因在48和96 h前20條通路中均有富集。而內質網蛋白質加工相關的基因只在48 h時富集（圖7）。

圖7 差異表達基因KEGG富集前20通路Fig.7 Top 20 pathways of KEGG enrichment analysis of differentially expressed genes

2.5 銀杏小孢子母細胞中響應秋水仙堿誘導的轉錄因子

對各個時間段的差異表達基因進行統計分析，發現在24、48 和96 h 分別鑒定到了12、41、54 個，共21 類轉錄因子，其中ERF、DOF、bHLH 和WRKY 轉錄因子家族廣泛存在于各個時間點（圖8）。具體來看，在秋水仙堿處理24 h 后，發現有4 類轉錄因子家族共12 個基因，包括ERF（9 個）、DOF（1 個）、bHLH（1 個）以及WRKY（1 個）參與表達調控；而在48 h，除ERF（5 個）、bHLH（9 個）、DOF（1 個）、WRKY（3 個）外，MYB、AP2、HD-ZIP、NAC、B3、HSF、NF-YC、RAV、Trihelix、WOX、YABBY 等13 類轉錄因子也高差異表達；在96 h，ERF、bHLH 轉錄因子數目明顯增加，分別為20 和11 個，而MYB、HSF、AP2、HD-ZIP、HSF、Trihelix 等轉錄因子依然參與調控，NF-YB、bZIP、MIKC_MADS、VOZ 等3 類轉錄因子僅高差異表達（圖8-A）。為了準確錨定到小孢子母細胞中響應秋水仙堿誘導的關鍵轉錄因子，用log2（FC）值對48 h鑒定到的41個轉錄因子進行了熱圖分析。結果如圖8-B 所示，處理組中共有23 個轉錄因子上調表達，綜合FPKM（fragments per kilobase of exon model per million mapped fragments）值（≥5），其中log2（FC）值大于2 的轉錄因子有3 個（BBM2、ERF18 和RA212），大于1 的有12 個，分別是AP2/ERF 轉錄因子家族的AIL1、ERF16、ERF79，NAC 轉錄因子家族的NAC25 和NAMB2，MYB-related 中的LHY-2 和ODO1，NF-YC 中的NFYC1，bHLH 中的PAR2，WRKY 家族的WRKY15，Trihelix 家族中的RAV1 以及熱脅迫響應蛋白家族HSF 中的HSFA2e。這些轉錄因子可能為與秋水仙堿誘導銀杏小孢子母細胞相關的重要候選基因。

圖8 各時間點差異表達轉錄因子類別（A）和48 h差異表達轉錄因子熱圖（B）分析Fig.8 Classification of transcription factor families in different treated time（A） and the heatmap of differentially expressed transcription factors（B）

2.6 銀杏小孢子母細胞中持續響應秋水仙堿的差異表達基因

在篩選到的23 個差異表達基因中，僅4 個基因得到功能注釋，分別為熱應激蛋白HSP11、轉錄因子DOF2.4、跨膜運動相關蛋白STC和磷饑餓響應相關蛋白PHL1。其中，HSP11在秋水仙堿處理48和96 h表達量較高，FPKM值分別為110.324 5和98.905 6；DOF2.4在對照組中的表達量始終高于處理組，24 h 兩者相差10 倍，最高96 h 表達量相差59 倍；而STC只在96 h 高表達，且CK96 的表達量（64.147 0 FPKM）高于CC96（22.313 7 FPKM）；PHL1在0 h 的CK 中表達量最高（0.287 2 FPKM），然后隨著秋水仙堿處理時間的延長而逐漸降低（圖9）。另外，兩個未知基因（unknown11和unknown19）在48 h 時高差異表達，可能是該過程中關鍵的候選基因。

圖9 23個重疊差異表達基因熱圖分析Fig.9 Heatmap analysis of 23 overlapped differentially expressed gene

總體來看，篩選到的23個差異表達基因并未全部注釋到Swissport數據庫中，可能由于所參考的基因組信息量較少。

3 討論

秋水仙堿作為抗有絲分裂藥物，其濃度、處理開始時間和持續時間均會影響染色體的加倍［16］。微管是幾乎所有細胞中都存在的細胞器，在形態發生、細胞運動和信號傳導等方面發揮重要作用，微管結構的變化還會影響細胞周期的正常進行。微管的結構亞基是100 kDa 的微管蛋白［17］，微管蛋白通過參與紡錘體、成膜體的形成影響染色體分離和胞質分裂，從而參與真核細胞的細胞周期［18］。微管蛋白的聚合-解聚受溫度、化學藥物等影響［19-21］。用秋水仙堿誘導2n 花粉時，花粉母細胞減數分裂和胞質分離異常與2n 花粉的產生密切相關，微管蛋白則調控減數分裂和胞質分離的正常進行［22］。本研究發現，在處理48 h 時，秋水仙堿顯著抑制了微管蛋白合成，嚴重影響了微管蛋白的聚合-解聚，遲滯減數分裂過程中染色體分離而進一步發育產生2n 花粉，這與李赟等［23］對銀白楊（P.alba）花粉誘導加倍的結果一致。表明秋水仙堿誘導銀杏2n 花粉的原理與其他樹種類似，均是通過影響細胞內的微管蛋白合成，進而干擾染色體的正常分離。

活性氧（reactive oxygene species，ROS）和丙二醛（MDA）是影響生活細胞的重要因子［24-25］，當植物受到脅迫時，細胞內ROS 和MDA 積累，植物細胞會誘導抗氧化酶（POD、SOD等）活性［26］，清除體內自由基。本研究生化指標結果顯示，48 h 時，處理組的總蛋白含量、SOD、POD 和MDA 含量都出現了明顯變化，這與郭鵬［27］的研究結果基本相似，說明48 h 是銀杏小孢子母細胞響應秋水仙堿的一個重要時間節點。

秋水仙堿誘導生殖細胞多倍化困難可能是細胞受脅迫損傷后相關基因表達調控所致［28］。本研究在48 h特異性地富集到了與微管、細胞或亞細胞運動相關下調表達基因的同時，也富集到了上調表達的內質網蛋白質加工相關基因，這表明在48 h時，秋水仙堿與微管蛋白二聚體結合可能抑制了微管的組裝，但未減少微管的解聚，導致微管減少［29］，進而減少了花粉細胞內染色體的聚集，降低2n花粉得率。

進一步尋找影響上述過程的關鍵因子至關重要。本研究在秋水仙堿誘導48 h 的花粉母細胞中鑒定到了41 個差異表達的轉錄因子，其中最大的三個家族是AP2/ERF（7）、MYB-related（2）和NAC（2）。Aya 等［30］研究發現，編碼AP2 類型的轉錄因子SMOS1 可引起細胞內微管取向異常。AP2 還能參與微管乙酰化，影響細胞定向運動和趨化性［31］。此外，細胞骨架解體是花粉絨氈層凋亡的重要表現［32］。研究發現，MYB80轉錄因子可以通過調節細胞程序性死亡時間調控絨氈層發育［33］。而細胞周期的研究結果表明，MYB-related 轉錄因子家族成員Cef1p蛋白功能的缺失會導致α-微管蛋白合成基因TUB1的mRNA 低效剪接，導致細胞周期停滯在G2/M 期［34］。這說明本研究篩選到的AP2/ERF、MYB-related 家族的轉錄因子可能在秋水仙堿誘導銀杏小孢子母細胞的過程中發揮著關鍵作用。

本研究中NAC 轉錄因子和熱應激蛋白HSPs 在秋水仙堿處理過程中均出現特異性表達。NAC轉錄因子作為植物中最大的轉錄因子家族之一［35］，已被證明可能是花粉減數分裂和四分體形成的基因調控模塊中的關鍵基因［36］。研究發現，在水稻花粉中，Osj10gBTF3（一種NAC 蛋白）可以通過與Hsp90 家族的分子伴侶OsHSP82協調作用，影響下游蛋白質的折疊和組裝，從而影響花粉發育［37］。而以往研究也發現，在性器官中高表達的Hsf11，在脅迫條件下會負調控應激轉錄因子HsfA2e，從而調節植物對脅迫的響應［38］。另外，HsfA2和Hsp17-CII 在番茄花藥發育過程中受到精細調控，并在熱應激條件下得到進一步誘導［39］，同時，HsfA2的抑制也會降低花粉在減數分裂階段受到脅迫時的生活力［40］。本研究通過對比分析，在48 h 的處理組中篩選到了高差異表達的熱應激響應基因HSP11、熱應激轉錄因子HSFA2e以及NAC轉錄因子家族中的NAC25和NAMB2，綜合已有研究，推測這些基因可能在銀杏小孢子母細胞響應秋水仙堿的基因調控模塊中發揮著重要作用，但這一假設仍有待進一步研究證實。

4 結論

本研究通過避光瓶浸法，用0.6%的秋水仙堿處理進入減數分裂粗線期的銀杏小孢子母細胞，在處理的0～96 h內對其微管骨架進行顯微觀察，并分析小孢子囊內可溶性蛋白、游離脯氨酸、丙二醛含量以及超氧化物酶和過氧化物酶活性等生理生化指標變化趨勢，發現48 h是秋水仙堿誘導銀杏小孢子母細胞2n雄配子的關鍵時間節點，并利用轉錄組測序技術篩選到了4 690個差異表達基因，其中熱應激響應基因HSP11在48 h 處理組中表達量增加，另外通過對比分析，在48 h篩選到了23 個上調表達的轉錄因子，其中最大的三個家族是AP2/ERF、MYB-related 和NAC 蛋白家族，這些均可作為后續進一步研究中的關鍵候選調控因子。