基于蟾蜍分泌腺转录组对蟾蜍二烯酸内酯合成途径关键基因的挖掘

张灵迂，王普庆，周罗静，高继海，侯飞侠

成都中医药大学药学院 西南特色中药资源国家重点实验室，四川 成都 611137

蟾酥是蟾蜍科动物中华大蟾蜍Bufobufo gargarizansCantor 或黑眶蟾蜍BufomelanostictusSchneider 皮肤腺体的干燥分泌物。蟾蜍皮肤存在2种外分泌腺：一种是黏液腺又称腺泡腺[1-2]，负责产生黏多糖和黏蛋白多糖，黏多糖自发的分泌能使皮肤保持湿润；另一种是浆液腺又称为合胞腺[1]，腺中充满肽类、胍衍生物、生物胺、甾体以及生物碱等内含物[3]，其中甾体类物质蟾蜍二烯酸内酯（bufadienolides，BDs）是蟾蜍抵御天敌的主要毒性物质，也是蟾酥主要的药效成分[4]。BDs 在细胞中作为一种Na+，K+-ATP 酶抑制剂[5]，可增加细胞内钙浓度，破坏细胞内钙稳态，并诱导不同类型的细胞凋亡，是重要的抗癌药物[6-7]。目前BDs 的主要获取途径是从蟾酥中提取，此方法依靠野生蟾蜍为主要原材料，势必会造成野生蟾蜍的衰退甚至灭绝。微生物转化和化学合成是替代从野生蟾蜍提取BDs的重要手段，但目前关于BDs 生物合成途径的研究较少，导致现阶段对其合成的分子机制仍不清楚。

目前BDs 合成途径研究显示，胆固醇是该类化合物的前体物质[8]，涉及的基因家族包括CYP450、UGTs、ATS 等；涉及的关键酶有SRD5β、3β-HSD和CYP11A1 等，研究人员从蟾蜍耳后腺中克隆了3β-羟基类固醇脱氢酶（3βHSD）和类固醇5β-还原酶（5βSRD）[9-10]，并推测这2 种酶是将BDs 母核A/B 环转化为顺式构型的关键酶，其中3βHSD 是类固醇生成组织中孕酮、雌激素和雄激素生物合成的重要步骤，可分别催化孕烯醇酮和脱氢表雄酮转化为孕酮和雄烯二酮；5βSRD 催化孕酮为孕烯醇酮，可以将孕酮还原为5β-二氢孕酮；CYP11A1 在工程酵母试验中被鉴定为切割胆固醇侧链的酶[11]。这些结果推动了BDs 类化合物合成途径的研究进展，但到目前为止，人们对其合成途径中诸多中间体转化所涉及的基因和酶的了解仍处于初级阶段，还需进一步挖掘和研究。

1 材料

中华大蟾蜍样品捕获于课题组自养2 年龄材料，经成都中医药大学侯飞侠副教授鉴定为中华大蟾蜍BufobufogargarizansCantor。

2 方法

2.1 样品采集和前处理

样品经无菌水洗净，乙醚麻醉后剪取皮肤，挑取背部疣粒和耳后腺体，置于磷酸缓冲盐（phosphate buffered saline，PBS）溶液中洗净，放入盛有OCT 包埋剂的包埋盒中，立即放入液氮速冻后置于−80 ℃冰箱中保存备用，实验设置3 个重复。

2.2 蟾蜍腺体切片制备和显微观察

冷冻组织置Thermo 冷冻切片机中切片，切片厚度10～15 μm。采用连续切片法，选取厚薄均匀、完整的切片，用75%～90%梯度酒精洗脱固定后染色。染色方法如下：苏木素染色2～3 min，蒸馏水洗去浮色，分化液3 min，伊红染色3 s，蒸馏水洗去浮色，梯度乙醇（90%～75%）洗脱，风干，光学显微镜下观察。

2.3 激光显微切割捕获分泌腺

冷冻切片在无酶水中洗去包埋剂，浸泡在无水乙醇：冰乙酸（19∶1）溶液中固定2～3 min，使用LeicaDFC7000T 激光捕获显微切割仪器进行目标材料捕获。实验前调试焦距和激光强度等参数，选择合适的绘图工具分别切下耳后腺中浆液腺（CE）、背部疣粒中粘液腺（CN）和浆液腺（CK）3 种分泌腺，每种分泌腺细胞数量不少于5 000 个。使用含细胞裂解液的无酶PCR 管收集捕获的细胞，瞬时离心，−80 ℃冷冻后用于转录组测序。

