银杏小孢子母细胞响应秋水仙碱低频多倍化的细胞学和转录组分析

柯 梦 邓厚银 郭文斌 李艳星 米跃骐 李 云 孙宇涵，*

（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 蛋白家族，这些均可作为后续进一步研究中的关键候选调控因子。