CN104946740A

CN104946740A - 一种海洋褐潮生物群落结构的检测方法

Info

Publication number: CN104946740A
Application number: CN201510230272.5A
Authority: CN
Inventors: 宋伦; 刘卫东; 王年斌; 韩家波; 孙明; 宋广军; 吴金浩; 张玉凤; 王昆; 王召会; 马志强; 李楠; 宋永刚; 李爱; 田金; 杨爽; 赵海勃; 刘桂英; 胡超魁
Original assignee: RESEARCH INSTITUTE OF OCEAN FISHERY SCIENCE LIAONING PROVINCE
Current assignee: RESEARCH INSTITUTE OF OCEAN FISHERY SCIENCE LIAONING PROVINCE
Priority date: 2015-05-07
Filing date: 2015-05-07
Publication date: 2015-09-30

Abstract

本发明公开一种海洋褐潮生物群落结构的检测方法，其利用宏基因组学高通量测序技术检测水域环境微微型生物18S rDNA，将得到的OTU和物种注释结合，再对OTUs进行多样性、均匀度、丰度、优势度指数分析，从而判断出海洋褐潮生物的优势程度；并有利于统计判断出褐潮暴发优势种优势度演变规律，本发明实施例设计并验证了针对褐潮生物设计的引物序列，通过对实际发生的褐潮暴发案例的验证，证明本发明所述方法可以鉴定海水中海洋褐潮的生物群落结构，为褐潮监测业务化提供技术支持，为海水健康养殖规避环境风险提供技术支撑。

Description

一种海洋褐潮生物群落结构的检测方法

技术领域

本发明涉及海洋微微型真核生物种类检测的技术领域，特别涉及一种海洋褐潮生物群落结构的检测方法。

背景技术

褐潮国内外相关研究概况：褐潮是由微微型海藻引起的有害藻华，也可以理解为微微型浮游生物的赤潮即为褐潮，与传统赤潮相比其藻密度极高、耐低光、持续时间较长。美国研究机构将这种有害藻华命名为“褐潮”，以区别一般的“赤潮”，并为国际相关机构所采纳(CosperEM，et al.,1990)。近几年，中国渤海褐潮频发，其主要致灾种类为抑食金球藻(Aureococcusanophagefferens)，属于微微型藻类，该微藻对养殖业的危害极大，会使贝类滞长甚至死亡。中国是继美国和南非之后第三个出现褐潮的国家，对中国来说这是一种新型的海洋生态灾害，2009年至今已连续5年在渤海海域暴发由微微型藻抑食金球藻引发的褐潮，大部分养殖区出现扇贝滞长或者死亡现象。2010年褐潮给河北省造成直接经济损失2.05亿元(中国海洋灾害公报，2010)。2012年褐潮继续向南部扩散，影响到山东沿海地区，2013年7月5日在秦皇岛-绥中附近海域发生的抑食金球藻褐潮持续时间较长，面积达1450平方公里，占渤海赤潮发生面积的77％(2013年北海区海洋环境公报，2013)。

传统对浮游生物鉴定和检测主要是依靠显微镜观察形态学差异进行分类，然而在实际应用中该方法存在着一定的局限性：首先，浮游生物，尤其是微型和微微型浮游生物，其个体微小，不易观察，实际操作困难，常规的镜检无法实现种类鉴定；其次，该方法依赖于经验丰富的分类学研究人员，需要有扎实的分类学功底，且费时费力，不适合大批量样品的研究；最后，许多浮游生物形态学差异并不明显，难以分辨容易造成人为误差。

美国早期对抑食金球藻的检测主要使用免疫荧光法(Anderson DM,et al.,1993)，但是检测周期较长、灵敏度较低。目前主要利用免疫流式细胞技术、定量聚合酶链反应、色素组分分析和分子生物学等技术。

免疫流式细胞技术(Stauffer BA,et al.,2008)利用共轭荧光素异硫氰酸酯(FITC单克隆抗体)单克隆抗体标记抑食金球藻细胞表面进行检测和计数。该方法可以在很大的丰度范围内(10³-10⁶cells ml^-1)对抑食金球藻的样品进行快速准确的检测和计数。但是这种技术在细胞浓度低于5000cells ml^-1时灵敏度较低。

定量聚合酶链反应(Popels LC,et al.,2003)，是基于抑食金球藻的一个独特的18SrDNA序列，开发一种TaqMan分子探针，在检测***中，实时定量PCR来检测较低细胞水平下抑食金球藻的浓度，灵敏度可以达到1cells ml^-1。

针对我国秦皇岛海域发生的褐潮，孔凡洲(Kong FZ，et al.,2012)等对采集到的浮游植物样品分级过滤后，进行样品色素组分分析，结果显示其主要组分有19-丁酰氧基岩藻黄素(Fucoxanthin,Fuco)、硅甲藻黄素(Diadinoxanthin,Diad)和叶绿素a(Chlorophylla,Chla)等，并且But-fuco色素的含量高，再根据藻华细胞数量高、细胞个体小、对贝类摄食具有抑制效应等，确定藻华原因种为抑食金球藻。