2.4 转录组建库和测序

采用DNBSEQ Low Input Smart-Seq Eukaryotic mRNA library 方法对捕获的分泌腺细胞进行建库，利用DNBSEQ 和PE100 方法测序。通过SOAPnuke软件清除原始数据里包含接头，未知碱基（N 含量大于5%）以及低质量的reads，获得Clean reads。

2.5 生物信息学分析

2.5.1 转录组数据的质控、比对和定量分析 利用FastQC、Fastp 软件对Clean reads 进行质量评估及剪切、过滤接头和低质量序列，得到高质量clean data。使用STAR-2.6.1b[17]软件将转录本比对在参考基因组中（ASM1485885v1），计算每个样品的比对率。用featureCounts（1.6.0）[18]计算定量获得表达矩阵。用R 软件包DESeq2 根据原始表达矩阵构建dds 矩阵，计算差异表达基因及对样品做主成分分析，使用伪发现率（false discovery rate，FDR）来校正多次检测的P值（FDR≤0.05），并使用表达水平差异的倍数的双对倍数（log2FC）来确定基因表达的显著性差异（log2FC≥2）。

2.5.2 差异基因的功能富集分析 注释序列用蛋白序列（GCF-014858855.1-ASM1485885v1- protein.faa）在eggnog-mapper[eggNOG-mapper (embl.de)]网站中注释，作为GO 富集的背景基因集，使用clusterProfiler（4.4.4）R 包[19]进行GO 富集分析。同时将蛋白序列上传到KAAS [KAAS-KEGG Automatic Annotation Server(genome.jp)] 数据库中，用单向最佳命中率（single-directional best-hit，SBH）方法，以非洲爪蟾Xenopuslaevis、热带爪蟾Xenopustropicalis、南非爪蟾Nanoranaparkeri和人Homosapiens作参考物种，预测功能并建立KEGG 背景基因集，使用clusterProfiler（4.4.4）R 包进行KEGG 富集分析。

太极虎遇到了同样的厄运，长剑总是在他身上不要命的地方浅浅的，一触即收。太极虎不再是太极虎，而是一只病猫。他全身上下像筛子一样冒着鲜血，人也因失血过多而恍忽，长剑在求生的直觉支持下做着毫无目标的胡乱挥动。

2.5.3 蟾酥药效成分合成相关基因的蛋白互作分析和表达量分析 根据差异基因的GO 和KEGG 富集分析结果，结合PubMed 文献库检索挑选出与甾体激素合成、类固醇合成代谢胆汁酸合成等相关功能和通路的基因，并在STRING[STRING: functional protein association networks (string-db.org)]网站中，以非洲爪蟾为参考背景建立蛋白互作网络，并用Cytoscape（3.7.1）[20]绘制网络图。

选取上述分析后存在互作关系的基因，关联到样品表达量矩阵中，使用微生信网站（https://bioinformatics.com.cn/plot_basic_ballon_plot_0 48）的小工具绘制表达量气泡图。

2.6 荧光定量PCR 验证

选取转录组中识别的9 个参与BDs 合成相关差异表达基因，进行荧光定量PCR（qRT-RCR）验证。以测序实验同批次样品的cDNA 为模板，使用GAPDH为内参基因，利用CDS 序列和Primer 5 软件设计qRT-RCR 特异性引物。采用SsoFastTMEva Green Su-peremix（Bio-Rad，USA）进行定量表达检测，采用2−ΔΔCt方法计算基因的相对表达量。

3 结果与分析

3.1 中华大蟾蜍分泌腺显微观察

中华大蟾蜍背部疣粒的分布呈现从脊椎向两侧外延逐渐减少的趋势（图1-a）。耳后腺显微切片观察发现其中的浆液腺形态呈椭圆形，几乎充满真皮层，腺体呈一字型紧密分布，两侧的浆液腺形态较小，中间的较大（图1-b），其外围分布大量的分泌细胞（图1-c）。背部疣粒切片可观察到疣粒中央散在分布有2～4 个较大的圆形浆液腺（图1-d），表皮下方分布2～6 个圆形的粘液腺（图1-e）。

图1 中华大蟾蜍形态及分泌腺显微结构图Fig.1 Morphology and secretory gland microstructure of Bufo bufo gargarizans

3.2 转录组测序与比对结果

