基于单细胞转录组测序技术筛选蟾蜍耳后腺免疫修复基因

关键词：单细胞转录组测序；中华大蟾蜍；耳后腺；免疫修复基因；实时荧光定量PCR

中图分类号：S865.9：Q781 文献标识号：A 文章编号：1001-4942（2024）11-0133-08

蟾酥是蟾蜍科动物中华大蟾蜍（Bufo bufo gargarizans Cantor）或黑眶蟾蜍（Bufo melanostictusSchneider）耳后腺和皮肤腺体的干燥分泌物，具有抗癌、强心、镇痛、抗炎、抑菌等多种药理作用。蟾酥是多种临床药物，如牙痛一粒丸、六神丸、麝香保心丸等的原材料，每年市场需求量极大，且价格逐年增长。随着野生中华大蟾蜍种群资源的减少，蟾酥药材需求量增大，人工驯养繁殖中华大蟾蜍获取蟾酥药材，成为行业热点。然而，人工饲养蟾蜍取酥后，因耳后腺损伤感染导致蟾蜍死亡率极高，大多蟾蜍在进行一次采酥后就会死亡，这对于蟾蜍资源来说是一种极大的浪费，且蟾酥产量也达不到市场需求。

高通量测序技术已被广泛用于农业研究，包括动植物全基因组测序、基因组重测序、转录组测序等。其中转录组测序技术通过与代谢组学等技术结合已被用于多种植物的机制研究。在转录组测序技术的基础上，单细胞转录组测序（single-cell RNA-sequencing，scRNA-seq）技术是可以在单个细胞水平上对RNA进行高通量测序和分析的一项新技术。该技术可以同时分析成千上万个细胞，能很好地将群体内单个细胞异质性区分开来，显著提高了我们对动物组织、器官和有机体复杂性的理解。单细胞测序技术已应用于多种动物研究，张卫东等在绒山羊皮肤组织中鉴定出17个主要细胞类群，筛选出3个在绒山羊毛乳头细胞中特异性表达的基因SOX2、FGF7与APOD;崔宝山等对牛体内雌性和雄性囊胚进行单细胞测序研究，发现二者的转录组信息存在巨大差异，并指出影响胚胎发育性别的特异性基因为FBP1、GADD45G、FHL2、FOSB昶WNT2B。

本研究对成年中华大蟾蜍耳后腺组织进行单细胞转录组测序，旨在鉴定蟾蜍耳后腺所有细胞，并对耳后腺细胞进行细胞亚群分类，从差异基因着手探索细胞异质性，筛选免疫基因，为提高中华大蟾蜍耳后腺取酥后存活率提供可能，为丰富蟾酥资源提供参考。

1材料与方法

1.1试验材料

本研究所用中华大蟾蜍采自山东康宇中药材有限公司中华大蟾蜍养殖基地，经山东中医药大学张永清教授鉴定为蟾蜍科蟾蜍属的两栖动物中华大蟾蜍（Bufo bufo gargarizans Cantor）。于2021年8月对未取酥和取酥后0、1、2、4、8、12、24、48、72、120、168 h的中华大蟾蜍进行耳后腺取样，液氮速冻后-80℃冰箱保存备用。

1.2 scRNA-seq文库构建

委托GeneDenovo生物技术有限公司（中国广州），利用Illumina测序平台的PE150测序模式进行高通量测序。

1.3测序数据分析

利用10×Genomics官方分析软件Cell Ranger对原始数据进行过滤、比对、定量，鉴定回收细胞，去除低质量的reads后对每个样本的reads数和测序质量进行初步统计，测序结果上传GSA数据库（GSA编号：GRA012164）。

1.4细胞亚群分类

使用Seurat软件进行细胞聚类分析，保留细胞基因数量在330～2300范围内、有效分子条形码（unique molecular identifier，UMls）总数小于5600且单细胞中线粒体基因表达量比例小于10.0%的细胞。利用Harmony算法对数据进行合并及批次效应矫正，对合并后的数据进行主成分分析（principal component analysis，PCA）降维，采用soft k-means clustering算法对降维后的数据进行聚类。基于细胞亚群分类的结果，进一步利用t-分布领域嵌入算法（t-distributed stochastic neighbor embedding，tSNE）对细胞亚群分类结果进行可视化分析。

1.5细胞亚群上调表达基因分析

采用Seurat的秩和检验分别对不同类细胞群进行基因差异表达分析，以目标亚群或对照亚群中基因在25%以上的细胞中有表达、P≤0.01、基因表达倍数lOg2FC≥0.36为指标，筛选上调表达基因。然后利用Gene Ontology（GO）和KEGG数据库对各聚类中基因的功能富集情况进行分析，包括GO功能分类注释、GO功能显著性富集分析及在差异基因中显著性富集的通路。根据GO和KEGG分析结果，筛选与蟾蜍耳后腺取酥后免疫修复相关的基因。

1.6实时荧光定量PCR分析