宏基因组学和高通量测序技术的出现，加速了微微型生物鉴定技术的步伐。该方法主要是通过从环境样品中提取全部生物的基因组DNA，构建宏基因组文库，研究环境样品中全部生物的遗传组成及其群落功能。具体的是对环境中全部生物的总DNA进行克隆，采用构建宏基因组文库和筛选等手段获得新的生理活性物质；或者根据rDNA数据库进行引物设计，通过***发生学分析得到该环境中生物的遗传多样性和分子生态学信息。高通量测序技术操作简单、成本低廉，能快速高通量得到特定的DNA片段，较之前的分子生物学方法能更全面展示生物群落的组成结构，极大地推动了浮游生物多样性的研究。

然而目前有关褐潮的鉴定方法都是针对抑食金球藻这一种微微型生物，而褐潮不仅仅是这一种，还有很多微微型生物都有可能发生褐潮，上述的定量聚合酶链反应和定量聚合酶链反应方法对于鉴定已知种类是有用的，但检测未知种类的褐潮却有很大的局限性，因为导致褐潮暴发的种类并非一种，及时了解自然水体中各种微微型生物的种类及比例变化规律是预报褐潮的关键，其中褐潮生物群落结构的检测显得尤为重要。

但目前没有针对褐潮生物设计的引物序列，也没有专门检测褐潮生物群落结构的方法。

褐潮其实是从赤潮中分离出来的一个生态学概念，特指微微型真核生物暴发产生的赤潮。传统意义上的赤潮种类鉴定方法一般用显微镜直接观察计数(粒径一般都大于3μm)，根据鉴定种类及其数量即可研究群落结构。由于褐潮种类为小于3μm的微微型生物，利用显微镜鉴定等一般方法无法鉴定其种类，上述所有鉴定方法都是针对抑食金球藻这一种微微型生物，而褐潮不仅仅是这一种，还有很多微微型生物都有可能发生褐潮，而利用DNA测序技术检测目标水域微微型生物18S rDNA，将得到的DNA片段序列输入GenBank进行检索比对，查找对应种类拉丁文及中文名，进而可以鉴定褐潮生物的种类，目前尚未见到针对褐潮生物设计的专有引物序列，也没有将18S rDNA测序比对方法利用在褐潮生物群落结构分析中，本发明是将这种方法进行扩展的同时，专门设计并验证了针对褐潮生物设计的引物序列，利用克隆测序出现的概率统计各种微微型生物的优势度、群落多样性指数、均匀度指数等参数。

2010.11.24公开的发明专利《一种浮游植物群落结构的检测方法》(申请号：201010184810)是利用变性梯度凝胶电泳(DGGE)技术研究浮游植物群落结构的检测方法，但该方法不能鉴定具体种类，只能根据电泳条带统计生物多样性等参数，却无法了解优势种类名称及其优势度等。而本方法可以鉴定出目标水域中褐潮生物的具体种类及优势度，根据具体优势度和多样性指数可以初步判断褐潮暴发的概率，一般当褐潮生物优势度大于0.2，并且多样性指数低于3.0或均匀度指数低于0.8或物种丰度指数低于1.0，基本可以判定即将暴发该种褐潮，而针对具体优势种类可以开展针对性的防治措施。

发明内容

本发明为了解决背景技术中存在的问题，本发明的目的是提供一种新的海洋褐潮生物群落结构的检测方法。该方法首次设计并验证了针对褐潮生物设计的引物序列(序列如SEQ IDNo:1～2)，将18S rDNA高通量测序比对技术运用到褐潮生物的鉴定中，并运用该技术对褐潮生物多样性、优势度、均匀度进行分析，进而了解褐潮生物群落的优势生物，建立了海洋褐潮生物群落结构的检测方法。

本发明的目的是通过如下技术方案来实现的：

一种海洋褐潮生物群落结构的检测方法，该方法包括利用SEQ ID No1～2所述的核苷酸序列分别作为上游引物和下游引物，以海水表层微微型生物样品中分离得到的18S rDNA为模板，进行PCR扩增反应的步骤。

本发明所述的海洋褐潮生物群落结构的检测方法，其优选的技术方案包括如下步骤：

(1)利用SEQ ID No1～2所述的核苷酸序列分别作为上游引物和下游引物，以海水表层微微型生物样品中分离得到的18S rDNA为模板，进行PCR扩增反应；

(2)对于PCR扩增产物，构建DNA文库，对于构建的DNA文库，使用Illumina MiSeq测序平台PE300测序方案进行测序；

(3)测序得到的原始数据经拼接、过滤，得到OTUs(Operational Taxonomic Units)；

(4)将每个OUTs和物种注释结合，对OTUs进行多样性、均匀度、丰度、优势度指数分析，同时对物种注释在各个分类水平上进行海洋褐潮生物群落结构的统计分析。

上述所述的海洋褐潮生物群落结构的检测方法，其优选的技术方案还包括如下步骤：统计同期水样中小型浮游植物种类和数量，根据测序样品中该中小型浮游植物与其他某种微微型真核生物种类间的序列数比例，测算出该种微微型真核生物的数量，参照《海洋赤潮监测技术规程》(HY/T 069-2005)中赤潮判别与生物量分级标准，微微型真核生物数量>10⁷个/升(赤潮生物体长<10μm)时即判为褐潮。其中所述统计同期水样中小型浮游植物种类和数量的方法为：取同期水样，利用显微镜鉴定可鉴定计数的中小型浮游植物优势种和数量。

