副溶血弧菌(Vibrio parahaemolyticus)是一种革兰阴性弧菌，广泛存在于海洋生态系统中，是引起人类海产品相关急性肠胃炎的主要病原菌^[1]。副溶血弧菌的毒力因子主要有直接耐热溶血素(thermostable direct hemolysin, TDH)、TDH相关溶血素(TDH-related hemolysin, TRH)、Ⅲ型分泌系统1 (type Ⅲ secretion system 1, T3SS1)、T3SS2、Ⅳ型分泌系统1(type Ⅵ secretion systems, T6SS1)、T6SS2等^[2]。其中，TDH和TRH是副溶血弧菌的主要毒力决定因子，均具有溶血活性、肠毒性和细胞毒性^[3-5]；T3SS1主要具有细胞毒性，而T3SS2主要与肠毒性相关^[6]；T6SS1和T6SS2均具有细胞黏附活性，而T6SS1还具有杀菌活性^[7-8]。此外，一些其他因子如胞外蛋白酶等也和副溶血弧菌的致病性相关^[9]。

副溶血弧菌具有很强的生物膜形成能力^[10]。生物膜是指细菌在固体表面附着生长时，通过胞外多糖(exopolysaccharides, EPS)、脂类、胞外DNA等基质黏连在一起的、具有一定立体结构的菌落集体^[11]。生物膜是细菌适应不利环境的一种生存策略，生物膜中的菌体细胞对抗生素等有害物质的杀伤、温度等生长参数的变化、宿主免疫防御机制等均具有很强的抵抗能力^[11]。生物膜形成需要EPS、Ⅳ型菌毛、鞭毛等结构的参与^[12-14]。在副溶血弧菌中，cpsA-K和scvA-O负责EPS的合成，二者均与生物形成呈正相关^[15-16]。副溶血弧菌表达2种Ⅳ型菌毛，即甘露糖敏感血凝素(mannose-sensitive haemagglutinin Type Ⅳ pili, MSHA)和几丁质调节菌毛(chitin-regulated pili, ChiRP)，分别由VP2692-2707和pilABCD编码，二者均与生物膜形成呈正相关^{[13, 17]}。鞭毛是细菌的“运动器官”，既有助于细菌克服表面张力到达固体表面而促进生物膜形成，又有助于菌体游离生物膜^[12]。副溶血弧菌有单一极生鞭毛和许多侧生鞭毛，分别介导菌体在液体中游动(swimming)和在固体表面群集性爬动(swarming)^[18-19]。此外，外膜蛋白、胞外DNA和膜融合蛋白(membrane fusion protein, Mfp)等也可影响副溶血弧菌生物膜形成^{[14, 20]}。

生物膜形成是一个被多种因素紧密调节的动态过程，如ToxR^[21]、OxyR^[22]、RpoN^[23]、CpsQ^[24]和AphA^[25]等调控子起正调控作用，而OpaR^[26]、QsvR^[14]、RobA^[27]和VP0610^[28]等起负调控作用；此外，环二鸟苷酸(c-di-GMP^[12])、非编码RNA(如Qrr2^[29])和一些环境因素如温度^[30]、盐度^[31]等对生物膜形成也具有调控作用。与副溶血弧菌浮游细胞中的表达水平相比，生物膜细菌中共有956个基因的转录水平发生显著性变化^[32]。然而，有关副溶血弧菌生物膜及其调节网络仍需进一步深入研究。本文利用转录组测序(RNA-seq)技术研究了副溶血弧菌在生物膜形成中间状态下相对于非生物膜形成条件下的基因转录组差异，共发现979个显著性差异表达基因，包括许多推定的调控子基因、毒力基因和生物膜形成相关基因。转录组分析提供了生物膜形成时基因表达的全景图，也为今后进一步研究生物膜形成调控机制提供重要信息。

1 材料与方法 1.1 材料 1.1.1 菌株