使用TB Green Premix Ex TaqⅡ和TB Green Fast Primx EX TaqⅡ试剂盒，对成年中华大蟾蜍未取酥和取酥后不同时间段的耳后腺RNA样本进行反转录后进行qRT-PCR分析。选择在蟾蜍生物体中稳定表达的β-actin作为qRT-PCR验证的内参基因，引物采用Primer3软件（ht-tps：//bioinfo.ut.ee/primer3-0.4.01）设计，由Bio-Sune（中国上海）合成。引物序列见表1。

2结果与分析

2.1测序数据分析

利用10×Genomics官方分析软件Cell Ranger对原始数据进行初步分析，共捕获8414个细胞用于文库构建和对端测序，结果如表2所示。共获得389788801个Reads，其中有效条形码占92.3%，有效分子条形码（unique molecular identifi-er，UMls）占99.9%，81.2%的Reads可定位于中华大蟾蜍基因组。滤除细胞基因数小于330和超过2300、UMls超过5600和线粒体基因在单细胞的表达超过10.0%的细胞后，7687个高质量细胞被保留并进行后续分析。

2.2细胞亚群分类

利用Seurat软件的典型相关分析函数对7687个高质量细胞进行聚类分析，通过tSNE非线性聚类的方法得到12个聚类（图1）。每个聚类中的细胞数从34～4053个不等，占细胞总数的0.46%～52.37%（图2）。对获得的12个聚类进行细胞亚群识别，结果（图3）显示，Cluster 0中有红细胞相关基因RHD、Alas2的表达，初步推测为红细胞：Clus-ter 1中COLIA1、SMOC2的表达水平高于其他亚群，推测为来源于结缔组织的成纤维细胞：Cluster 3中类固醇通路相关基因，如HSD1784、AKRID1、cyp2k4、SULT3A1和Cyp2c29存在高表达，推测可能是类固醇分泌细胞：Cluster 4和Cluster 7中存在高表达的内皮细胞marker基因SELE、VWF和CAV1，推测来源于导管的上皮细胞；Cluster 10中T细胞、B细胞抗原受体信号转导和B细胞发育相关基因PTPRC、BLNK显著表达，推测为来源于结缔组织的浆细胞。Cluster 6中也存在基因PTPRC、BLNK的表达，但不如Cluster 10显著，因此，选择通过Clus-ter 10筛选蟾蜍耳后腺免疫修复基因。

2.3细胞亚群上调表达基因分析

为了解各细胞亚群的分子表达特征，采用Seurat的秩和检验对12个亚群中的上调基因进行筛选，并根据GO和KEGG分析其潜在功能。在上调基因的最小细胞百分比大于25%、平均表达量的10g2FC≥0.36的条件下，确定每个聚类中上调表达基因数量（图4），然后选择每聚类中丰度前5的基因绘制热图（图5）。结果显示，Clus-ter 10中的上调基因数为217个，涉及POU2AF1、IRF4、Cd79b和IGHC4等高表达基因。其中，POU2A F1主要在B细胞中表达，调控B细胞的发育和稳态，还是肺气肿的关键调节因子，通过调节B细胞稳态参与慢性阻塞性肺疾病的进展；IRF4是相对分子质量为52 kDa的一种转录因子，其表达受限于免疫细胞，参与B细胞、T细胞、浆细胞、巨噬细胞等多种细胞的分化。

2.4 Cluster 10的GO和KEGG分析

对Cluster 10的上调基因进行GO功能注释与分析，图6展示GO分析的前20个GO terms，均为描述基因参与的生物过程（biological process），主要涉及淋巴细胞活化（lymphocyte ac-tivation）、淋巴细胞分化（lymphocyte differentia-tion）、免疫系统过程（lmmune system process）、白细胞活化（leukocyte activation）和正向调节生物过程（positive regulation of biological process）等，这些生物过程均与免疫系统有关。

对Cluster 10的KEGG富集分析显示，前20个途径包括了代谢（ metabolism）、环境信息处理（environmental information processing）、有机系统（organismal systems）及人类疾病（human diseases）（图7A），其中有机系统及人类疾病为主要途径；KEGG气泡图（图7B）前5个富集途径为B细胞受体信号通路（B cell receptor signaling pathway）、Fc epsilon RI信号通路（Fc epsilon RI signaling pathway）、趋化因子信号通路（chemokine signaling pathway）、原发性免疫缺陷（primary immunodefi-ciency）和Fc γ R介导的吞噬作用（Fc gamma R-mediated phagocytosis）。吞噬作用是最重要的免疫保护机制之一，Fc γ R可以介导刺激信号进入细胞，引起下游信号的改变，并导致肌动蛋白聚合和吞噬。而B细胞可作为免疫效应细胞直接参与体液免疫应答。

上述分析表明Cluster 10可能在机体的免疫应答和免疫功能中发挥重要作用，这与细胞亚群的鉴定结果相一致。

2.5 Cluster 10主要差异基因qRT-PCR分析