实验测得9 个样品共计clean reads 63.01 Gb，平均GC 含量44.45%，平均Q2096.85%，与参考基因组（ASM1485885v1）的平均比对率为74.87%（表1）。背景基因集注释了26 064 条GO terms，24 823 KEGG通路，注释和比对质量较好，可用于后续实验。

表1 中华大蟾蜍分泌腺转录组测序及比对数据Table 1 Transcriptome sequencing and comparison data of secretory gland of Bufo bufo gargarizans

3.3 差异基因表达分析及富集结果

3.3.1 差异基因表达基因分析 主成分分析显示各生物学重复的样本聚集，可得知每组样本之间具有特异性，表明了数据的可靠性（图2-a）。3 种分泌腺中共检测到4 550 个差异基因（DEGs），其中样品CE较CK 上调62 个基因，下调254 个基因（图2-b）；CE 较CN 上调1 062 个基因，下调713 个基因（图2-c）；CK 较CN 上调1 744 个基因，下调712 个基因（图2-d）；CE 和CK 共有上调基因880 个。

图2 样品间差异表达基因数量统计及主成分分析图Fig.2 Number statistics and principal component analysis of differentially expressed genes among samples

3.3.2 差异基因富集分析 将差异基因进行KEGG和GO 富集分析以深入了解差异基因的生物学功能。结果显示，样品CK 和CE 共有316 个DEGs，其中CK 较CE 上调DEGs 254 个（P＜0.05，log2FC≥2），GO 富集主要涉及含吡啶化合物的代谢（GO:0072524）、脂肪酸分解代谢（GO:0009062）、羧酸分解代谢（GO:0046395）和过氧化物酶体（GO:0005777）；KEGG 通路主要富集在甾类激素生物合成（ko00140）、过氧酶体（ko04146）、原代脂肪酸代谢（ko00120）和花生四烯酸途径（ko00590）等途径。CK 较CE 下调62 个DEGs（P＜0.05，log2FC≥2），GO 富集功能主要在含吡啶化合物的代谢过程（GO:0072524）、雌二醇-17β-脱氢酶活性（GO:0004303）、有机酸分解代谢过程（GO:0016054）、羧酸分解代谢过程（GO:0046395）和甲基转移酶活性（GO:0008168）等；在KEGG富集主要涉及缬氨酸、亮氨酸和异亮氨酸的降解（ko00280）、甾类激素生物合成（ko00140）、类固醇生物合成（ko00100）等。

样品CE 和CN 共存在1 775 个DEGs，其中CE 较CN 上调1062 个DEGs（P＜0.05，log2FC≥2），GO 富集主要涉及类固醇的生物合成（ GO:0006694 ）、类固醇代谢过程的调节（GO:0019218）、花生四烯酸的分泌（GO:0050482）和泛素样蛋白连接酶活性（GO:0061659）等；KEGG富集涉及通路有甘氨酸、丝氨酸和苏氨酸代谢（ko00260）、脂肪酸降解（ko00071）和乙醛酸和二羧酸代谢（ko00630）等。CE 较CN 下调713 个DEGs（P＜0.05，log2FC≥2），GO 功能富集主要涉及粘液层（GO:0070701）、蛋白-O-糖基化（GO:0006493）、细胞对凝集素的反应（GO:1990858）和先天免疫反应激活细胞表面受体信号通路（GO:00022200）等；KEGG 富集主要涉及矿物质的吸收（ko04978）、甲状腺激素信号通路（ko04919）、细胞粘结（ko04510）和JAK-STAT 信号通路和轴突导向（ko04360）等。

样品CK 与CN 共有2 457 个DEGs，其中CK较CN 上调1 744 个DEGs（（P＜0.05，log2FC≥2），GO 富集主要涉及在有机酸分解代谢过程（GO:0016054）、类固醇生物合成的过程（GO:0006694）、17-00 甾酮还原酶活性（GO:0072555）和c21 甾体激素生物合成过程（GO:0006700）等；KEGG 富集主要涉及过氧酶体通路（ko04146）、甾类激素生物合成（ko00140）、P40 细胞色素对外源的代谢作用（ko00980）和氨基苯甲酸酯降解（ko00627）等。CK 较CN 下调712 个DEGs（P＜0.05，log2FC≥2），GO 功能富集主要有蛋白质糖基化（GO:0006486）、刺激C 型凝集素受体信号通路（GO:0002223）、外源凝集素响应答（GO:1990858）和体液免疫应答（GO:0006959）等；KEGG 主要富集在ECM-受体交互通路、细胞粘附（ko04510）和中性粒细胞外陷阱形成（ko04613）等通路。