本发明还公开了本发明所述的检测方法在在判断即将或已经暴发褐潮中的应用。所述应用的具体方法为：根据海洋褐潮生物群落结构的统计分析，当第一优势种单种优势度超过0.20并且多样性指数低于3.0或均匀度指数低于0.8或物种丰度指数低于1.0时，认为即将或已经暴发褐潮。优选的方案中，所述应用的具体方法为：根据海洋褐潮生物群落结构的统计分析，当第一优势种单种优势度超过0.20并且多样性指数低于3.0、均匀度指数低于0.8以及物种丰度指数低于1.0时，认为即将或已经暴发褐潮。

本发明所述OTUs(Operational Taxonomic Units)即为分类单位，在本发明中97％一致性的序列聚类成为OTUs，每个OTUs代表一个物种。

本发明利用宏基因组学高通量测序技术检测水域环境微微型生物18S rDNA，将得到的OTUs(分类单位，Operational Taxonomic Units)和物种注释结合，再对OTUs进行多样性、均匀度、丰度、优势度指数分析，同时对物种注释在各个分类水平上进行群落结构的统计分析，用于褐潮暴发预测。同时镜检同期水样中小型浮游植物种类和数量，根据OTUs比例，分别计算微微型真核生物数量，微微型真核生物数量>10⁷个/升时即判为褐潮。比较采用本发明的检测方法检测实际发生的褐潮前期、中期和后期的海洋褐潮生物群落结构的统计分析结果与镜检结果显示，本发明的检测结果以及由此确定的判断标准与实际发生的此次褐潮的情况完全相符，说明本发明方法具有较好的实用性。

本发明相对于现有技术的有益效果：本发明首次设计并验证了针对褐潮生物设计的引物序列，利用18S rDNA高通量测序比对方法分析褐潮优势种的相对丰度，建立了海洋褐潮生物群落结构的检测方法。本发明的方法可以鉴定出目标水域中褐潮生物的具体种类及优势度，根据具体优势度和多样性指数可以初步判断褐潮暴发的概率，一般当褐潮生物优势度大于0.2，并且多样性指数低于3.0或均匀度指数低于0.8或物种丰度指数低于1.0，基本可以判定即将暴发该种褐潮，而针对具体优势种类可以开展针对性的防治措施。本发明的方法有望应用于对褐潮演变趋势进行预判，为褐潮监测业务化提供技术支持，为海水健康养殖规避环境风险提供技术支撑。

附图说明

图1为多重序列比对结果；

图2为不同样品的Tags和OTUs数目统计；

图3为根据样品间的OTUs聚类分析结果绘制的韦恩图(Venn Graph)；

图4为研究OTUs之间的******，并使用KRONA对物种注释结果进行可视化展示分类级别；

图5为每个样品在各分类水平上的序列构成柱形图；

图6为在门水平上的物种相对丰度柱形图；

图7为物种丰度聚类图；

图8为单样品(LDW1.1样品)中特定物种分类树；

图9为多样品(LDW1.1、1.2、1.3、2.1、2.2、2.3)中特定物种分类树；

图10为用Qiime软件(Version 1.7.0)绘制得到的样品稀释曲线；

图11为用Qiime软件(Version 1.7.0)绘制得到的Chao1指数曲线；

图12为样品的Shannon指数曲线；

图13为Rank Abundance曲线；

图14Beta多样性指数热图；

图15为基于门水平的PCA分析结果；

图16为样品PCoA分析；

图17为样品Bray Curtis距离聚类树。

具体实施方式

下述非限定性实施例可以使本领域的普通技术人员更全面地理解本发明，但不以任何方式限制本发明。以下实施例中的实验方法如无特别说明均为常规方法，所用材料、试剂等均可从生物或化学试剂公司购买。

以下实施例中所用材料、试剂的商品来源为：Tris饱和酚(分析纯，北京百奥森泰生物技术有限公司)；氯仿(分析纯，上海将来实业股份有限公司)，β-巯基乙醇(上海试一化学试剂有限公司)，无水乙醇(分析纯，天津市津东天正精细化学试剂厂)，蛋白酶K(生工生物工程(上海)股份有限公司)，Taq DNA聚合酶和dNTP(TaKaRa)，引物合成由生工生物工程(上海)股份有限公司完成，建库试剂盒(NEBUltra^TMDNA Library Prep Kit forIllumina)，PCR扩增产物的测序以及数据处理采用Illumina MiSeq测序平台PE300测序方案(诺禾致源公司)；本发明涉及的引物以及序列如表1。

表1.本发明涉及的引物及其序列

实施例1

海洋褐潮生物群落结构的检测方法，包括以下a～e的步骤：

a、样品采集与预处理

采集2-5L表层海水，首先经200目的筛绢过滤，然后用0.22μm孔径的硝酸纤维素膜收集微微型、微型及小型浮游生物，将滤膜转移至10mL无菌离心管中，置于-20℃或-80℃冷冻，用于微微型真核生物的分析。

b、DNA的提取

(1)裂解藻细胞：上述过滤滤膜移入1.5mL灭菌离心管中，用灭菌剪刀将其剪碎，加入500μL CTAB(十六烷基三甲基溴化铵法)缓冲液、1μL巯基乙醇、5μL蛋白酶K混和均匀，放入55℃空气浴中1～2h，直到藻液呈透明状；