副溶血弧菌RIMD2210633于1996年从日本一名旅行者腹泻液中分离而来，属于O3:K6型大流行菌株^[17]。该菌株具有多种毒力基因位点，主要包括T3SS1、Vp-PAI (包含T3SS2和2个拷贝的tdh基因)、T6SS1、T6SS2等，并具有很强的生物膜形成能力^{[10, 17]}。

1.1.2 主要试剂

M肉汤培养基(3.74% Difco marine broth 2216)和HI肉汤培养基(2.5% Bacto heart infusion)购自BD Bioscience；TRIzol Reagent购自Invitrogen公司；2×Taq PCR MasterMix、SuperReal荧光定量预混试剂彩色版(SYBR Green)、FastKing一步法除基因组cDNA第一链合成预混试剂等购自天根生化科技(北京)有限公司；结晶紫购自上海罗恩试剂。

1.1.3 主要仪器

NanoDrop 2000超微量分光光度计购自Thermo Scientific；测序平台Illumina HiSeq来自Illumina；实时荧光定量PCR仪购自Roche。

1.2 培养方法

取10 μL甘油保存的菌种接种于5 mL的HI肉汤中，置37 ℃下200 r/min培养12 h。按1:50 (体积比)稀释转接至5 mL新鲜的HI肉汤中，置于37 ℃下200 r/min培养至OD₆₀₀为1.4 (称为第2轮培养物)。按1:50 (体积比)接种至10 mL的M肉汤中，充分混匀后，分装至24孔细胞培养板中，每孔1 mL，置于30 ℃下150 r/min培养48 h，同时收集生物膜菌体和液体中的游离菌体，供后续的分子生化试验用。

1.3 RNA提取、文库构建与测序

试验组的细菌培养方法如1.2所述。试验组中同时包含生物膜菌体和游离菌体，在培养过程中，生物膜菌体和游离菌体可不断变化、交换，因此本研究将试验组的培养条件(即1.2中所描述的培养方法和条件)定义为生物膜形成中间状态。将第二轮培养物按1:1 000 (体积比)稀释接种至10 mL的M肉汤中，置于37 ℃下200 r/min连续培养(此条件无法形成生物膜)，收集对数中期(OD₆₀₀为1.2)的细菌培养物^[33]，作为对照组。收集试验组和对照组的菌体(每组3个生物学重复)，并溶解在TRIzol中，运送至苏州金唯智生物技术有限公司(GENEWIZ Biotechnology Co. Ltd.)进行后续处理，包括RNA提取、RNA样品质量检测、rRNA去除、mRNA富集、cDNA文库构建、文库纯化、文库定量和文库测序等。

1.4 测序数据质量评估、序列比对

采用高通量测序平台Illumina HiSeq 300PE对文库进行测序，利用软件Bcl2fastq (v2.17.1.14)进行图像碱基识别(base calling)，获得原始测序数据(pass filter fata)。采用软件FastQC (v0.10.1)分析测序数据质量，碱基的质量值(quality scores, Q)以–10log₁₀(e)计算，e为错误率。一般情况下，Q为13 (Q13)的错误率为5%，Q20的错误率为1%，Q30的错误率为0.1%。此外，接头序列或者3'-末端碱基质量过低会对结果产生负面影响，因此需要使用Cutadapt (version 1.9.1)软件对低质量数据进行过滤。采用bowtie2 (v2.2.6)软件将过滤后测序片段(clean data)与参考基因组(副溶血弧菌RIMD2210633)进行比对分析。对完全匹配的reads (total mapped reads)进行统计定位，区分基因(gene)和基因间隔区(intergenic)。对比对到基因组上的total mapped reads进行统计，计算reads数，并取log₂值。

1.5 差异表达基因注释与功能分析

