工作空间

文章信息

徐军, 王瑞, 李锐, 高海侠, 张瑞良, 王文静, 赵霞, 李琳
鼠伤寒沙门菌baeSRacrB双基因缺失株相关耐药基因的RNA-seq分析
生物技术通报, 2019, 35(12): 76-84

XU Jun, WANG Rui, LI Rui, GAO Hai-xia, ZHANG Rui-liang, WANG Wen-jing, ZHAO Xia, LI Lin
RNA-seq Analysis of Drug-resistant Genes Associated with baeSR and acrB Double Gene Deletion Strains of Salmonella typhimurium
Biotechnology Bulletin, 2019, 35(12): 76-84

文章历史

收稿日期:2019-05-27

鼠伤寒沙门菌baeSRacrB双基因缺失株相关耐药基因的RNA-seq分析
徐军, 王瑞, 李锐, 高海侠, 张瑞良, 王文静, 赵霞, 李琳     
安徽农业大学动物科技学院,合肥 230036
摘要:应用转录组测序技术(RNA-seq)筛选鼠伤寒沙门菌baeSRacrB相关的靶基因,为鼠伤寒沙门菌耐药机制的深入研究奠定基础。采用RNA-seq技术分析baeSRacrB基因缺失后表达量变化的相关耐药基因,通过GO富集和KEGG通路探索其主要功能,qRT-PCR验证部分差异基因。以|log2(Fold Change)| > 1且q value < 0.005为条件,其中以CR △acrB和CR为比较组,共筛选到1 320个差异表达基因,894个下调,426个上调;以CR △baeSRacrB和CR为比较组,共筛选到1 377个差异表达基因,972个下调,405个上调。GO富集分类结果显示差异基因主要集中在跨膜运输、细胞黏附、纤毛或鞭毛运动等功能;KEGG数据库显示,差异基因主要富集在代谢途径、ABC转运系统、双组份信号转导系统、鞭毛组装、β-内酰胺抗性等通路。其中,ompRcsgDfliAfliG可通过调节细菌生物膜形成能力调控其对多种抗生素的耐药性;emrAmdtCyejE属于外排泵相关基因,可介导细菌对黏菌素、新生霉素等的耐药性;STM3031可赋予沙门菌头孢曲松抗性;pmrF可参与调节细菌对多粘菌素的抗性。鼠伤寒沙门菌acrBbaeSR基因缺失可影响相关耐药基因的表达。
关键词鼠伤寒沙门菌    baeSR    acrB    RNA-Seq    相关耐药基因    
RNA-seq Analysis of Drug-resistant Genes Associated with baeSR and acrB Double Gene Deletion Strains of Salmonella typhimurium
XU Jun, WANG Rui, LI Rui, GAO Hai-xia, ZHANG Rui-liang, WANG Wen-jing, ZHAO Xia, LI Lin     
College of Animal Science and Technology, Anhui Agriculture University, Hefei 230036
Abstract: The transcriptome sequencing technology(RNA-seq)was used to screen the target genes related to baeSR and acrB of Salmonella typhimurium, which laid a foundation for further study on the drug resistance mechanism of S. typhimurium. RNA-Seq was applied to screen the related drug-resistant genes with varied expressions after the deletion of baeSR and acrB genes, GO enrichment and KEGG pathways to explore their functions, and qRT-PCR to verify some differential genes. Under the condition of |log2(Fold Change)| > 1 and q value < 0.005, a total 1 320 differentially expressed genes were screened by CR △acrB/CR as comparison group, of which 894 were down-regulated and 426 were up-regulated. A total 1 377 differentially expressed genes were screened by CR △baeSRacrB/CR, of which 972 genes were down-regulated and 405 genes were up-regulated. GO enrichment classification results showed that the differential genes were mainly concentrated in transmembrane transport, cell adhesion, cilia or flagellum movement and other functions. KEGG database showed that the differential genes were mainly concentrated in metabolic pathways, ABC transport system, two-component signal transduction system, flagellum assembly, beta-lactam resistance and so on. Among them, ompR, csgD, fliA and fliG regulated their resistances to various antibiotics by regulating the ability of bacterial biofilm formation. emrA, mdtC and yejE belonged to efflux pump-related genes, which mediated the resistance of bacteria to colistin and neomycin. STM3031 endowed Salmonella with ceftriaxone resistance, and pmrF regulated bacterial resistance to polymyxins. In conclusion, the deletion of acrB and baeSR genes in S. typhimurium may affect the expressions of related resistance genes.
Key words: Salmonella typhimurium    baeSR    acrB    RNA-Seq    related drug-resistant genes    