(3)抽提：加入300μL饱和酚，300μL CI(氯仿∶异戊醇＝24：1)轻柔颠倒5min，使之呈乳浊液；14000×g 4℃离心10min，小心用微量吸头取出上清液(约500μL)至无菌离心管中，加入500μL CI混匀，14000×g 4℃离心10min；

(4)沉淀：小心用微量吸头取出上清液(约400μL)至无菌离心管中，加二倍体积的预冷无水乙醇和1/10体积的NaAc，充分混匀，放入-20℃冰箱30min；14000×g 4℃离心10min，小心弃上清；

(5)洗涤：DNA沉淀用预冷的70％酒精洗两次，每次洗涤后都需离心，小心弃上清；

(6)溶解：将沉淀真空干燥或自然晾干，溶于20-40μLTE缓冲液中，作为DNA模板，置于-20℃待用。

c、PCR扩增

以步骤b制备得到的藻类DNA作为模板，以V4-R/F、V9-R/F、C4-R/F分别作为引物，进行PCR扩增反应，PCR反应体系为50μL，其中包括5U/μL Taq DNA聚合酶0.2μL(Takara)、10×PCR buffer(10mmoL/L Tris-HCl，500mmol/L KCl，2.0mmol/L MgCl₂和0.01％gelatin，pH8.3)4.5μL，dNTPs(200mmol/L)4.5μL，上下游引物10μmol/L各1μL、20ng的模板DNA 2μL，以无菌水补充体系至50μL。

PCR扩增反应程序为：94℃预变性3min；94℃变性30s，58℃退火30s，72℃延伸2min，33个循环；72℃延伸10min。

PCR扩增产物使用1.5％浓度的琼脂糖凝胶进行电泳检测；

d、文库构建和DNA片段测序：

使用建库试剂盒进行文库的构建，构建好的文库经过Qubit定量和文库检测，合格后，使用Illumina MiSeq测序平台进行上机测序；

e、数据分析

(1)测序得到的原始数据，存在一定比例的干扰数据，为了使信息分析的结果更加准确、可靠，首先对原始数据进行拼接、过滤，得到有效数据。然后基于进行OTUs(Operational Taxonomic Units)(97％一致性的序列聚类成为OTUs，每个OTUs代表一个物种)和物种注释结合，从而得到每个样品的OTUs和分类谱系的基本分析结果。再对OTUs进行丰度、多样性指数等分析，同时对物种注释在各个分类水平上进行群落结构的统计分析。

(2)根据得到的每个OTUs比例，计算多样性指数(H′)，进而对褐潮生物群落结构多样性进行分析；所述多样性指数(H′)的计算公式为：

H^{'} = - Σ_{i = 1}^{s} P_{i} \log_{2} P_{i}

其中P_i为第i个OTUs数量占所有OTUs数量的比例。多样性指数是表示群落内种类多样性的程度，用来判断微微型生物群落或生态***的稳定性指标。多样性指数越低，说明微微型生物群落越不稳定，暴发褐潮的风险越高。

(3)均匀度指数(E)的计算公式为：

E＝H′/log₂ ^S

其中H′为多样性指数，S为种类数。均匀度指数表示微微型生物群落中全部物种个体密度的分配状况，均匀度指数越低，说明微微型生物群落越不均匀，暴发褐潮的风险越高。

(4)物种丰度指数(d_Ma)的计算公式为：

d_{Ma} = \frac{(S - 1)}{\log_{2} N}

其中S为种类数，N为总OTUs数量。物种丰度指数表示微微型生物群落中物种多寡程度，丰度指数越低，说明微微型生物群落种类越少，越不稳定，暴发褐潮的风险越高。

(5)物种优势度(Y)的计算公式为：

Y = \frac{n_{i}}{N} f_{i}

其中n_i为第i种的OTUs数，N为总OTUs数。f_i为第i种在各样品中出现的频率。物种优势度表示微微型生物群落中某一物种在其中所占的优势程度，优势度越高，说明暴发该优势种褐潮的风险越高。通过计算上述指数可以为褐潮预判提供重要的预警辅助信息。

(6)数量统计：镜检同期水样中小型浮游植物种类和数量，根据OTUs比例，分别计算微微型真核生物数量，参照《海洋赤潮监测技术规程》(HY/T 069-2005)中赤潮判别与生物量分级标准，微微型真核生物数量>10⁷个/升(赤潮生物体长<10μm)时即判为褐潮。

实施例2

目前针对18S rDNA真核生物分子鉴定技术所使用扩增引物均为真核生物通用引物，为有针对性的检测褐潮生物种类，设计并验证了褐潮藻类专有引物。

一.引物设计

从NCBI数据库下载179条藻类18S rRNA序列，包含11个门，23个纲，179种。应用软件MEGA4对179条序列进行多重序列比对，比对结果见图1。

根据比对结果选取基因区域F(V4)632-653，R(V4)1054-1076，F(V9)，1385-1405，R(V9)1737-1759，其序列如表2。

表2.本发明的藻类18S rRNA可变区及其碱基序列

使用软件DNAMAN检测两引物的互补性，FASTPCR进行评估引物的属性，如表3。

表3.引物属性

为验证设计引物的优越性，选取真核生物通用引物用于比较分析，通用引物序列C4，如表1中F-C4(上游引物)和R-C4(下游引物)

自行设计的V4、V9引物(序列如表1)与已有C4引物对同一样本提取得到的DNA模板进行PCR扩增、测序，对于测序结果采用如下方法进行数据处理。对于样品的采集与处理、样品DNA的提取、PCR扩增、测序方法参见实施例1的步骤a～d。

二.测序数据处理

采用Illumina MiSeq/HiSeq测序平台得到的下机数据(Raw Data)存在一定的低质量数据，会干扰分析的结果，因此在进一步分析前，需要对下机数据进行预处理，具体处理步骤如下：

1)数据拆分：根据Barcode序列将下机数据拆分为不同样品数据，并截去Barcode序列和PCR扩增引物序列；