使用Htseq软件(V0.6.1)和FPKM (fragments per kilo bases per million reads)方法计算基因表达量。使用Bioconductor软件包的DESeq2 (V1.6.3)进行基因差异表达分析。当基因表达变化在2倍以上且Q value (fdr, padj)≤0.05时，则认为是差异显著性表达基因(differentially expressed gene, DEG)。以差异基因的FPKM值为表达水平，做层次聚类(hierarchical clustering)分析。利用gene ontology (GO)数据库分析DEG的分子功能(molecular function)、细胞组分(cellular component)及参与的生物过程(biological process)。通过KEGG (Kyoto encyclopedia of genes and genomes)数据库分析DEG参与的最主要的代谢通路(pathway)。利用蛋白质COG数据库(cluster of orthologous groups of proteins)对DEG进行注释和分类分析，可以预测DEG编码产物的功能。

1.6 实时定量PCR (quantitative real-time PCR, qPCR)

提取试验组和对照组中副溶血弧菌的总RNA，利用FastKing一步法除基因组cDNA第一链合成预混试剂盒将其逆转录成cDNA，最后用SuperReal荧光定量预混试剂彩色版试剂盒进行qPCR分析。以16S rRNA基因的表达量为内参，采用经典的2^‒ΔΔCt法对靶基因的转录水平进行相对定量^[34]。所用引物如表 1所示。qPCR试验至少重复2次，每次3个生物学重复，试验结果用平均值±标准差(standard deviation, SD)表示，使用配对Student's t检验分析是否具有统计学差异，以P<0.01表示具有统计学意义。

2 结果与分析 2.1 RNA-seq数据分析

本研究一共测了6个Illumina文库，包括3个对照组(control)和3个试验组(test)，每个文库都含有超过1 400万的原始测序片段(raw reads，表 2)，文库测序结果已保存至National Center for Biotechnology Information (NCBI)数据库(登录号：PRJNA874225)中。利用Cutadapt (version 1.9.1)软件过滤去除污染及接头序列，获得合格序列(clean reads)。Clean reads的Q20和Q30含量分别在98.46%–98.67%和95.32%– 95.89%之间，GC含量在46.87%–47.75%之间(表 2)，接近副溶血弧菌RIMD2210633的平均GC含量(45.4%)^[17]。Clean reads与参考基因组的比对结果显示(表 3)，单一匹配率(uniquely mapped)在97.275%–98.321%之间、多重匹配率(multiple mapped)在0.962%–1.818%之间、总匹配率(total mapped)在98.908%–99.389%之间，这说明过滤后的reads结果可靠，满足后续分析要求。

2.2 差异表达基因分析

使用Bioconductor软件包的DESeq2 (V1.6.3)对转录组数据进行分析，以差异表达变化2倍以上且Q value (fdr, padj)≤0.05为筛选依据，如图 1A所示，一共筛选出979个DEGs，其中379个被下调，600个被上调。以log10 (FPKM+1)值对DEGs进行聚类分析，如图 1B所示，红色表示高表达基因，蓝色表示低表达基因，颜色从蓝到红，表示基因表达水平逐渐增高。聚类分析能将相同功能或密切联系的基因聚集成类，不同的颜色的区域代表不同的聚类分组信息，同组内表达模式相近的基因可能具有相似的功能。

2.3 差异基因的GO富集分析

利用GO数据库对DEGs的功能进行分类统计，如图 2所示，在分子功能(molecular function)方面，富集结果主要集中在铁硫簇结合相关(28个)、结构分子活性(7个)、质子转运ATP合酶活性(5个)、氨基酸结合(3个)、谷氨酸合成酶活性(3个)和mRNA结合(3个)；在细胞组分(cellular component)方面，富集结果主要集中在细胞膜相关(99个)、完整膜组件(73个)、鞭毛相关(5个)及线粒体基质相关(4个)；在生物过程(biological process)方面，富集结果主要集中在跨膜转运(10个)、三羧酸循环(9个)、钠离子转运(8个)、氧化还原反应(7个)、氨基酸(精氨酸、组氨酸、l-谷氨酸盐、谷氨酸盐)合成代谢(21个)、ATP合成耦合质子转运(4个)、碳水化合物转运(4个)、糖醛酸循环(3个)、氧化磷酸化(3个)、糖酵解(3个)、C4-二羧酸运输(3个)、冷应激(3个)、丙酸分解代谢(3个)及四氢嘧啶合成(3个)。