鼠伤寒沙门菌(Salmonella typhimurium)是世界上分布最广泛的食源性病原菌之一,可广泛感染人和动物,引起食源性、感染性腹泻和败血症等,严重威胁养殖业发展和人类健康[1]。目前,抗生素是治疗鼠伤寒沙门菌感染的最有效手段。然而,长期不合理使用和滥用抗生素使得病原菌对抗生素产生耐药性,对兽医临床和动物性食品安全带来极大挑战[2]。多重耐药外排泵系统是引起细菌产生耐药性的最常见原因之一,同时也是细菌生物膜形成和毒力所必需的,因其作用底物广泛且与细菌的多种生理功能相关而被人们广泛关注[3]。其中AcrAB-TolC是革兰阴性菌中最重要的多重耐药外排泵,也是目前研究最清楚的三联外排泵系统,其表达水平受多种调控基因的调节[4]

双组份信号转导系统(Two-component system,TCS)普遍存在于细菌中,可调节包括运动性、毒力、新陈代谢、抗生素耐药性和应激反应等过程[5]。BaeSR双组分系统首次在大肠杆菌中发现,随后相继在沙门菌和鲍曼不动杆菌中发现[6-8]。Lin等[9]证明在沙门菌中BaeSR可直接与STM3031的启动子区域结合,从而调节其对头孢曲松的耐药性。Rosner等[10]发现当某些有毒代谢物在细胞中积累时,BaeSR和CpxAR系统被激活,进而激活AcrD和MdtABC多重耐药外排泵,调控大肠杆菌对多种抗生素的敏感性。本实验室前期发现,acrB的存在可能会掩盖其他外排泵的功能,从而阻碍baeSR对其他多重耐药外排泵的调控[11]。因此,BaeSR和AcrB对于研究抗生素耐药性产生的机制至关重要。

随着测序技术的不断进步,RNA-Seq技术已经成为细菌耐药性研究中的一种常用实验手段。细菌在抗生素作用条件下,本身存在的多种耐药相关基因的表达量会发生一定的变化,利用转录组学研究方法,可以探究细菌耐药前后发生差异表达的基因,从而分析细菌产生耐药性的原因[12]。RNA-Seq技术不仅可以对已知的耐药机制进行确证,还可以从筛选得到的差异表达基因中寻找新的耐药机制,从而揭示其功能及其可能在代谢过程中发挥的作用。本研究通过对鼠伤寒沙门菌CR、CR△acrB及CR△baeSRacrB基因缺失株进行转录组学测序分析,以期筛选与耐药相关差异表达基因并对其表达量变化、发挥的功能进行分析,为进一步探讨鼠伤寒沙门菌产生耐药性的机制,寻找潜在或新的药物靶标,遏制细菌耐药性的产生和传播提供理论依据。

1 材料与方法 1.1 材料

鼠伤寒沙门菌ATCC13311体外诱导环丙沙星耐药株CR由本实验室保存,CR△acrB、CR△baeSRacrB基因缺失株均由本实验室构建。LB肉汤购自上海生工生物有限公司,RNA提取试剂盒、RNA纯化试剂盒均购自QIAGEN公司,TransStart Green qPCR SuperMix购自北京全式金生物技术有限公司。

1.2 方法 1.2.1 样品采集

将鼠伤寒沙门菌CR、CR△acrB和CR△baeSRacrB分别划线接种于LB固体平板,置于37℃培养16 h后,分别挑取单菌落接种于LB液体培养基中,37℃ 180 r/min过夜培养;次日,1:100转接于新鲜LB液体培养基中,生长至对数生长期时收集菌液,于4℃ 5 000 r/min离心10 min后弃上清,并用PBS洗3次,收集菌体,液氮速冻,-80℃保存备用。

1.2.2 RNA的提取与检测

由北京诺禾致源科技股份有限公司对所有样品进行RNA抽提,用琼脂糖凝胶电泳对RNA降解程度以及是否有污染进行检测,用Nanodrop(OD260/OD280)对RNA的纯度进行检测,用Qubit对RNA的浓度进行精确定量检测,最后使用Agilent 2100对RNA的完整性进行精确检测。