3.4 蟾酥药效成分生物合成基因的筛选

BDs 生物合成途径很可能是复杂的网络结构[20]，目前研究的途径主要有两种，一是通过甲羟戊酸→胆固醇→孕烯醇酮→孕酮→强心甾类[21]，二是胆固醇通过胆汁酸合成途径合成蟾毒内酯类化合物[4]，结合相关文献，在样品CE、CK 上调差异基因的功能富集结果中检索出包括类固醇、睾酮、雌二醇、前列腺素、雄甾酮、花生四烯酸、石胆酸、胆汁酸在内的28 种物质的代谢、合成及运输等功能，发现94 个基因参与调控，见表2。

表2 蟾酥药效成分合成相关功能基因展示Table 2 Display of functional genes related to synthesis of drug constituents of bufonis venenum

3.4.1 相关基因的蛋白互作网络分析 为进一步探究上述基因之间的关系（图3），利用近缘物种非洲爪蟾X.laevis作参考数据集，将上述基因在STRING数据库中计算，结果显示94 个基因中有34 个基因存在互作关系，如图4 所示，主网络富集的KEGG途径是初级胆汁酸生物合成（XLA00120）；缬氨酸、亮氨酸和异亮氨酸的降（XLA00280）、丙氨酸代谢（XLA00640）及脂肪酸生物合成（XLA01040）等。其中的HUB 蛋白有：酰基-CoA 氧化酶1（acyl-CoA oxidase 1, acox1）、甾醇载体蛋白2（sterol carrier protein 2，scp2）、酰基辅酶A 硫代酯酶8（acyl-CoA thioesterase 8，acot8）及羟基类固醇 17-β 脱氢酶 4 型（hydroxysteroid 17-beta dehydrogenase 4，Hsd17b4.L）等。进一步根据STING 网站计算的KEGG 途径，发现基因AMACR、HSD17B4、SCP2和ACOT8参与调控初级胆汁酸合成（图5），基因LOC122926609参与调控石胆酸结合，结合相关文献，推测BDs 合成途径之一是通过胆汁酸初级合成过程，石胆酸可能是该途径中BDs 类化合物的重要前体。

图3 蟾酥药效成分合成相关基因蛋白互作网络图Fig.3 Interaction network diagram of gene proteins related to synthesis of effective component of Bufonis venenum

图4 BDs 合成相关基因与表达量关联分析气泡图Fig.4 Bubble map of BDs synthesis-related genes correlated with expression levels

图5 BDs 生物合成途径推测Fig.5 Predicted map of BDs biosynthetic pathway

3.4.2 相关基因的表达量分析 为探求“2.4”项中的基因在样品CE 组、CK 组和CN 组的表达情况，将上述存在互作关系的基因和参与石胆酸结合的基因LOC122926609均关联到表达量矩阵中。分析结果显示，大部分基因在CE 组和CK 组的表达量大于在CN组；通过与表达量数据关联及前面的蛋白互作分析，可以侧面验证这些基因可能是BDs 下游合成途径（胆汁酸合成途径）的关键基因。据前文描述，发现基因AMACR、HSD17B4、SCP2和ACOT8参与调控初级胆汁酸合成，将该途径调出，并将上述基因表达气泡附上，如图5 所示，推测该途径是BDs 下游合成途径，上述4 个基因是调控该途径的关键基因。

3.5 荧光定量PCR 验证

本研究为验证转录组数据的可靠性，从34 个基因中选取9 个与BDs 合成相关差异表达基因进行qRT-PCR 验证。分析目的基因在9 个组别中的表达差异，与转录组测序结果进行比较，验证实验结果可靠性。结果显示这9 个基因的基因表达量差异与测序结果具有一致性，结果见图6。

图6 差异基因qRT-PCR 验证Fig.6 qRT-PCR verification of differential genes

4 讨论

4.1 激光捕获显微切割技术在中药有效成分分子机制研究中的应用

激光捕获显微切割在肿瘤细胞的分离、病理组织的分离及大脑细胞的分离等方面技术非常成熟[22-23]，通常与其他技术如转录组学、蛋白组学、代谢组学等联合[24-27]，在空间角度分析样品间的功能及成分分析上提供技术支撑[28]。但此技术在中药中应用极少，而中药材恰恰具有不同部位不同功效的特殊性，LCM 技术的应用可把药用部位的研究精准到细胞类群上，推进药用动植物的表观与其基因间的联系。蟾蜍皮肤腺体作为蟾酥分泌合成的重要组织，目前的研究仅停留在组织水平上，造成组织内细胞功能平均化。