2.4 差异基因KEGG富集

利用KEGG数据库对DEGs参与的pathway进行注释解析，以确定DEGs参与的主要代谢途径和信号通路。富集筛选标准为Q value<0.05，并将富集最显著(Q value值最小)的前30条展示于图 3中。DEGs主要集中在细胞代谢通路上，共有551个基因，其中占比最高的通路是代谢途径(172个)、次级代谢物的生物合成(80个)、不同环境中微生物代谢(63个)、抗生素的生物合成(58个)和氨基酸的生物合成(40个)；其次是ABC转运蛋白，共有52个DEGs；此外，和致病性相关的DEGs有27个；有机系统相关的DEGs最少，只有13个。

2.5 COG注释分析

利用COG数据库对DEGs进行注释和分类分析，以预测DEGs编码的蛋白或蛋白集合的功能，如图 4所示，从数量多少的角度分析，DEGs主要集中在氨基酸转运与代谢、一般功能预测基因、能量产生与转换以及未知功能基因，其余依次为碳水化合物的转运与代谢、无机离子的转运与代谢、转录、信号转导、细胞壁/膜/包膜的生物生成、翻译/核糖体结构及生物生成、翻译后修饰/蛋白质转换/伴侣分子、细胞运动、脂质转运与代谢、胞内运输/分泌/囊泡运输、辅酶转运与代谢、复制/重组/修复、次要代谢物的生物合成/运输/分解、核苷酸的转运与代谢、防御机制以及细胞周期控制/细胞分裂/染色体分化，且在绝大多数的功能分类中，上调基因占大多数。

2.6 关键基因分析

在979个DEGs中，虽然未知功能基因和一般功能预测基因约占1/5，但是也包含许多生物膜形成相关基因、关键毒力基因、推定调控子基因、外膜蛋白基因等。其中，生物膜形成相关基因(表 4)包括10个推定的c-di-GMP代谢基因、1个侧生鞭毛蛋白基因(lafA)、13个极生鞭毛合成相关基因、6个荚膜多糖合成相关基因、6个Scv基因、2个ChiRP合成相关基因、3个MSHA合成相关基因和Mfp基因位点(cpsQ-mfpABC)；关键毒力基因(表 5)包括14个T3SS1基因(包括exsACDE)、14个Vp-PAI基因(1个tdh2和13个T3SS2基因)、3个T6SS1基因和6个T6SS2基因；45个隶属于LysR、GntR、TetR/AcrR、AraC、MarR和MerR等家族调控子的编码基因(表 6)；15个已知或推定的外膜蛋白基因(表 7)。

2.7 qPCR验证

为验证RNA-seq数据的可靠性，本研究从DEGs中选取12个基因(包括8个上调基因和4个下调基因)作为研究靶标，进行qPCR试验，如图 5所示。虽然qPCR结果与RNA-seq数据在变化倍数上有一定的差异，但基因的表达趋势是一致的，即均为上调或者下调，这说明转录组测序结果是可靠的。

3 讨论与结论