1.2.3 文库构建与测序

检测合格后的样品通过Ribo-zero试剂盒去除rRNA,加入fragmentation buffer对样品进行片段化处理以及合成双链cDNA,并用AMPure XP beads对其纯化。对纯化后的双链cDNA进行末端修复、加A尾并连接测序接头,再用AMPure XP beads进行片段大小选择。最后进行PCR扩增,并用AMPure XP beads纯化PCR产物,得到最终的cDNA文库,将构建好的文库使用HiSeq测序平台进行测序。

1.2.4 差异基因表达分析

采用HTSeq软件对各样品进行基因表达水平分析,比对结果使用FPKM(Fragments per kilobase of exon per million fragments mapped)法进行计算,差异表达基因的筛选标准为|log2(Fold Change)| > 1且q value < 0.005。

1.2.5 差异基因GO和KEGG富集分析

采用GOseq方法对所有筛选到的差异表达基因进行GO富集分析,明确差异表达基因在Gene Ontology中的分布状况。以KEGG数据库中Pathway为单位,应用超几何检验,找出在差异表达基因中显著性富集的Pathway,进一步确定差异表达基因所参与的代谢通路及发挥的生物学意义。

1.2.6 差异表达基因的验证

随机选取7个耐药相关差异表达基因,根据NCBI中公布的鼠伤寒沙门菌的基因组序列(CP009102.1),以沙门菌GAPDH为内参基因,利用Oligo 6.0设计所用引物,引物名称及序列见表 1。采用SYBR Green Ⅰ法进行mRNA表达量分析,数据采用2-△△Ct处理。

表 1 引物序列
2 结果 2.1 测序数据质量评估及基因组比对结果

检测样品中CR、CR△acrB和CR△baeSRacrB的质量浓度分别为1 188,372和688 ng/μL;OD260/280分别为2.077、2.11和2.02,OD260/230分别为2.055、1.01和0.76,满足建库测序的要求。从转录组测序数据分析报告中得知,CR、CR△acrB和CR△baeSRacrB分别获得11 693 470、13 342 112和108 24 816条原始序列,对原始数据进行过滤后,剔除低质量数据的Clean reads分别为11 311 842、12 959 780和10 500 966,其错误率分别为0.03%、0.02%和0.02%。通过Bowtie2将过滤后的测序序列进行基因组定位分析,结果显示鼠伤寒沙门菌CR、CR△acrB和CR△baeSRacrB的比对结果分别是95.74%、97.60%和97.19%(表 2)。该结果证明本次测序准确性较高,后续可直接用该部分数据进行分析。

表 2 数据产出质量统计
2.2 基因差异表达分析

采用HTSeq软件对各样品进行基因表达水平分析,以同时满足|log2(Fold Change)| > 1且q value < 0.005为标准对差异表达基因进行筛选。以CR△acrB、CR为比较组,共筛选到1 320个差异表达基因,其中894个基因表达下调,426个基因表达上调(图 1-A);以CR△baeSRacrB、CR为比较组,共筛选到1 377个差异表达基因,其中972个基因表达下调,405个基因表达上调(图 1-B),从上述结果中筛选部分耐药相关差异表达基因(表 3)。

A:CR△acrB/CR差异基因分布;B:△baeSRacrB/CR差异基因分布;红色:上调差异基因;绿色:下调差异基因;蓝色:无显著性差异表达基因 图 1 差异基因分布火山图
表 3 耐药相关差异表达基因
2.3 差异表达基因的GO富集分析

采用GOseq将差异表达基因注释在Gene Onto-logy中,GO分为分子功能(Molecular function)、生物过程(Biological process)和细胞组分功能(Cellular component)3个部分。根据试验目的筛选差异基因后,富集分析差异基因在Gene Ontology中的分布状况。以CR△acrB、CR为比较组得到30个GO功能注释(图 2-A),其中分子功能类4个,无显著富集;生物过程16个,主要集中在跨膜运输、细胞定位、细胞黏附等功能;细胞组分功能10个,主要集中在细胞膜及其相关组成成分。以CR△baeSRacrB、CR为比较组得到30个GO功能注释(图 2-B),其中分子功能类5个,无显著富集;生物过程16个,主要集中在细胞运动、跨膜运输、细胞黏附、纤毛或鞭毛运动、细胞定位等功能;细胞组分功能9个,主要集中在细胞膜及其膜的组成部分。