2)PE Reads拼接：将拆分的数据使用FLASH(V1.2.7，http://ccb.jhu.edu/software/FLASH/)对每个样品的reads进行拼接，得到的拼接序列为原始Tags数据(Raw Tags)；

3)Tags过滤：拼接得到的Raw Tags，需要经过更严格的过滤处理后，得到高质量的Tags数据(Clean Tags)。参照Qiime(V1.7.0，http://qiime.org/index.html)的Tags质量控制流程，进行如下操作：

a)Tags截取：将Raw Tags从连续低质量值(默认质量阈值为<＝3)碱基数达到设定长度(默认长度值为3)的第一个低质量碱基位点截断；

b)Tags长度过滤：Tags经过截取后得到的Tags数据集，进一步过滤掉其中连续高质量碱基长度小于Tags长度75％的Tags；

以上各步骤的数据产出统计结果见表4。

表4.数据预处理统计及质量控制

说明：Raw PE表示原始下机的PE reads数目；Raw Tags是指拼接得到的Tags序列数目；Effective Tags是指Tags过滤低质量和短长度后，并最终用于后续分析的数据；Base是指最终Effective Data的碱基数目；AvgLen指Effective Tags的平均长度；Q20和Q30是指Effective Tags中碱基质量值大于20(测序错误率小于1％)和30(测序错误率小于0.1％)的碱基百分比；GC(％)表示是Effective Tags中的中GC碱基的含量的百分比；Effective(％)表示Effective Tags的数目与Raw data中的PE Reads数目的百分比。LDW1.1、LDW1.2、LDW1.3为同一个环境样品，LDW2.1、LDW2.2、LDW2.3为同一个环境样品，1.1(2.1)、1.2(2.2)、1.3(2.3)所使用扩增引物为自行设计的V4、V9和真核细胞通用引物C4引物。

三.OTU分析及物种注释