目前对副溶血弧菌生物膜的研究多集中在调控基因或环境因素、天然化合物、化学合成物等对生物膜形成的调控上，也有研究通过组学分析生物膜菌体和游离菌体的基因表达谱差异。然而，还未见比较副溶血弧菌在生物膜形成中间状态和非生物膜形成条件下的基因表达差异情况的报道。本研究以非生物膜形成条件为参照，采用RNA-seq分析了副溶血弧菌在生物膜形成中间状态下的基因转录谱，结果发现共有979个基因的转录水平发生了显著性变化，包括379个上调基因和600个下调基因(图 1)，这些基因涉及小分子物质转运与代谢、能量代谢、合成代谢、基因表达、细胞运动、防御等多种细胞通路(图 2‒图 4)。之前的研究发现，相对于游离细菌，在生物膜细菌中共发现了956个显著性差异表达基因，包括537个上调基因和427个下调基因，且多集中在氨基酸转运与代谢基因、一般功能预测基因、未知功能基因、无机离子的转运与代谢基因、细胞运动相关基因等^[32]。虽然2个RNA-seq研究发现的DEGs数量差别不大、涉及到的细胞通路也有重叠，但是上调和下调基因的数量差异巨大，说明培养方法和取样方式影响试验结果。

副溶血弧菌表达单一极生鞭毛和许多侧生鞭毛^{[19, 35]}。极生鞭毛基因位点位于Ⅰ号染色体上，持续表达^[35-36]；侧生鞭毛基因位点位于Ⅱ号染色体上，只在特定的条件下(如接触固体表面时)表达^{[19, 36]}。侧生鞭毛系统只有一个鞭毛蛋白基因，即lafA，而极生鞭毛系统有6个鞭毛蛋白基因，分别为flaA、flaB、flaC、flaD、flaE和flaF^{[19, 35]}。本研究发现lafA和flaA被上调，而flaB、flaD和flaF被下调，此外还有其他9个极生鞭毛相关基因被上调或下调(表 4)，这说明极生鞭毛和侧生鞭毛在生物膜形成过程中的作用以及不同极生鞭毛蛋白基因对生物膜形成的贡献强度可能存在差异。一般认为侧生鞭毛有利于副溶血弧菌成熟生物膜的形成^[12]，但仍然缺乏相应的证据，有待于后续进一步研究。另外，鞭毛基因表达受群体感应系统(quorum sensing, QS)核心调控子AphA和OpaR的调控^[37-39]。低密度时，AphA促进极生鞭毛和侧生鞭毛基因的表达，而高密度时，OpaR起抑制作用^[37-39]。此外，一些其他的调控子，如ToxR^[21]、VpaH^[40]、OxyR^[22]、LafK^[41]、SwrT^[42]、SwrZ^[42]、H-NS^[43]和CalR^[44]等，对swimming、swarming和生物膜形成等均具有调控作用。可见，通过调控极生鞭毛和侧生鞭毛基因的表达可能是调控副溶血弧菌生物膜形成的重要机制之一。

EPS是生物膜基质的主要成分，占生物膜基质干重的90%以上。cpsA-K和scvA-O负责副溶血弧菌的EPS合成，cps或scv突变株的生物膜形成能力均弱于野生株的^[15-16]。本研究发现若干scv基因被显著诱导，而未见显著性表达的cps基因(表 4)。生物膜基质中EPS含量高低与副溶血菌的生物膜形成不相关^[20]，本研究结果似乎支持了这一结论。此外，2个ChiRP基因和3个MSHA基因也被显著性上调(表 4)。MSHA和ChiRP均能促进生物膜形成，但机制不同，前者主要介导菌体黏附到固体表面，而后者主要促进菌体凝集^{[13, 17]}。然而，MSHA似乎并不影响成熟生物膜的形成，因为MSHA基因突变株仍能形成成熟的生物膜^[13]。MSHA也可介导副溶血弧菌黏附到真核细胞表面，提示MSHA可能和致病性有关^[45]。Mfp基因(cpsQ-mfpABC)转录也被显著性上调(表 4)，这和Mfp与生物膜形成能力呈正相关的结论相一致^[14]。另外，共有10个推定的c-di-GMP代谢相关基因的转录水平发生显著性变化，其中7个被上调、3个被下调(表 4)，说明在生物膜形成过程中伴随着c-di-GMP代谢池的改变。菌体细胞内c-di-GMP浓度越高，生物膜形成能力越强。然而，除scrG外^[46]，其余9个基因是否参与c-di-GMP代谢，还有待于进一步确证。本研究还发现了6个显著性下调的荚膜多糖(capsule polysaccharide, CPS)合成相关基因(表 4)，这和CPS与生物膜形成呈负相关的结果一致^[14]。