A:CR△acrB/CR GO富集柱状图;B:CR△baeSRacrB/CR GO富集柱状图;“*”表示显著富集 图 2 GO功能富集
2.4 差异表达基因的KEGG Pathway富集分析

在生物体内,不同基因相互协调行使其生物学功能,通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢和信号转导途径。以CR△acrB、CR为比较组(图 3-A),结果显示1 071个差异表达基因在数据库中注释到88个信号通路,其中581个差异基因表达下调,490个差异基因表达上调;以CR△baeSRacrB、CR为比较组(图 3-B),结果显示1 086个差异表达基因在数据库中注释到86个信号通路,其中605个差异基因表达下调,481个差异基因表达上调。差异基因主要富集在代谢途径、次生代谢产物合成、ABC转运系统、双组分信号转导系统、鞭毛组装、β-内酰胺抗性等通路。

A:CR△acrB/CR差异基因KEGG富集散点图;B:CR△baeSRacrB/CR差异基因KEGG富集散点图;点的大小表示此pathway中差异表达基因个数多少,点的颜色对应于不同的Qvalue范围 图 3 KEGG通路富集
2.5 差异表达基因的qRT-PCR验证

随机筛选7个差异表达基因,以GADPH为内参基因,与RNA-Seq分析结果相比,qRT-PCR的检测结果与其基本一致(图 4)。随机选取的7个差异表达基因在qRT-PCR和RNA-Seq中的诱导表达变化趋势相同,从而证实了RNA-Seq分析结果的可靠性。

图 4 部分差异表达基因qRT-PCR验证结果
3 讨论

近年来,从分子水平上揭示沙门菌的耐药机制,已成为沙门菌最重要的研究方向之一[13]。鼠伤寒沙门菌的耐药性主要依赖于外排泵系统、双组分信号转导系统、生物膜等因素。抗菌药物耐药相关外排泵可分为5大家族,其中RND家族主要以三聚体形式转运底物,通常由染色体基因编码,其结构由膜融合蛋白、RND转运蛋白和外膜因子组成。在革兰阴性菌中,RND家族外排泵作为最主要的多重耐药外排泵通过从细菌细胞中泵出多种抗微生物化合物来赋予其多重耐药性,在沙门菌耐药性方面发挥着不可替代的作用[14]AcrB因其作用底物广泛且在多种革兰阴性菌中均存在其同源物,被认为是对人类健康最重要的RND系统[15]。双组分信号转导机制可以调节细菌对环境的适应性从而影响其对抗生素的敏感性,如BaeSR可通过激活RND型外排泵AcrD和MdtABC的转录改变沙门菌的多要耐药性[10, 16-17]。本实验室致力于鼠伤寒沙门菌耐药机制的研究,通过自杀质粒介导的同源重组系统构建acrBbaeSR的基因缺失株,并对其进行运动性、生物膜的形成、耐药性等方面的研究,结果显示,acrBbaeSR基因缺失导致生物膜形成能力下降,耐药性减弱,表明基因acrBbaeSR可调控鼠伤寒沙门菌的耐药性。为进一步探究其耐药机理,对CR、CR△acrB以及CR△baeSRacrB进行转录组测序。

转录组测序结果显示CR△acrB/CR、CR△baeSRacrB/CR分别筛选到1 320、1 377个差异表达基因。CR△acrB/CR中894个基因表达下调,426个基因表达上调;CR△baeSRacrB/CR中972个基因表达下调,405个基因表达上调。通过对差异表达基因进行GO功能分类发现,差异表达基因主要参与跨膜运输、细胞黏附、纤毛或鞭毛运动以及细胞定位等生物学功能,生物膜以及膜的组成部分等细胞组分功能,这些生物过程均与细菌耐药性有着密切联系。先前得到的acrBbaeSR基因缺失可以导致生物膜形成能力下降、运动性降低、耐药性减弱等试验结果与其基本一致。通过KEGG富集分析和以往文献报道对转录组学结果进行深入分析,筛选出CR△acrB和CR△baeSRacrB共同变化的耐药相关差异表达基因,包括STM3031pmrFfliAOmpRcsgDemrAmdtCyejE等。