为了研究样品的物种的组成多样性信息，用Uparse软件(Uparse v7.0.1001，http://drive5.com/uparse/)对所有样品的全部Effective Tags序列聚类，默认提供以97％和95％的一致性(Identity)将序列聚类成为OTUs(Operational Taxonomic Units)结果。

在OTUs构建过程中，对不同样品的Effective Tags数据，低频数的Tags数据和Tags注释数据等信息进行初步统计，统计结果的展示见图2。图2中，横坐标是样品名；第一纵坐标(Tags Number)是Tags数目；Total Tags(红色)：指过滤后得到的拼接序列总数；Taxon Tags(蓝色)：指用于构建OTUs并且获得分类信息的Tags数目；Unclassified Tags(茶绿色)指用于构建OTUs但没有获得分类信息的Tags数目；Unique Tags(橙色)：指频数为1，并且无法被聚类到OTUs的Tags数目；第二纵坐标(OTUs Number)是OTUs的数目，OTUs(紫色)：指最终得到的OTUs数目。

按所有样品间序列最小值对OTUs聚类结果进行标准化处理后，根据OTUs聚类分析结果，分析不同样品(组)之间的OTUs的共有、特有信息，并绘制韦恩图(Venn Graph)，结果见图3。同一OTUs中的序列被视为是来源于某一个相同分类单元的序列，作为一个假定的分类单元。Uparse构建OTUs时会选取代表性序列(依据其算法原则，筛选的是OTUs中出现频数最高的序列)，将这些代表性序列集合用RDP Classifier(Version 2.2，http://sourceforge.net/projects/rdp-classifier/)与Silva数据库(http://www.arb-silva.de/)进行物种注释分析(设定阈值为0.6～1)，研究OTUs之间的******，并使用KRONA对物种注释结果进行可视化展示(图4)，圆圈从内到外依次代表不同的分类级别，扇形的大小代表不同OTU注释结果的相对比例。

根据物种注释，统计每个样品在各分类水平，即界、门、纲、目、科、属、种(Kingdom,Phylum,Class,Order,Family,Genus,Species)上的序列数目，物种序列构成柱状图见图5。图5中，横坐标(Sample Name)是样品名；纵坐标(Sequence Number)表示注释到该水平的序列数目；柱状图自上而下的灰度顺序对应于右侧的图例灰度顺序。根据物种注释结果，选取在门(Phylum)分类水平上最大相对丰度排名前十的门，生成的物种相对丰度分布柱形图见图6。图6中，横坐标(Sample Name)是样品名；纵坐标(Relative Abundance)表示相对丰度；Others表示最大相对丰度(某个门在所有样品中相对丰度的最大值)最高的10个门之外的所有门的相对丰度的和。V4引物扩增LDW1和LDW2两个样品，柱形图中直观的体现浮游植物(硅藻diatomeae，绿藻Chlorophyta)所占的比例，高于其他两种引物，因此初步判定V4引物优于其他两种。

根据所有样品在属水平的物种注释及丰度信息，选取丰度排名前35的属及其在每个样品中的丰度信息绘制热图，并从分类信息和样品间差异两个层面进行聚类，便于结果展示和信息发现，从而找出研究样品中聚集较多的物种或样品，结果展示见图7。图7中，横坐标为样品信息，纵坐标为物种注释信息，图中左侧的聚类树为物种聚类树；上方的聚类树为样品聚类树；中间热图对应的值为每一行物种相对丰度经过标准化处理后得到的Z值，即一个样品在某个分类上的Z值为样品在该分类上的相对丰度和所有样品在该分类的平均相对丰度的差除以所有样品在该分类上的标准差所得到的值。

基于每个样品的物种分类结果，筛选特别关注的物种(默认是在种水平上，选择每个分类中各样品最大相对丰度前20的分类所对应的结果进行展示)进行物种分类树统计，其中LDW1.1的物种分类树见图8。图8中，不同颜色的圆圈表示不同的分类水平，对应左侧图例；圆圈的大小表示该分类的相对丰度大小；分类名下方的两个数字均表示相对丰度百分率，前者表示该分类占所有分类中的百分率，后者则表示该分类占所选取的分类中的百分率。粗体字的分类表示在该样品中不存在该分类注释，但是存在于其他分析样品中。

多个样品的物种分类树见图9。图9中，圆圈中不同灰度的扇形表示不同的样品，扇形的大小表示样品在该分类上相对丰度的比例大小；分类名下方的数字表示所有样品在该分类上的平均相对丰度百分率，有两个数字，前者表示占所有物种的百分率，后者表示占所选取物种的百分率。

四.样品复杂度分析(Alpha Diversity)

Alpha Diversity用于分析样品内(Within-community)的群落多样性，主要包含三个指标：稀释曲线(Rarefaction Curves)，物种丰富度(Species Richness Estimators)和群落多样性(Community Diversity Indices)。用Qiime(Quantitative Insights Into Microbial Ecology)软件(Version 1.7.0)对样品复杂度指数进行计算并绘制了相应的曲线。横坐标为从某个样品中随机抽取的测序条数；纵坐标为Chao1指数值，用来评估群落物种数；不同的样品使用不同的颜色和线型的曲线表示。

(1)稀释曲线

用Qiime软件(Version 1.7.0)计算并绘制Rarefaction Curve，即稀释曲线。稀释曲线是从样品中随机抽取一定测序量的数据，统计它们所代表物种数目(也即是OTUs数目)，以数据量与物种数来构建的曲线。稀释曲线可直接反映测序数据量的合理性，并间接反映样品中物种的丰富程度，当曲线趋向平坦时，说明测序数据量渐进合理，更多的数据量只会产生少量新的OTUs，见图10。图10中，横坐标为从某个样品中随机抽取的测序条数；纵坐标为基于该测序条数能构建的OTU数量，用来反映测序深度情况；不同的样品使用不同的灰度和线型的曲线表示。

Species Richness Estimators是用来估计群落样品中包含的物种总数。Chao1指数是广泛使用的指数之一，其计算公式：

S_chao1：估计群落样品中包含的物种总数

S_obs：观察到的物种数

n₁：仅有一条序列的OTUs的数量

n₂：仅有两条序列的OTUs的数量

类似于稀释曲线，绘制数据量与Chao1指数的曲线，见图11。图11中，横坐标为从某个样品中随机抽取的测序条数；纵坐标为Chao1指数值，用来评估群落物种数；不同的样品使用不同的灰度和线型的曲线表示。

群落多样性(Community Diversity Indices)包含样品中的物种组成的丰富度(Richness)和均匀度(Evenness)两个评估指标，广泛采用的计算参考是Shannon指数(Shannon's diversityindex)，也叫香农-维纳(Shannon-Wiener)或香农-韦弗(Shannon-Weaver)指数，它的计算考虑到样品中的分类总数(Richness)，和每个分类所占的比例(Abundance)。其计算公式：

H_{shannon} = - Σ_{i = 1}^{S_{obs}} \frac{n_{i}}{N} \ln \frac{n_{i}}{N}

S_obs：观察到的OTUs数量

n_i：某个种的OUT中的个体数

N：群落中的个体数

当样品中只有一个物种时Shannon指数达到最小值0，Shannon值越大说明群落多样性越高，当样品中有两个以上物种且每个物种丰度为1时Shannon指数达到最大。类似于稀释曲线，绘制数据量与Shannon指数的曲线。在Shannon曲线中，当曲线趋向平坦时，说明测序数据量足够大，可以反映样品中绝大多数的微生物信息，见图12。图12中，横坐标为从某个样品中随机抽取的测序条数；纵坐标为Shannon指数，用来估算群落多样性的高低；不同的样品使用不同的颜色和线型的曲线表示。

对不同样品在不同一致性(Identity)阈值水平(默认提供97％，95％两个阈值)下的AlphaDiversity分析的指数进行统计。其中，97％水平下聚类成为一个OTU的序列被认为可能是源自于同一个种的(Species Boundary)的序列；95％水平下聚类成为一个OTU的序列被认为可能是源自于同一个属的(Genus Boundary)的序列。表5为Alpha Indices统计表。

表5 Alpha Indices统计表

将样品中的OTUs按相对丰度(或者包含的序列数目)由大到小排序得到对应的排序编号，再以OTUs的排序编号为横坐标，OTUs中的相对丰度(也可用该等级OTU中序列数的相对百分含量)为纵坐标，将这些点用折线连接，即绘制得到Rank Abundance(丰度等级)曲线，见图13。图13中，横坐标为按OTUs丰度排序的序号(Rank)，纵坐标为对应的OTUs的丰度(Abundance)，不同的样品使用不同的颜色和线型的曲线表示。Rank Abundance曲线可直观的反应样品中包含的分类丰富度和均匀度，即在水平方向，分类的丰度由曲线的宽度来反映，分类的丰富度越高，曲线在横轴上的跨度越大；在垂直方向曲线的平滑程度，反映了样品中分类的均匀程度，曲线越平缓，物种分布越均匀。

五.多样品比较分析(Beta diversity)

Beta Diversity是对不同样品的微生物群落构成进行比较。首先根据所有样品的物种注释结果和OTUs的丰度信息，将相同分类的OTUs信息合并处理得到物种丰度信息表(ProfilingTable)，进一步计算Bray-Curtis距离。最后，通过多变量统计学方法主坐标分析(PCoA，Principal Co-ordinates Analysis)，主成分分析(PCA,Principal Component Analysis)，非加权组平均聚类分析(UPGMA，Unweighted Pair-group Method with Arithmetic Means)等分析，进一步从结果中挖掘各样品间微生物群落结构的差异和不同分类对样品间的贡献差异。

可以用Bray-Curtis距离作为Beta多样性距离作为衡量两个样品间的相异系数的指标，其值越小，表示这两个样品在物种多样性方面存在的差异越小。以Bray-Curtis距离绘制的热图(Heatmap)展示结果见图14。图14中，方格中的数字是样本两两间间的相异系数，相异系数越小的两个样品，物种多样性的差异越小；方格中的数字值代表bray curtis距离。

主成分分析(PCA，Principal Component Analysis)，是一种应用方差分解，对多维数据进行降维，从而提取出数据中最主要的元素和结构的方法。PCA能够提取出最大程度反映样品间差异的两个坐标轴，从而将多维数据的差异反映在二维坐标图上，进而揭示复杂数据背景下的简单规律。如果样品的群落组成越相似，则它们在PCA图中的距离越接近。展示结果见图15。图15中，横坐标表示第一主成分，百分比则表示第一主成分对样品差异的贡献值；纵坐标表示第二主成分，百分比表示第二主成分对样品差异的贡献值；图中的每个点表示一个样品，同一个组的样品使用同一种颜色表示。

对于一系列样品的特征(比如不同的物种注释)组成的多维向量，多维向量里的某些元素本身没有区分性，比如某个物种分类在所有的样品中的丰度都为1，或者近似于1，那么这个特征本身就没有区分性，用它区分不同样品类别，贡献会非常小，很可能无法较好区分样品间的差异。为了挖掘那些变化大的元素，即方差大的那些特征，并除掉那些变化不大的特征，从而使特征留下的都是“精品”，因此我们进行了PCoA分析。

主坐标分析(PCoA，Principal Co-ordinates Analysis)，是一种与PCA类似的降维排序方法，通过一系列的特征值和特征向量排序从多维数据中提取出最主要的元素和结构。区别在于PCA是基于样品的相似系数矩阵来寻找主坐标，而PCoA是基于距离矩阵来寻找主坐标。我们基于Bray curtis距离来进行PCoA分析，并选取贡献率最大的主坐标组合进行作图展示。如果样品的距离越接近，则它们在PCoA图中的距离越接近，结果展示见图16，其中横坐标表示一个主成分，纵坐标表示另一个主成分，百分比表示主成分对样品差异的贡献值；图中的每个点表示一个样品，同一个组的样品使用同一种颜色表示。在PCoA图中样品距离越接近，表示物种组成结构越相似，因此群落结构相似度高的样品倾向于聚在一起，群落差异很大的样品则会远远分开。使用软件：Qiime软件(Version 1.7.0)中的beta_diversity_through_plots.py的PCoA分析模块作图。

为了研究不同样品间的相似性，还可以通过对样品进行聚类分析，构建样品的聚类树。在环境生物学中，UPGMA(Unweighted Pair-group Method with Arithmetic Mean)是一种较为常用的聚类分析方法，它最早便是用来解决分类问题的。UPGMA的基本思想是：首先将距离最小的2个样品聚在一起，并形成一个新的节点(新的样品)，其分支点位于2个样品间距离的1/2处；然后计算新的“样品”与其它样品间的平均距离，再找出其中的最小2个样品进行聚类；如此反复，直到所有的样品都聚到一起，最终得到一个完整的聚类树。

以Bray-Curtis距离矩阵做UPGMA聚类分析，并将聚类结果与各样品在门水平上的物种相对丰度整合展示，展示结果见图17。图17中，左侧是样品相似性聚类树，右侧的是各样品在门水平上的物种相对丰度分布图。

六.鉴定浮游植物种类统计

通过LDW1、LDW2两个样品三种引物测序分析，平均鉴定种数50，分别鉴定浮游植物种类数见表6，可以看出，自主设计的V4引物鉴定的浮游植物种类较其他两种引物更具有优越性。

表6设计引物鉴定的浮游植物种类比较

引物名称	LDW1	LDW2
			V4(自主设计)	60	60
V9(自主设计)	32	55
			C4(通用引物)	38	52

实施例3

2013年6-7月，大连交流岛养殖池塘和近岸海域海水呈金黄色，经实施例1所述的方法对此海域进行褐潮的预警分析，经过对目标海域中取1L表层海水，样品采集与预处理后进行裂解藻细胞、对其进行DNA的提取、按照实施例2中，使用V4-F/R进行PCR、试剂盒建库、采用Illumina MiSeq测序平台PE300测序方案上机测序，分析数据信息，得到OTUs信息以及进一步得到藻类种类和丰度等信息。最后按照实施例1中步骤e的方法进行统计分析，检测结果如表7-9。本实施例中每个检测样本均做3个平行，第一优势种检测结果标准误在±0.04内。

表7.2013年6月20日褐潮前期检测结果

2013年6月20日褐潮前期生物群落统计结果：计算多样性指数为3.902，均匀度指数为0.888，物种丰度指数为1.524。同期利用显微镜鉴定角毛藻(镜检可鉴定计数)丰度为60°10⁴个/L，根据测序种类间的序列数比例，测算出各种类丰度见表7。

表8.2013年7月3日褐潮中期检测结果

2013年7月3日褐潮中期生物群落统计结果：计算多样性指数为2.949，均匀度指数为0.775，物种丰度指数为0.991。同期利用显微镜鉴定角毛藻(镜检可鉴定计数)丰度为300°10⁴个/L，根据测序种类间的序列数比例，测算出各种类丰度见表8。同时对水体中优势种进行镜检计数，发现优势度明显(严重区域达到0.95以上)，水色较重区域主要优势种密度超过10⁷个/L，已经达到赤潮(微微型浮游生物赤潮即为褐潮)警戒的密度(赤潮监测技术规程(HY\T069-2005，赤潮生物体长<10μm时，形成赤潮生物浓度>10⁷个/L)。

表9.2013年7月15日褐潮后期检测结果

2013年7月15日褐潮后期生物群落统计结果：计算多样性指数为3.706，均匀度指数为0.872，物种丰度指数为1.372。同期利用显微镜鉴定刚毛根管藻(镜检可鉴定计数)丰度为10°10⁴个/L，根据测序种类间的序列数比例，测算出各种类丰度见表9。

分析褐潮优势种出现在微型生物群落中的比例，可以对海洋褐潮生物群落结构进行分析，其对褐潮演变趋势进行预判有显著意义，此次褐潮暴发的前中后期的检测结果从表7～9的统计可知：在褐潮暴发前期(2013年6月20日)检测结果显示，第一优势种Aureococcusanophagefferens单种优势度为0.20，褐潮中期(2013年7月3日)第一优势种Aureococcusanophagefferens单种优势度为0.36，褐潮后期(2013年7月15日)第一优势种Aureococcusanophagefferens单种优势度降至0.18。从多样性指数和丰度指数变化趋势也可看出，褐潮前后变化比较明显。褐潮前后多样性指数均在3.7以上，均匀度指数均在0.87以上，物种丰度指数1.3以上，褐潮时多样性指数、均匀度指数、物种丰度指数均下降较大幅度。根据实施例3研究结果，第一优势种单种优势度超过0.20并且多样性指数低于3.0或均匀度指数低于0.8或物种丰度指数低于1.0，可以认为即将或已经暴发褐潮(优势生物浓度>10⁷个/L才能确认为褐潮)。

应用本申请实施例1所述的方法作出的检测结果，与实际发生的此次褐潮的情况完全相符，证明了上述方法不仅可以作为海洋褐潮生物群落结构的分析方法，而且可作为判断褐潮暴发预警的重要依据。

Claims

1.一种海洋褐潮生物群落结构的检测方法，其特征在于，该方法包括利用SEQ ID No1～2所述的核苷酸序列分别作为上游引物和下游引物，以海水表层微微型生物样品中分离得到的18S rDNA为模板，进行PCR扩增反应的步骤。

2.根据权利要求1所述的检测方法，其特征在于，包括如下步骤：

(3)测序得到的原始数据经拼接、过滤，得到OTUs；

3.根据权利要求2所述的检测方法，其特征在于，还包括如下步骤：统计同期水样中小型浮游植物种类和数量，根据测序样品中该中小型浮游植物与其他某种微微型真核生物种类间的序列数比例，测算出该种微微型真核生物的数量，参照《海洋赤潮监测技术规程》HY/T069-2005中赤潮判别与生物量分级标准，微微型真核生物数量>10⁷个/升时即判为褐潮。

4.权利要求1或2所述的检测方法在判断即将或已经暴发褐潮中的应用。

5.根据权利要求4所述的应用，其特征在于，根据海洋褐潮生物群落结构的统计分析结果，当第一优势种单种优势度超过0.20并且多样性指数低于3.0或均匀度指数低于0.8或物种丰度指数低于1.0时，认为即将或已经暴发褐潮。