有报道解析了蟾蜍3 个部位（耳后腺、背部皮肤、腹部皮肤）的基因表达谱[29]，并阐述了耳后腺及背部皮肤富集了毒素成分相关的功能和通路。本研究从中华大蟾蜍耳后腺和背部皮肤腺体中分离出浆液腺和粘液腺，发现浆液腺整体上调差异基因的富集除基础的代谢外，主要有甾体激素的合成、胆固醇代谢、胆汁酸合成等；另外，在“2.4”项中发现存在互作关系的 34 个基因除个别基因LOC122926609、VEGFA等在CE 中表达量高于CK外，其余大部分基因在CE 和CK 中表达相当，推测背侧浆液腺和耳后腺浆液腺都具有合成BDs 类成分的潜力。

4.2 BDs 类化合物合成途径关键基因推测

BDs 类成分广泛用于各种癌症，在肝癌、胰腺癌和乳腺癌等多种癌细胞治疗具有较强的细胞毒性作用[30-33]，是癌症治疗的潜在药物，该化合物是一类甾体强心苷，在C-17 位与六元不饱和内酯环（α-吡喃酮环）相连[4]。目前对这类成分的研究，主要还是药理及化学合成方向，生物合成的研究较少。但早在1998 年，研究者就对动植物中合成的BDs类成分的情况作出统计，在植物中，主要由景天科、风信子科、鸢尾科、楝树科、毛茛科及檀香科植物合成；动物则主要存在于蟾蜍、萤火虫、颈槽蛇属及某些哺乳动物中[34]，不过后续更加具体的工作进展缓慢。近几年开始兴起对该类化合物的合成途径关键基因的验证，如验证基因酶SRD5β 和3βHSD在孕烯醇酮-孕酮-二氢孕酮转化等[9-10]。除此之外，对BDs 合成路线的推测的研究也有一定的进展，一是参考洋地黄苷（五元内酯环）的合成途径：甲羟戊酸→胆固醇→孕烯醇酮→孕酮→强心甾类，该途径的特点是胆固醇到孕烯醇酮的转化过程中，胆固醇侧链断开，后由丙二酰辅酶A 提供酰基形成五元内酯环[21]；二是胆固醇经胆汁酸初步合成途径合成石胆酸，以石胆酸为中间体经系列反应合成BDs[11]。

在本研究中对BDs 合成路线的推导，更偏向于初级胆汁酸途径。本课题组在CE、CK 样品上调的差异基因的富集分析中发现多个基因参与调控胆汁酸合成，且基因AMACR、HSD17B4、ACOT8及SCP2存在较强的互作关系，并在CE 和CK 中表达趋势一致，且并未发现有基因参与调控孕烯醇酮的转化。除此之外，查阅文献得知，同时给蟾蜍和红海葱接种孕烯醇酮-20-14C，孕烯醇酮是红海葱中BDs 的良好前体，但在蟾蜍中并不能得到此结果，推测孕烯醇酮不是蟾蜍体内合成BDs 的重要中间体[35]；另外有研究从蟾蜍毒液分离出7-α 羟基胆固醇和7-β 羟基胆固醇[35]，这2 种化合物是胆汁酸生物合成的重要中间体，该结果表明蟾蜍BDs 类化合物的合成经由胆固醇-胆汁酸盐途径[36]，由此，推测BDs 生物合成经由胆汁酸初级合成途径合成石胆酸，再经系列反应合成BDS 类成分母核。基因AMACR、HSD17B4、ACOT8和SCP2及LOC122926609该途径的关键基因，作者已开展对上述基因功能的验证，以期能探索 BDs 类成分的生物合成路线。

BDs 是一类重要的抗癌化合物，研究其生物合成途径是开发相关临床药物的重要基础，而其合成途径关键基因的挖掘是推进该途径的重要研究历程。本研究从微观的腺体细胞着手，探究中华大蟾蜍两种不同功能的外分泌腺差异基因的表达，通过生物信息学等手段推测出BDs 下游合成途径的关键基因，旨在探究BDs 生物合成的分子机制。

利益冲突所有作者均声明不存在利益冲突