STM3031已被证实可赋予沙门菌头孢曲松抗性,本团队对acrB基因缺失后可造成STM3031表达量显著下调,acrBbaeSR基因缺失后表达量进一步下降,与Lin等[9]实验结果一致。KEGG Pathway富集分析显示pmrF属于氨基糖和核苷酸糖代谢途径,可参与调节细菌对多黏菌素的抗性[18]。细菌生物膜的形成与临床感染密切相关,鞭毛在生物膜的形成初期起着重要的介导黏附作用,并为生物膜的形成提供动力[19]。因此,生物膜和鞭毛运动都是引起细菌产生耐药性的重要原因。本研究中OmpR基因表达量显示上调,csgDfliAfliCfliG等相关基因均显示不同程度的下调。其中OmpR属于双组份信号转导系统,可激活外膜孔蛋白相关基因、脂多糖修饰相关基因等的转录[20]。CSG区由csgBAcsgDEFG组成,可编码菌毛等相关基因,csgD是重要的生物膜转录调节因子,在细胞黏附、聚集和生物膜形成方面发挥至关重要的作用[21]。朱春红等[22]通过构建的不同肠炎沙门菌的fliC基因缺失株,证明其在生物膜形成过程中发挥着重要作用。本研究结果证实acrBbaeSR可通过调控沙门菌的鞭毛运动和生物膜形成能力改变其对抗生素的敏感性。

外排泵在细菌多重耐药性方面发挥至关重要的作用,本研究筛选到多个外排泵相关基因均与以往报道相符。其中emrA属于主要易化因子家族(MFS),可与emrBtolC结合形成三联体系;mdtC属于RND家族外排泵,BaeSR可激活mdtC的表达;yejE属于ABC转运系统,均可介导菌株对黏菌素、新生霉素的耐药性[23]。因此,多种外排泵相关基因的表达量发生变化,可能是acrBbaeSR调控鼠伤寒沙门菌产生耐药性的主要原因。在后续研究中,我们将对STM1545fliAcsgDpmrf进行EMSA验证,进一步明确鼠伤寒沙门菌双组分系统BaeSR和AcrB的耐药调控机制。

4 结论

本研究初步探讨了鼠伤寒沙门菌acrBbaeSR基因缺失后的表达变化,发现菌株对抗生素的敏感性可能通过AcrB和BaeSR介导的多重耐药外排泵基因的表达、细胞膜的通透性、核苷酸代谢以及鞭毛运动等多种耐药机制调控。