在关键毒力因子方面，有14个T3SS1基因表达水平具有显著性差异，包括11个下调基因和3个上调基因，下调基因中包含exsACDE (表 5)。ExsA是AraC家族调控子，能直接激活T3SS1基因的转录；ExsD能结合ExsA，阻止ExsA对T3SS1基因启动子的结合；ExsC能结合ExsD，从而阻止ExsD对ExsA的结合^[47]。ExsE能与ExsC相互作用，从而阻止ExsC与ExsD的结合，因而负调控T3SS1基因的表达^[48]。可见，T3SS1的表达受ExsA-ExsC-ExsD-ExsE调控系统的严格调控。许多其他调控子如H-NS^[49]、OpaR^[50]、AphA^[50]、QsvR^[50]等均是通过直接调控exsACDE的转录来调控T3SS1基因的表达。此外，tdh2、13个T3SS2基因以及若干T6SS1和T6SS2基因的转录均被显著性下调(表 5)。T6SS可以调控细菌生物膜形成，比如荧光假单胞菌(Pseudomonas fluorescens)和瓜类果斑病菌(Acidovorax citrulli)的T6SS基因突变株的生物膜形成能力均显著性低于野生株^[51-52]；溶藻弧菌(V. alginolyticus)的磷酸酶PppA (编码基因位于T6SS基因位点内)对c-di-GMP合成和生物膜形成能力均具有抑制作用^[53]。然而，T3SS1、TDH、T3SS2、T6SS1和T6SS2是否影响副溶血弧菌的生物膜形成，还有待于进一步研究。

本研究还发现了45个显著性差异表达的推定调控子基因(表 6)，包含一些已知功能的基因，如swrZ^[42]、phoBR^[54]、exsA^{[47-48, 55]}、aphA^{[25, 37]}、fis^[56]和cpsQ^[14]，但是其余38个调控子基因的功能还是完全未知的，对这些基因功能的深入研究，将有助于解析副溶血弧菌生物膜形成调控网络。DEGs中还包含15个假定外膜蛋白基因，其中5个下调、10个上调(表 7)，表明在生物膜形成过程中重塑了菌体表面的外膜蛋白。然而，这些蛋白是否能在外膜中形成孔蛋白通道，以及它们是否影响副溶血弧菌生物膜形成，还需要进一步研究。

总之，本研究利用RNA-seq技术分析了副溶血弧菌在生物膜形成中间状态和非生物膜生长条件下的基因表达谱差异，发现979个显著性差异表达基因，包括主要毒力基因、关键生物膜形成相关基因、外膜蛋白基因、大量的未知功能基因和一般功能预测基因等。在差异表达基因中，绝大多数的毒力基因位点包括T3SS1、Vp-PAI、T6SS1、T6SS2和CPS等相关基因均是被显著下调的，而一些生物膜形成相关基因位点如Scv多糖基因、Mfp合成基因等是被显著上调的。DEGs中还包含许多假定的调控子基因和若干c-di-GMP代谢基因，表明副溶血弧菌的生物膜形成受到复杂调控网络的紧密调控。此外，DEGs中还有大量与氨基酸代谢、离子代谢、核苷酸代谢以及能量代谢等有关的基因，可见物质和能量代谢与生物膜形成的关系也值得进一步研究。然而，RNA-seq数据仅提供了生物膜形成时基因表达的整体概况，未来还需要更多的研究来解析副溶血弧菌生物膜形成的分子调控机制。