对Cluster 10中的差异基因，采用2^-ΔΔCt法计算相对表达量，并采用t-检验分析差异显著性，筛选出MYOIF和CYBB两个基因，对其在采酥后不同时间段的中华大蟾蜍耳后腺组织中的表达进行qRT-PCR分析。结果（图8）显示，两基因相对表达量在取酥后12 h内并未有显著变化，在24～48 h之间表达水平逐步升高，之后缓慢下降，推测两个基因与免疫修复具有正向调节的关系。

3讨论与结论

当动物机体受到损伤或感染时，免疫系统会自发抵抗病原微生物的入侵，并调节机体稳定。蟾蜍皮肤湿润，皮肤分泌物是其抵抗微生物侵染、保护其免受环境和致病性损伤的第一道防线。而当蟾蜍取酥后，耳后腺受到损伤，皮肤分泌功能受到影响，使其防御能力大大降低，极易受到各种病菌的感染导致死亡。蟾蜍耳后腺相对于其他组织在相应甾族化合物的合成代谢上具有明显优势，已在其分泌物中鉴定出多种抗菌肽（Buforins，Peroniins，Brevinins等）。本研究使用10×Genomics系统对蟾蜍耳后腺进行scRNA-seq，鉴定蟾蜍耳后腺中的所有细胞亚群类型，研究细胞异质性，确定免疫细胞亚群，并筛选蟾蜍耳后腺免疫修复相关基因。

在scRNA-seq结果中，Cluster 10中GO富集分析的前20个GO terms均为基因参与的生物过程，其中MYOIF和CYBB参与了大部分生物过程，包括免疫系统过程、白细胞活化、正向调节生物过程、免疫应答、细胞活化、免疫系统过程的调节等。在KEGG分析通过Pathway显著性富集确定差异基因参与的最主要生化代谢途径和信号转导途径，其中MYOIF参与人类疾病中致病性大肠杆菌感染通路，CYBB参与有机系统中的白细胞经内皮细胞迁移通路。

通过qRT-PCR技术对MYOIF和CYBB两个基因在蟾蜍耳后腺采酥后不同时间的相对表达量进行验证，发现在采酥12h后其表达量随时间的延长逐渐增加，在48h时表达量最高，推测这两个基因对蟾蜍机体具有正向的免疫调控作用。MYOIF是一种非常规的肌球蛋白，通过控制肌动蛋白动力学在促进巨噬细胞介导的吞噬作用中发挥作用。MYOIF通过将AP2A1/αTAT1乙酰转移酶复合物募集到α-微管蛋白上，使α-微管蛋白乙酰化。乙酰化α-微管蛋白在关键信号分子Syk和CARD9从细胞膜向细胞质的易位中起关键作用，这是抗真菌信号转导的重要步骤。通过降低MYOIF表达，可以减少由IgE交联和Mas相关的G蛋白偶联受体X2（MRGPRX2）刺激诱导的人肥大细胞脱颗粒，且MYOIF对于AKT激活和DRP1磷酸化至关重要，是肥大细胞脱颗粒的重要调节因子。MYOIF取代VAV1的C端SH3结构形成的融合基因VAV1-myolf的表达可诱导CD4+T细胞淋巴瘤，CD4+、CD69+细胞增加，诱导共刺激受体（ICOS）上调，IL-10分泌增加。CYBB编码黄细胞色素b558（也称为gp91phox或NOX2）的β-链，而黄细胞色素b558是吞噬细胞（如粒细胞、单核细胞和巨噬细胞）中烟酰胺腺嘌呤二核苷酸磷酸（NADPH）氧化酶复合物（phox）的必需元素。CYBB在多种免疫细胞中表达，对IFN-γ、LPS和TNF-α等炎症细胞因子产生反应，是吞噬细胞杀微生物氧化酶系统的主要组成部分。通过CRISPR/cas9介导的小鼠CYBB敲除，使转基因小鼠产生杀微生物氧化酶缺陷，对自发微生物感染的敏感性增加，包括葡萄球菌、克雷伯氏菌以及真菌白色念珠菌。

本研究基于scRNA-seq技术，对中华大蟾蜍耳后腺取酥后免疫修复基因进行筛选。通过细胞亚群分析，确定Cluster 10为免疫细胞亚群。通过GO和KEGG对Cluster 10上调表达基因的分析，鉴定出多个参与免疫作用的生物过程与通路，并对筛选出的候选基因进行qRT-PCR验证。根据qRT-PCR分析结果推测MYOIF和CYBB基因为蟾蜍耳后腺免疫修复正向调控基因。本研究可为后续通过生物技术手段对蟾蜍养殖进行人工干预提供数据支持，为挖掘蟾蜍耳后腺免疫修复基因，研究蟾蜍耳后腺取酥后的修复机制，提高蟾蜍采酥后成活率提供理论依据，同时为其他种蟾蜍耳后腺单细胞组学研究提供参考信息。