参考文献
[1]
Cevallos-Almeida M, Martin L, Houdayer C, et al. Experimental infection of pigs by Salmonella Derby, S. Typhimurium and monophasic variant of S. Typhimurium: comparison of colonization and serology[J]. Veterinary Microbiology, 2019, 231: 147-153. DOI:10.1016/j.vetmic.2019.03.003
[2]
Sharma C, Rokana N, Chandra M, et al. Antimicrobial resistance: its surveillance, impact, and alternative management strategies in dairy animals[J]. Frontiers in Veterinary Science, 2018, 4: 237. DOI:10.3389/fvets.2017.00237
[3]
Neuberger A, Du D, Luisi BF. Structure and mechanism of bacterial tripartite efflux pumps[J]. Research in Microbiology, 2018, 169(7-8): 401-413. DOI:10.1016/j.resmic.2018.05.003
[4]
Zwama M, Yamaguchi A. Molecular mechanisms of AcrB-mediated multidrug export[J]. Research in Microbiology, 2018, 169(7-8): 372-383.
[5]
Sandeep T, et al. Two-component signal transduction systems of pathogenic bacteria as targets for antimicrobial therapy: an overview[J]. Frontiers in Microbiology, 2017, 8: 1878. DOI:10.3389/fmicb.2017.01878
[6]
Baranova N, Nikaido H. The baeSR two-component regulatory system activates transcription of the yegMNOB(mdtABCD)transporter gene cluster in Escherichia coli and increases its resistance to novobiocin and deoxycholate[J]. Journal of Bacteriology, 2002, 184(15): 4168-4176.
[7]
Guerrero P, Collao B, Morales EH, et al. Characterization of the BaeSR two-component system from Salmonella Typhimurium and its role in ciprofloxacin-induced mdtA expression[J]. Archives of Microbiology, 2012, 194(6): 453-460. DOI:10.1007/s00203-011-0779-5
[8]
Lin MF, Lin YY, Yeh HW, et al. Role of the BaeSR two-component system in the regulation of Acinetobacter baumannii adeAB genes and its correlation with tigecycline susceptibility[J]. BMC Microbiology, 2014, 14(1): 119.
[9]
Lin YH, Yang TC, Hu WS. Characterization of stm3030 and stm3031 genes of Salmonella enterica serovar Typhimurium in relation to cephalosporin resistance[J]. Journal of Microbiology, Immunology and Infection, 2019, 52(2): 282-288. DOI:10.1016/j.jmii.2018.09.008
[10]
Rosner JL, Martin RG. Reduction of cellular stress by TolC-Dependent efflux pumps in Escherichia coli indicated by BaeSR and CpxARP activation of spy in efflux mutants[J]. Journal of Bacteriology, 2013, 195(5): 1042-1050. DOI:10.1128/JB.01996-12
[11]
王文静, 徐军, 高海侠, 等. 鼠伤寒沙门菌baeSR基因缺失株的构建及其对抗生素敏感性分析[J]. 微生物学通报, 2019. DOI:10.13344/j.microbiol.china.181072
[12]
裘鹏, 郭玉茹, 张泽辉, 等. 基于转录组学技术分析细菌耐药机制的研究进展[J]. 中国家禽, 2019, 41(4): 45-49.
[13]
Correia S, Michel Hébraud, Patrícia Poeta, et al. How combined multicomparative proteomic approaches can improve the understanding of quinolone resistance in Salmonella Typhimurium[J]. Future Microbiology, 2018, 13(4): 403. DOI:10.2217/fmb-2017-0250
[14]
Reygaert WC. An overview of the antimicrobial resistance mechanisms of bacteria[J]. AIMS Microbiology, 2018, 4(3): 482-501.
[15]
Buckner MMC, Blair JMA, et al. Beyond antimicrobial resistance: evidence for a distinct role of the AcrD efflux pump in Salmonella biology[J]. MBio, 2016, 7(6): e01916.
[16]
Wood CR, Mack LE, Actis L A. An update on the acinetobacter baumannii regulatory circuitry[J]. Trends in Microbiology, 2018, 26(7): 560-562. DOI:10.1016/j.tim.2018.05.005
[17]
Nishino K, Nikaido E, Yamaguchi A. Regulation of multidrug efflux systems involved in multidrug and metal resistance of Salmonella enterica serovar typhimurium[J]. Journal of Bacteriology, 2007, 189(24): 9066-9075. DOI:10.1128/JB.01045-07
[18]
Steinbuch KB, Fridman M. Mechanisms of resistance to membrane-disrupting antibiotics in Gram-positive and Gram-negative bacteria[J]. Med ChemComm, 2016, 7(1): 86-102.
[19]
Wolska KI, Grudniak AM, Rudnicka Z, et al. Genetic control of bacterial biofilms[J]. J Appl Genet, 2015, 21: 1-14.
[20]
Seo SW, Gao Y, Kim D, et al. Revealing genome-scale transcrip-tional regulatory landscape of OmpR highlights its expanded regulatory roles under osmotic stress in Escherichia coli K-12 MG1655[J]. Scientific Reports, 2017, 7(1): 2181. DOI:10.1038/s41598-017-02110-7
[21]
Newman SL, Will WR, et al. The curli regulator CsgD mediates stationary phase counter-silencing of csgBA in Salmonella Typhim-urium[J]. Molecular Microbiology, 2018, 108(1): 101-114. DOI:10.1111/mmi.13919
[22]
朱春红, 董立伟, 张江英, 等. 鞭毛素fliC基因的缺失对肠炎沙门氏菌生物膜形成的影响[J]. 中国预防兽医学报, 2016, 38(4): 303-306. DOI:10.3969/j.issn.1008-0589.2016.04.11
[23]
Wang Z, Bie P, Cheng J, et al. The ABC transporter YejABEF is required for resistance to antimicrobial peptides and the virulence of Brucella melitensis[J]. Scientific Reports, 2016, 6: 31876. DOI:10.1038/srep31876