Benchmark of 16S rRNA gene amplicon sequencing using ...

As the official Illumina protocol adopted the V3–V4 (V34) region, these two regions are widely used in gut microbiota studies [11,12,13,14,15,16 ... Skiptomaincontent Advertisement SearchallBMCarticles Search Benchmarkof16SrRNAgeneampliconsequencingusingJapanesegutmicrobiomedatafromtheV1–V2andV3–V4primersets DownloadPDF DownloadPDF Research OpenAccess Published:10July2021 Benchmarkof16SrRNAgeneampliconsequencingusingJapanesegutmicrobiomedatafromtheV1–V2andV3–V4primersets ShoichiroKameoka1,2,DaisukeMotooka1,3,4,SatoshiWatanabe2,RyuichiKubo2,NicolasJung1,YukiMidorikawa2,NatsukoO.Shinozaki2,YuSawai2,AyaK.Takeda2&ShotaNakamura1,3,4,5  BMCGenomics volume 22,Article number: 527(2021) Citethisarticle 3777Accesses 12Citations 11Altmetric Metricsdetails AbstractBackground16SrRNAgeneampliconsequencing(16Sanalysis)iswidelyusedtoanalyzemicrobiotawithnext-generationsequencingtechnologies.Here,wecomparedfecal16Sanalysisdatafrom192JapanesevolunteersusingthemodifiedV1–V2(V12)andthestandardV3–V4primer(V34)setstooptimizethegutmicrobiotaanalysisprotocol.ResultsQIIME1andQIIME2analysisrevealedahighernumberofunclassifiedrepresentativesequencesintheV34datathanintheV12data.Thecomparisonofbacterial compositiondemonstratedthatatthephylumlevel,ActinobacteriaandVerrucomicrobiaweredetectedathigherlevelswithV34thanwithV12.Amongthesephyla,weobservedhigherrelativecompositionsofBifidobacteriumandAkkermansiawithV34.Toestimatetheactualabundance,weperformedquantitativereal-timepolymerasechainreaction(qPCR)assaysforAkkermansiaandBifidobacterium.WefoundthattheabundanceofAkkermansiaasdetectedbyqPCRwasclosetothatinV12data,butwasmarkedlylowerthanthatinV34data.TheabundanceofBifidobacteriumdetectedbyqPCRwashigherthanthatinV12andV34data.ConclusionsTheseresultsindicatethatthebacterialcompositionderivedfromtheV34regionmightdifferfromtheactualabundanceforspecificgutbacteria.WeconcludethattheuseofthemodifiedV12primersetismoredesirableinthe16SanalysisoftheJapanesegutmicrobiota. PeerReviewreports BackgroundThe16SrRNAgene(16S)isconservedinmostbacteriaandarchaeaandisapproximately1,500 bplongwithninedifferenthypervariableregions(V1–V9).Ampliconsequencingtargetingthe16SrRNAgene(16Sanalysis)iswidelyusedtoanalyzemicrobiotausingnext-generationsequencing(NGS)technologies[1,2].Becauseofthelimitationofshort-readNGStechnologies,variousuniversalprimerstargetingthepartialsequencesinhypervariableregions(e.g.V1–V2,V1–V3,V3–V4,V4,etc.)havebeendevelopedformicrobiomeanalysis.Furthermore,manyexperimentalandanalyticalvariationsinthe16Sanalysisprotocolhavebeenreported[3,4].Currently,thereisnogoldstandardmethod,andthestandardizationofvariousapplicationsisanimportantissueinmetagenomicanalysis.Forexample,theDNAextractionprocessisawell-knowninfluencingfactorformetagenomicanalysis.Carefulcomparisonsbetweenthebead-beatingandenzymaticlysismethodshaveindicatedhowthisextractionprocessaffectsthemetagenomicdata[5,6,7].Oneofthemostcontroversialissuesin16Sanalysisprotocolvariationsistheselectionofthehypervariableregion(s)totarget.Historically,theV1–V2(V12)regionhasbeenemployedinmanyreportsofgutmicrobiotausingpastNGStechnologies,forexample,454pyrosequencing[8,9,10].AstheofficialIlluminaprotocoladoptedtheV3–V4(V34)region,thesetworegionsarewidelyusedingutmicrobiotastudies[11,12,13,14,15,16,17,18].Claessonetal.showedthattheV34primer-paircombinationcausesamplificationartifactsandhasadeviatingcompositioncomparedtootherregions,includingV12[19].Incontrast,Chenetal.reportedthatV34ismoresuitableforgutmicrobiotaanalysisthanV12becauseithasahigherpotentialtodetecttheorderBifidobacteriales[20].However,Kimetal.havedevelopedanewversionoftheV1forwardprimer(27Fmod)withimprovedBifidobacterium-detectioncomparedtothatoftheprimerusedbyChenetal.(27F-YM)[21].Thisproblemdoesnotonlyexistforthegutmicrobiota.ComparativestudiesofenvironmentalwaterreportedthattheV34orV4primersetsareoptimal[22,23].Asimilarconclusion,theV34regionismoresuitablethanV12,wasreportedforhumanvaginalmicrobiomeanalysis[24].Meanwhile,acomparativestudyofthehumanoralmicrobiomereportedthatV1–V3ismoresuitablethanV34[25].Asshownintheseexamples,thechoiceofoptimalvariableregionsmightvarydependingontheanalysistarget,thespecificityoftheprimers,theGCcontentsintheselectedregion,andthebacterialcompositionsofdifferentsamples.Dataprocessingtoolsareanotherfactorthatimpactstheinterpretationofthemicrobiome.QIIME(QuantitativeInsightsIntoMicrobialEcology),oneofthemostpopularbioinformaticstoolsfor16Sanalysis,comprisesversion1(QIIME1:qI)releasedin2011andversion2(QIIME2:qII)releasedin2018[26,27].qIIhasbeencompletelyredesignedasanalgorithmforclusteringampliconsequencevariants(ASVs)fromoperationaltaxonomicunits(OTUs)inthepreviousversionandisreportedtoallowformoreaccurateclusteringofASVs.However,insomecases,qIisstillbeingutilizedtocomparethevastamountofdatathathavebeenproducedinthepast.Inthisstudy,wefocusedonselectingthehypervariableregionsof16SrRNAforhumangutmicrobiomeanalysis.Toourknowledge,therehavebeennoreportscomparingthegutmicrobiotabetweentheV12regionwith27FmodandV34region.ToidentifywhichprimersetthatisbettersuitedforanalyzingtheJapanesegutmicrobiome,wecomparedtheintestinalmicrobiomedatafrom192Japanesevolunteersusingbothprimersets,V12(27Fmod)andV34.Inaddition,toevaluatethosedatawiththepreviouslypublishedstudies,acomparisonbetweenQIIMEversionsqIandqIIwasalsoconducted.Participants&methodsParticipantsAtotalof192healthyJapanesevolunteerswasrandomlyselectedfromthe1,644samplesintheMykinsolibraryrecruitedbetweenJuly2015andAugust2016.TheageandgenderdistributionsoftheselectedparticipantsareshowninTable 1S.Fecalsampling,DNAextraction,andsequencingWeperformedtheV1–V2regionsequencingasreportedinWatanabeetal.[28].Wecollectedfecalsamplesusingbrush-typecollectionkitscontainingguanidinethiocyanatesolution(TechnoSurugaLaboratory,Shizuoka,JPN),transportedatambienttemperature,andstoredat4 °C.DNAwasextractedfromthefecalsamplesusingtheDNeasyPowerSoilKit(QIAGEN,Hilden,DEU)accordingtothemanufacturer’sprotocol.TheampliconsofV12werepreparedusingtheforwardprimer(16S_27Fmod:TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAGRGTTTGATYMTGGCTCAG)andthereverseprimer(16S_338R:GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGTGCTGCCTCCCGTAGGAGT).TheampliconsofV34werepreparedusingtheforwardprimer(16S_341F:TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG)andthereverseprimer(16S_805R:GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC)withaKAPAHiFiHotStartReadyMix(Roche,Basel,CHE).Thesequencinglibrarieswerepreparedaccordingtothe16SlibrarypreparationprotocolprovidedbyIllumina(Illumina,SanDiego,CA,USA).DualindexadaptersforsequencingontheIlluminaMiSeqplatformwereattachedusingtheNexteraXTIndexkit(Illumina,SanDiego,CA,USA).Eachsequencinglibrarywasdilutedto5ng/µL.Wemixedequalvolumesofthelibrariestogiveafinalconcentrationof4nM.TheDNAconcentrationofthemixedlibrarieswasmeasuredbyquantitativereal-timepolymerasechainreaction(qPCR)withtheKAPASYBRFASTqPCRMastermix(Roche,Basel,CHE)usingprimer1(AATGATACGGCGACCACC)andprimer2(CAAGCAGAAGACGGCATACGA).Theselibrariesweresequencedina250-bppaired-endrunforV12usingtheMiSeqReagentKitv2(500cycles),andina300-bppaired-endrunforV34usingtheMiSeqReagentKitv3(600cycles).BioinformaticsanalysisTheQIIME1analysiswasperformedasreportedinWatanabeetal.[28]andthedetailsaredescribedinthefollowingstrategy:thepaired-endsequencingreadswereclusteredby97 %nucleotideidentity,andtaxonomicinformationwasassignedusingtheGreengenesdatabase(v13.8)[29]usingtheQIIMEpipeline(v1.8.0)[26].Thedataprocessingandassignmentwereperformedthefollowingsteps:(1)joiningofpaired-endreads;(2)qualityfilteringwithanaccuracyofQ30(> 99.9 %)andareadlength > 300 bp;(3)clusteringofOTUswith97 %identityusingUCLUST(v1.2.22q)[30];(4)assigningtaxonomicinformationtoeachOTUusingtheRDPclassifier[31]withthefull-length16SrRNAgenedatafromGreengenes(v13.8)todeterminetheidentityandbacterialcomposition.TheQIIME2analysiswasperformedusingthefollowingstrategy:constructanASVtablewithpaired-endsequencingreads,andthen,taxonomicinformationwasassignedusingGreengenes[29]usingtheQIIME2pipeline(version2020.2)[27].ThedataprocessingandassignmentbasedontheQIIME2pipelinewereperformedthefollowingsteps:(1)DADA2[32]forjoiningpaired-endreads,filtering,anddenoising;(2)assigningtaxonomicinformationtoeachASVusinganaiveBayesclassifierintheQIIME2classifierwiththe16SrRNAgeneV2regiondata forV12sequencingandV34regiondataforV34sequencingfromGreengenestodeterminetheidentityandbacterialcomposition.Totesttwo-groupdifferencesinthepercentageofanalyzablereadnumbersbetweenV12andV34,wecalculatedp-valuesusingtheWilcoxonsigned-ranktest.MicrobiomeanalysisAlphadiversitywasassessedatadepthof30,000readsusingtheOTU/ASVtablesderivedfromtheqIandqIIanalysis,respectively.Then,theWilcoxonsigned-ranktestwasperformedtotestthetwogroupsdiversitydifferences.BasedontheOTU/ASVtablesatthegenuslevel,Bray-CurtisdissimilaritybetweentheV12andV34regionswascalculated.ToassessthebetadiversitiesoftheV12andV34data,weperformedapermutationalmultivariateanalysisofvariance(PERMANOVA).TheWilcoxonsigned-ranktestwasperformedtotestdifferencesinrelativecompositionatthephylumandgenuslevelsbetweentheV12andV34regions,wecalculatedp-valuesusingtheWilcoxonsigned-ranktest.qPCRanalysisqPCRwasperformedusingKAPASYBRFastandLightCycler480SystemII(RocheDiagnosticsK.K.,Rotkreuz,CHE.)underthefollowingconditions:95 °Cfor30 s(95 °Cfor5 s,60 °Cfor30 s)×45cycles.TheprimersusedwereBacteroides(forward:CAATCGGAGTTCTTCGTGATATCTA;reverse:GTTGTGAAAGTTTGCGGCTCA),Faecalibacterium(forward:TGTAAACTCCTGTTGTTGAGGAAGATAA;reverse:GCGCTCCCTTTACACCCA),Bifidobacterium(forward:CGCGTCYGGTGTGAAAG;reverse:CCCCACATCCAGCATCCA),Akkermansia(forward:CTGAAGAACTCGGCACCCTT;reverse:CTTCTTCAGCTTCGGCAGGA),andbacterial16SrRNA(forward:ACTCCTACGGGAGGCAGCAGT;reverse:TATTACCGCGGCTGCTGGC).Tocalculatetherelativeabundanceofeachbacteria,thedatawerenormalizedbysubtractingthe16SrRNAcyclethreshold(Ct)valueforeachrespectivesamplefromtheCtvaluesforthetargetbacteriatocalculateΔCtvalues,whichareexpressedas2^[Ct(16SPCR)-Ct(targetPCR)],respectively[33,34].ResultsComparisonofthe16SProfilesUsingV1–V2andV3–V4PrimersWerecruited192volunteers,collectedtheirstoolsamples,andstoredtheminaguanidineisothiocyanate-basedreagentatambienttemperature.DNAextractionwascarriedoutusingabead-basedmethodaspreviouslyreportedinsimilarstudies[5,6,7].ThestandardprotocolwasusedforthepreparationofanIlluminalibrarytoallowforthesequencingof16Samplicons.Paired-endsequencingwasperformedwitha250-bplengthfortheV1–V2region(V12)andwitha300-bplengthfortheV3–V4region(V34).TheaveragereadcountsforV12andV34were44,442and47,220,respectively(Table 1).TheaveragesequencequalitiesofV12wereaboveQ30by249 bpforread1andby220 bpforread2.ThoseofV34wereaboveQ30by244 bpforread1andby204 bpforread2(Fig. 1S-A).WeanalyzedtheserawsequencedatawiththestandardqIandqIIpipelines.ThepreprocessingforqI(merge,adaptertrim,andprimercheckstep)yieldedanalyzablereadcountsforV12andV34of40,941(91.3 %;V12qI)and37,611(79.6 %;V34qI),respectively(Table 1;Fig. 1S).Incontrast,thepreprocessingforqII(filtering,denoising,merging,non-chimericstep)yieldedanalyzablereadcountsof35,609(80.1 %;V12qII)and31,014(66.0 %;V34qII),respectively.AlthoughtheanalyzablereadsforV34werefewerthanthoseforV12(p = 7.46e-62(qI).p = 1.27e-60(qII).Fig 1S-B,C),bothdatasetssatisfiedthetypicalreadamount(1,000–50,000reads)aftermergingandqualityfiltering(Table 1)[35].ThequalityofthemergedreadsdidnotdiffergreatlybetweenV12andV34(Fig. 1S-D).Basedontheseresults,weconcludedthatthesequencingdataobtainedwereofsufficientqualityandquantityformicrobiomeanalysisusingbothversionsoftheQIIMEpipeline. Table116SsequencingprofilecomparisonbetweenV12andV34FullsizetableDifferenceinobservedOTU/ASVWenextevaluatedthealphadiversitiesoftheV12andV34data.AsshowninFig. 1,at30,000reads,themedianOTUcountforV12qI(1,078)waslowerthanthatforV34qI(1,382.8).Incontrast,themedianvalueofASVcountforV12qII(238.65)washigherthanthatofV34qII(147).Tounderstandthereasonforthisdifference,weinvestigatedeachOTU/ASVassignedbyqIandqII(Fig. 2SandTable2S).AssignedOTU/ASVratesatthephylumlevelwere99.46 %(V12qI),97.76 %(V34qI),99.84 %(V12qII),and98.67 %(V34qII).Atthegenuslevel,theseratesdecreasedto69.92 %(V12qI),64.49 %(V34qI),63.04 %(V12qII),and66.51 %(V34qII).ThemajorphylaassignedwereFirmicutes,Bacteroidetes,Proteobacteria,andActinobacteria,asiscommonlyobservedinnormalgutflora.WecomparedbothindividualandtotalOTUcountsforthesemajorphylaandfoundthatintheoutputfromqI,theOTUcountsofBacteroidetes,Proteobacteria,andActinobacteriaforV34werehigherthanthoseforV12(Fig.3S-A).Conversely,intheoutputfromqII,theASVcountsofallmajorphylaforV34werelowerthanthoseforV12(Fig. 3S-B).Otherthanthesemajorphyla,therewereunclassifiedOTUs/ASVswherethehighesttaxonomiclevelassignedwasonlytothekingdomsArchaeaorBacteria.TheuseofqIIreducedthecountsoftheseunclassifiedASVstonearlyzerocomparedtotheunclassifiedOTUsofqI(Fig. 3S).Whetherclassifiedorunclassified,alotoftherepresentativesequencesofOTUsarefilteredoutbyDADA2.Welabeledthosesequencesasfiltered-outsequences.Thepercentageoffiltered-outsequenceswashigherinV34(94.3 %)thanV12(87.9 %)(Table3S).TheseresultsindicatethatboththeV12andV34dataweregreatlyaffectedbytheversionofQIIMEusedforanalysisandthatthecombinationoftheV34primersetandqIIprocessingresultedinthedetectionoffewerASVscomparedtothatwithV12qII. Fig.1Observedoperationaltaxonomicunits(OTUs)andampliconsequencevariants(ASVs)plotateachsequencingdepthforV12(red)andV34(blue).A qIandB qII.Statisticalanalysis(Wilcoxonsigned-ranktest)wasperformedatadepthof(A:qI)29,998and(B:qII)30,000.Doubleasterisksindicatestatisticalsignificance(p 28(green),passable:28>Q>20(orange)andpoor:20>Q(red).(B,C)Theratioofthenumberofanalyzablereadstothenumberofrawreadsisshown.Doubleasterisksindicatestatisticalsignificance(p<0.01).Figure2S.Percentageofclassifiedoperationaltaxonomicunits(OTUs)andampliconsequencevariants(ASVs).Ateachclassificationlevel,theresultsofcalculatingthepercentageofOTUs/ASVsassignedtoataxonomyinrelationtothetotalnumberofOTUs/ASVsareshown.Figure3S.(A,B)Differencesintheoperationaltaxonomicunit(OTU)andampliconsequencevariant(ASVs)countsfortheindicatedphylumbetweentheV12andV34regions.Doubleasterisksindicatestatisticalsignificance(p<0.01).(C,D)ComparisonofthetotalOTU/ASVnumbersfortheindicatedphylum.UnclassifiedOTUs/ASVsareincludedink__Bacteria;__,k__Bacteria;p__,k__Archaea;__,andUnclassified;__.Figure4S.BarchartoftheindividualbacterialcompositionsusingV12andV34atthephylumlevelusing(A)qI(upperpanel:V12;lowerpanel:V34)and(B)qII(upperpanel:V12;lowerpanel:V34).Figure5S.RelativecompositionofBacteroidetes,Firmicutes,andProteobacteriausingV12andV34.Figure6S.BarchartoftheindividualbacterialcompositionsusingV12andV34fortheindicatedphylausing(A)qI(upperpanel:V12.lowerpanelV34)and(B)qII(upperpanel:V12;lowerpanel:V34).Figure7S.BarchartoftheindividualbacterialcompositionsusingV12andV34fortheindicatedgenerausing(A)qI(upperpanel:V12;lowerpanel:V34)and(B)qII(upperpanel:V12;lowerpanel:V34).Figure8S.BarchartofbacterialrelativeabundanceusingaDNAmockcommunitykindlyprovidedbyNITE(NationalInstituteofTechnologyandEvaluation,Tokyo,JPN).Thiscommunityismadefromanequalmixofgenomicdatafromthe10indicatedstrains.qPCRwasperformedtargetingtherplLgeneofeachbacteriaandnormalizedbytotalbacteriausingthemeasuredvalueof16SrRNAgene.rRNAcopynumindicatesthepercentageofcopynumberofthe16SrRNAgene.V12andV34indicatetheresultsof16Sanalysis.Figure9S.Scatterplotoftheindicatedgeneratocomparethe16Sanalysisby(A,B)qI,(C,D)qII,andqPCR.Theidentityline(y=x)isindicated.Thegrayareaindicatesthe95%confidenceintervalforeachregressionline.Figure10S.Scatterplotoftheindicatedgeneratocomparethe16SanalysisbyqIIandqPCR.Theidentityline(y=x)isindicated.Thegrayareaindicatesthe95%confidenceintervalforeachregressionline.Figure11S.AlignmentsandsimilarityofCronobacterandOTUrepresentativesequencewithAkkermansiaforthe16SrRNAgene.(A)V12(B)V34of16SrRNAgene.Akkermansiaindicatesreferencesequences,whicharederivedfromAkkermansiamuciniphilastrainJCM30893.CronobactersequencesarederivedfromCronobactersakazakiistraincro360A2.Observed_OTUsequenceisarepresentativesequenceassignedtoAkkermansia,whichwasmatchedtoCronobacterthroughBLAST.ThevaluesofsimilarityindicatepercentidentitybetweenthereferencesequenceandthequerysequencecalculatedbyBLAST.Additionalfile2:Table1S.Distributionoftheselectedparticipants. Table2S.PercentageofclassifiedOTUs/ASVsateachclassificationlevel.Table3S.OTUs/ASVsnumbersandthepercentageoffiltered-outsequences.Table4S.ListsofbacteriathatshowedstatisticaldifferencesatthephylumlevelbetweenV12andV34.Table5S.ListsofbacteriathatshowedstatisticaldifferencesatthegenuslevelbetweenV12andV34.Table6S.AveragecompositionsofthegeneraBacteroidesandFaecalibacterium.Additionalfile3.RelativeabundanceofAkkermansia,Bifidobacterium,Bacteroides,andFaecalibacteriumbyqPCR.Tocalculatetherelativeabundanceofeachbacteria,thedatawerenormalizedbysubtractingthe16SrRNAcyclethreshold(Ct)valueforeachrespectivesamplefromtheCtvaluesforthetargetbacteriatocalculateΔCtvalues,whichareexpressedas2^[Ct(16SPCR)-Ct(targetPCR)],respectively.Rightsandpermissions OpenAccessThisarticleislicensedunderaCreativeCommonsAttribution4.0InternationalLicense,whichpermitsuse,sharing,adaptation,distributionandreproductioninanymediumorformat,aslongasyougiveappropriatecredittotheoriginalauthor(s)andthesource,providealinktotheCreativeCommonslicence,andindicateifchangesweremade.Theimagesorotherthirdpartymaterialinthisarticleareincludedinthearticle'sCreativeCommonslicence,unlessindicatedotherwiseinacreditlinetothematerial.Ifmaterialisnotincludedinthearticle'sCreativeCommonslicenceandyourintendeduseisnotpermittedbystatutoryregulationorexceedsthepermitteduse,youwillneedtoobtainpermissiondirectlyfromthecopyrightholder.Toviewacopyofthislicence,visithttp://creativecommons.org/licenses/by/4.0/.TheCreativeCommonsPublicDomainDedicationwaiver(http://creativecommons.org/publicdomain/zero/1.0/)appliestothedatamadeavailableinthisarticle,unlessotherwisestatedinacreditlinetothedata. ReprintsandPermissionsAboutthisarticleCitethisarticleKameoka,S.,Motooka,D.,Watanabe,S.etal.Benchmarkof16SrRNAgeneampliconsequencingusingJapanesegutmicrobiomedatafromtheV1–V2andV3–V4primersets. BMCGenomics22,527(2021).https://doi.org/10.1186/s12864-021-07746-4DownloadcitationReceived:11March2021Accepted:25May2021Published:10July2021DOI:https://doi.org/10.1186/s12864-021-07746-4SharethisarticleAnyoneyousharethefollowinglinkwithwillbeabletoreadthiscontent:GetshareablelinkSorry,ashareablelinkisnotcurrentlyavailableforthisarticle.Copytoclipboard ProvidedbytheSpringerNatureSharedItcontent-sharinginitiative KeywordsMicrobiota16SrRNANext-generationsequencing DownloadPDF Advertisement BMCGenomics ISSN:1471-2164 Contactus Submissionenquiries:[email protected] Generalenquiries:[email protected]