Transcriptome Sequencing and Positive Selected Genes Analysis of Bombyx mandarina

14 Pages·2015·0.59 MB·English

RESEARCHARTICLE
Transcriptome Sequencing and Positive
Bombyx
Selected Genes Analysis of
mandarina
TingcaiCheng☯,BohuaFu☯,YuqianWu,RenwenLong,ChunLiu,QingyouXia*
StateKeyLaboratoryofSilkwormGenomeBiology,SouthwestUniversity,Chongqing,China
☯Theseauthorscontributedequallytothiswork.
* [email protected]
a11111
Abstract
ThewildsilkwormBombyxmandarinaiswidelybelievedtobeanancestorofthedomesti-
catedsilkworm,Bombyxmori.Silkwormsareoftenusedasamodelforstudyingthemecha-
OPENACCESS nismofspeciesdomestication.Here,weperformedtranscriptomesequencingofthewild
silkwormusinganIlluminaHiSeq2000platform.Weproduced100,004,078high-quality
Citation:ChengT,FuB,WuY,LongR,LiuC,XiaQ
(2015)TranscriptomeSequencingandPositive readsandassembledtheminto50,773contigswithanN50lengthof1764bpandamean
SelectedGenesAnalysisofBombyxmandarina. lengthof941.62bp.Atotalof33,759unigeneswereidentified,with12,805annotatedinthe
PLoSONE10(3):e0122837.doi:10.1371/journal.
Nrdatabase,8273inthePfamdatabase,and9093intheSwiss-Protdatabase.Expression
pone.0122837
profileanalysisfoundsignificantdifferentialexpressionof1308unigenesbetweenthemid-
AcademicEditor:ErjunLing,InstituteofPlant dlesilkgland(MSG)andposteriorsilkgland(PSG).Threesericingenes(sericin1,sericin
PhysiologyandEcology,CHINA
2,andsericin3)wereexpressedspecificallyintheMSGandthreefibroingenes(fibroin-H,
Received:November3,2014 fibroin-L,andfibroin/P25)wereexpressedspecificallyinthePSG.Inaddition,32,297Sin-
Accepted:February15,2015 gle-nucleotidepolymorphisms(SNPs)and361insertion-deletions(INDELs)weredetected.
Published:March25,2015 Comparisonwiththedomesticatedsilkwormp50/Dazaoidentified5,295orthologous
genes,amongwhich400mighthaveexperiencedortobeexperiencingpositiveselection
Copyright:©2015Chengetal.Thisisanopen
accessarticledistributedunderthetermsofthe byKa/Ksanalysis.Thesedataandanalysespresentedhereprovideinsightsintosilkworm
CreativeCommonsAttributionLicense,whichpermits domesticationandaninvaluableresourceforwildsilkwormgenomicsresearch.
unrestricteduse,distribution,andreproductioninany
medium,providedtheoriginalauthorandsourceare
credited.
DataAvailabilityStatement:Allrawdataare
availablefromtheGenBankdatabase(accession Introduction
numbers:SRX738967andSRX738979).
ThewildsilkwormBombyxmandarinabelongstoLepidopteraBombycidae.Itiswidelyaccept-
Funding:Thisworkwassupportedbygrantsfrom
edasthatB.mandarinaiswidelyacceptedasanancestorofthedomesticatedsilkworm,Bom-
theNationalBasicResearchProgramofChina(973
program,2012CB114600),theNationalHigh-tech byxmori.B.mandarinahavebeendomesticatedforatleast5000yearstoincreaseandimprove
R&DProgram(863program,2011AA100306),the cocoonyield.Despitetheshortdomesticationhistory,manytraitsaredifferentbetweenB.
NationalNaturalScienceFoundationofChina(no. mandarinaandB.moriincludingbodysizeandcolorinthelarvalstage,sizeandsilkqualityof
31001034),andtheFundamentalResearchFunds
cocoons,fightbehaviorandegglayingintheadultstage.Thus,B.mandarinaandB.moriare
fortheCentralUniversities(XDJK2013C004).The
goodmodelsforstudyingspeciesdomestication.
fundershadnoroleinstudydesign,datacollection
Recently,basedonnext-generationsequencingresults,somestudieshavefoundthatpheno-
andanalysis,decisiontopublish,orpreparationof
themanuscript. typicchangesareinvolvedingeneticdivergences.Bombyxresequencingidentifiedabout16
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 1/14
TranscriptomeAnalysisoftheWildSilkworm
CompetingInterests:Theauthorshavedeclared millionsingle-nucleotidepolymorphisms(SNPs),311,608indelsand35,093structuralvaria-
thatnocompetinginterestsexist. tions(SVs)indomesticatedandwildsilkworms.Atotalof1041regions,covering2.9%ofthe
genomeand354protein-codinggenesweredetectedwithselectionsignals[1].Inaddition,347
SNPswereidentifiedinthesilkwormmitochondriagenome.Acytochromebgeneshows
strongpositiveselectioninthedomesticatedgroup[2].Copynumbervariation(CNV)isalso
animportantdomesticationmechanisminsilkworms.About319CNVshavebeenidentified
andabout49%aredistributedonuncharacterizedchromosomes.Approximately61%of
CNVsdirectlyoverlapwithsegmentalduplicationsinsilkworms.ThegenesinCNVsaremain-
lyrelatedtoreproduction,immunity,detoxificationandsignalrecognition[3].Forexample,
thecopynumberofcarotenoid-bindingprotein(CBP)genevariesfrom1to20amongB.mori
strains.IncontrasttoB.mori,B.mandarinawasfoundtopossessasinglecopyofCBPlacking
aretrotransposoninsertion,regardlessofhabitat.TheCBPgeneisevolutionarilyconservedin
thelepidopteranlineage,showingthatdomesticationcangeneratesignificantdiversityingene
copynumberandstructureoverarelativelyshortevolutionarytime[4].Thedomestication
mechanismmayinvolveregulatoryelements.Aputativeexpressionenhancerlocatedinthein-
tronoftheThgeneregulatestranscription.Thesitemightcontributetothebodycolortransi-
tionfromB.mandarinatoB.mori[5].Inadditiontogeneticdivergence,epigeneticdivergence
isalsoimportantinsilkwormdomestication.Comparativemethylomicsbetweendomesticated
andwildsilkwormsshowedtwiceasmanymethylatedcytosinesindomesticatedsilkwormsas
intheirwildcounterparts,suggestingatrendtowardincreasingDNAmethylationduringdo-
mestication.Genesthataredemethylatedandhavelowexpressionindomesticatedsilkworms
haveexperiencedselectivesweep,indicatingapossiblecorrelationwiththeenlargementofsilk
glandsindomesticatedsilkworms[6].
Next-generationsequencing-basedRNA-Seqanalysisprovidesopportunitiesfordenovoas-
semblyofgenomereference-freespecies[7].Thismethodprovidesinformationongeneex-
pression,generegulationandspeciesevolution.Transcriptomeanalysishasbeenwidely
reportedinsomespeciessuchasFormicaexsecta[8],Deliaantiqua[9],andPolistescanadensis
[10].Recently,RNA-Seqhasalsobeenusedasanefficientmethodtostudyadaptationtohigh-
elevationenvironmentsordomestication.Basedondenovoassembledtranscriptsandidentifi-
cationoforthologousgenes,nonsynonymoussite/synonymoussite(Ka/Ks)analysiscanpro-
videinsightsintotheprocessofadaptiveevolutionordomestication[11,12].
Inthepresentstudy,weperformedtranscriptomesequencingforthemiddlesilkgland
(MSG)andtheposteriorsilkgland(PSG)fromB.mandarinausinganIlluminaHiSeq2000se-
quencingplatform.Bydenovoassembly,wegeneratedanumberofunigenes.Bycomparison
withdomesticatedsilkworm(p50/Dazao)genesandKa/Ksanalysis,weidentifiedorthologous
genesandgenesthatmighthaveexperiencedortobeexperiencingpositiveselection.Function-
alannotationofthesegenesprovidedinformationaboutthedomesticationmechanism
insilkworms.
MaterialsandMethods
SamplecollectionandRNAextraction
Samplesofwildsilkwormswerecollectedatthesuburb(Lat/Lon:29.60°N106.28°W)of
Chongqingcity,China.Nospecificpermissionswererequiredfortheselocations.Thefield
studies(Lat/Lon:29.60°N106.28°W)didnotinvolveendangeredorprotectedspecies.Wedis-
sectedandcollectedtwotissues:middlesilkglands(MSGs)andposteriorsilkglands(PSGs)
fromasinglefifth-instarlarva.Tissueswereimmediatelyfrozenandstoredinliquidnitrogen.
TotalRNAswereextractedusingTRIzolReagent(Invitrogen)andtreatedwithDNase.The
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 2/14
TranscriptomeAnalysisoftheWildSilkworm
qualityandquantityofpurifiedRNAwereexaminedusingAgilentBioanalyzer2100(Agilent
Technologies)andQubitRNAAssayKit(Invitrogen,http://products.invitrogen.com).
LibraryconstructionforRNA-seqandsequencing
RNAlibrarieswereconstructedusingIlluminaTruSeqRNApreparationkitsfollowingthe
manufacturer’sinstructions.LibrarieswerecheckedandquantifiedusingAgilentBioanalyzer
2100(AgilentTechnologies)andQubitdsDNABRAssayKit(Invitrogen,http://products.
invitrogen.com).Librariesweresequencedfor100bppaired-endreadsusingtheIllumina
HiSeq2000platform.DataanalysisandbasecallingwereperformedbyIlluminainstrument
software.RawdatapresentedinthispublicationhavebeendepositedintheNCBIShortRead
Archive(http://www.ncbi.nlm.nih.gov/sra/)andareaccessiblethroughSRAaccessionnum-
bersSRX738967andSRX738979.
Sequencedataanalysisanddenovoassembly
AdaptersequenceswereremovedbyTrimomatic0.32(http://www.usadellab.org/cms/?page=
trimmomatic).Low-qualitysequencereadswereeliminatedbyscanningreadswitha4-base
slidingwindow,cuttingwhentheaveragequalityperbasedroppedbelow15[13].Resultingse-
quencereadsbelow50bpwereremoved.Beforeassembly,rawreadswereassessedforquality
usingFASTQC(v0.10.1)software.WealsoremovedribosomalRNAreadsbycomparingwith
ribosomalRNAdatacollectedusingBowtie2[14].Denovoassemblyofcleanedreadswascar-
riedoutwithTrinitysoftware(http://trinityrnaseq.sf.net).
Geneannotation
Openreadingframes(ORFs)werepredictedbytheTransDecodertool[15]usingdefaultpa-
rameters.AssembledunigeneswereannotatedusingTrinotateannotationpipeline(http://
trinotate.sourceforge.net/).AssembledunigeneswerefirstusedtoBlastPalignmentagainstthe
Swiss-protdatabase(downloaded03/07/2014),usingHMMER3.1[16]toidentifyproteindo-
mainsbysearchingthePfam_Adatabase(downloaded03/07/2014).Resultswereloadedusing
GeneOntology(GO)termsclassifiedwithWEGO[17].Assembledunigeneswerealsoaligned
tononredundant(nr)proteins(http://www.ncbi.nlm.nih.gov/)andCOGdatabases(http://
www.ncbi.nlm.nih.gov/COG/).WeusedtheonlineKEGGAutomaticAnnotationServer
(KAAShttp://www.genome.jp/kegg/kaas/)forKEGGpathway(http://www.genome.jp/kegg/)
annotationwithKO(KEGGOrthology)assignments.Pathwayenrichmentanalysisofpositive
selectiveunigeneswasperformedwithcalculatingmethodsasearlyreports[18].Pathways
withp-value(cid:1)0.01weredefinedasenriched.
Analysisofgeneexpression
Toanalyzetissue-specificexpressionpatterns,cleanedreadsweremappedtoassembleduni-
geneswithBowtieusingrun_RSEM_align_n_estimate.plfromTrinitysoftware.RSEM1.2.11
(http://deweylab.biostat.wisc.edu/rsem/)softwarewasusedtocalculateexpressionlevelsas
fragmentsperkilobaseofexonmodelpermillionmappedreads(FPKM).Unigenesdifferen-
tiallyexpressedbetweenMSGandPSGweredetectedwiththeedgeRbioconductorpackage
[19].Iffalsediscoveryrate(FDR)waslowerthan0.01,thep-valuewaslowerthan0.05andthe
highestFPKMofunigenewasfourtimesthelowestFPKM,aunigenewasconsidered
differentiallyexpressed.
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 3/14
TranscriptomeAnalysisoftheWildSilkworm
ReadsmappingandSNPidentification
CleanedreadsweremappedtotheassembledunigenesofthewildsilkwormbyTophat2
(http://ccb.jhu.edu/software/tophat/index.shtml).SAMtools(toolsforalignmentsinSAMfor-
mat)softwarewasusedtodetectputativeSNPs[20].Basequalitywasnolessthan20andcov-
eragewasnolessthan5forSNPdetection.
Analysisofpositiveselections
Areciprocalbesthitmethodwasusedtoidentifyputativeorthologsbetweentwospeciesusing
BLASTPofBLAST-2.2.28+[21].Wedownloadeddomesticatedsilkworm(p50/Dazao)CDS
sequencesfromSilkDBV2.0(http://silkworm.genomics.org.cn/).ORFsofassembledunigenes
werecomparedtosilkwormCDSsequencesusingBLASTPreciprocallytofindorthologous
pairs.CodingregionofbothsequenceswerealignedbyClustalw2.1(http://www.clustal.org/
clustal2/).Amanualcheckwasalsoconductedtocorrectpotentialerrors.Theratioofthe
numberofnonsynonymoussubstitutionspernonsynonymoussite(Ka)tothenumberofsyn-
onymoussubstitutionspersynonymoussite(Ks)wasusedtotestforpositiveselection.AKa/
Ksratiogreaterthan1wasevidenceofpositiveselectionandaKa/Ksratiolessthan1indicated
purifyingselection[22].KaKs_Calculator[23]softwarewasusedtoestimateKa,KsandKa/Ks
ratioforeachorthologouspairusingtheYNmethod[24].SequenceswithKs>0.1wereex-
cludedtoavoidpotentialparalogs[25].
Results
RNA-Seqanddenovoassembly
LibrariesofMSGandPSGweresequencedusingIlluminapaired-endsequencingtechnology.
Intotal,120,381,200rawreadsweregeneratedfromwildsilkworms(S1Table).Removing
adapters,low-qualitysequencesandribosomalRNAyielded100,004,078cleanreadswithGC
percentage46.50%(S1Table).Trinitysoftwarewasusedtoassemblecleanreadsinto50,773
transcriptswithN50length1764bpandmeanlength941.62bp.Clusteringgenerated33,759
unigeneswithN50length1437bpandmeanlength762.20bp.Table1showsanoverviewof
theassembledtranscriptsandunigenes.
Table1. SummaryoftranscriptsforBombyxmandarina.
Number(percentage)
Length(bp) Transcript(bp) Unigene(bp)
200–300 14,829(29.21%) 12,428(36.81%)
300–500 11,338(22.33%) 8,488(25.14%)
500–1000 9,588(18.88%) 5,694(16.87%)
1000–2000 8,332(16.41%) 4,117(12.20%)
2000+ 6,686(13.17%) 3,032(8.98%)
Totallength 47,808,990 25,730,955
Count 50,773 33,759
GCpercentage 39.91% 39.99%
Mediancontiglength 478 374
Averagecontig 941.62 762.20
N50length 1,764 1,437
doi:10.1371/journal.pone.0122837.t001
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 4/14
TranscriptomeAnalysisoftheWildSilkworm
Fig1.Geneontology(GO)annotation.
doi:10.1371/journal.pone.0122837.g001
SequenceAnnotation
Forfunctionalannotation,ORFswerepredictedfromassembledunigenesbytheTransDeco-
dertoolandannotatedusingTrinotateannotationpipelinewithPfam,UniProt/Swiss-Prot,
andGO.WealsoannotatedunigenesbyaligningwithNr,KEGG,andCOGdatabasesfor
12,805unigenessignificantlymatchedinNr,8273matchedinPfam,and9093similartopro-
teinsintheSwiss-Protdatabase(S2Table).E-valuedistributionshowedthat35.01%ofNrda-
tabaseannotatedunigeneshadanE-valueequalto0(S1Fig).
Functionalannotation
GOcategoriesarewidelyusedtoclassifygenefunctions[26].Atotalof9571unigeneshadaGO
annotationinthreemaingroups(S2Table,Fig1).Withincellularcomponent,14level-2catego-
rieswereidentified.Thetopthreegroupswerecell,cellpart,andorganelle.Withinmolecular
functions,15level-2categorieswereidentified.Thetopthreegroupswerebinding,catalyticac-
tivity,andtransporteractivity.Inbiologicalprocess,23level-2categorieswereidentified.The
topthreegroupswerecellularprocess,metabolicprocess,andbiologicalregulation.
Atotalof6245unigeneswereclassifiedinto25functionalcategories(S2Table,Fig2).Inad-
ditiontothelargestgroupofgeneralfunctionpredictiononly(20.45%),thetopfourgroups
werereplication,recombination,andrepair(8.98%),posttranslationalmodification,protein
turnover,chaperones(7.99%),translation,ribosomalstructureandbiogenesis(7.29%),tran-
scription(6.52%),andaminoacidtransportandmetabolism(5.75%).
Atotalof5893unigeneswereassignedto253biologicalpathwaysintheKEGGdatabase.
Predictedpathwaysweredividedintofivebiologicalclasses(S2Fig):organismalsystems
(29.66%),metabolism(24.90%),environmentalinformationprocessing(17.53%),geneticin-
formationprocessing(14.77%),andcellularprocesses(13.14%).
Geneexpressionanalysis
Atotalof1308unigenesweredifferentiallyexpressed(FDR(cid:1)0.01,p-value(cid:1)0.05,absolute
valueoflogFC(cid:3)2)betweenMSGandPSGinwildsilkworms(S3Fig)with883up-regulated
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 5/14
TranscriptomeAnalysisoftheWildSilkworm
Fig2.Clustersoforthologousgroup(COG)annotation.
doi:10.1371/journal.pone.0122837.g002
expressedinMSGand425up-regulatedinPSG(S3Table).Thesilkglandistheonlylocation
thatsynthesizesandsecretessilkproteins.TheMSGsynthesizessericinandthePSGsynthe-
sizesfibroinproteins(Table2).Threeunigenesencodingsericinproteinsweretissue-specific
totheMSG.Fourunigenesencodingfibrohexamerin-likeproteinsweretissue-specificwith
threeunigenesshowingsignificantlyhighexpression.Tissue-specificunigeneswerealsoin-
volvedinproteinsynthesisandsecretionsuchasgenesencodingsilkglandfactor3,forkfead
transcriptionfactorG1,ABCtransporters,andyellowproteins.InthePSG,threeunigenesen-
codingfibroinproteins,fibroin/P25,fibroin-Lchainandfibroin-Hchain,weresignificantly
specificallyexpressed.Twounigenesencodingglycine-tRNAligaseandalanine-tRNAligase
showedsignificantlyhighexpressioninthePSG.Thesegeneproductsprovideglycinesandala-
ninesforefficientsynthesisoffibroin.Otherspecificallyexpressedunigenesencodedtranscrip-
tionregulationfactors,heatshockproteinsandsugartransporters.
Table2. SummaryofexpressionrelatedtosilkproteinsbetweenMSGandPSG.
WBm_unigene length FPKM_MSG FPKM_PSG logFC PValue FDR Bm_CDS Ka/Ks Putativefunction
comp4655_c0_seq1 562 21.07 1580.95 6.035 1.29E-13 9.94E-12 Fibroin-H
comp11395_c0_seq1 1362 1298.87 137786.13 6.534 3.26E-15 3.24E-13 BGIBMGA009393 0.6473 fibroin-L
comp8545_c1_seq1 1353 182.40 31454.66 7.235 2.73E-17 3.76E-15 BGIBMGA001347 0.1938 fibroin/P25
comp14614_c0_seq5 2173 8339.87 6.45 -10.532 4.54E-27 4.99E-24 sericin1
comp14333_c2_seq1 1846 7606.27 6.60 -10.366 1.51E-26 1.03E-23 BGIBMGA011901 1.2050 sericin2
comp17080_c0_seq1 480 5175.15 6.27 -9.877 1.03E-24 4.19E-22 BGIBMGA012002 1.8625 sericin3
comp13620_c0_seq1 3130 65.03 408.04 2.453 3.71E-04 4.75E-03 BGIBMGA006216 alanine-tRNAligase
comp4656_c0_seq1 2624 121.56 986.25 2.825 5.35E-05 9.04E-04 BGIBMGA007637 0.0198 glycine-tRNAligase
comp17114_c0_seq1 1334 530.01 0.43 -10.453 3.77E-25 1.74E-22 BGIBMGA009261 fibrohexamerin
comp5260_c0_seq1 1067 3.18 0.05 -5.883 1.63E-06 3.95E-05 BGIBMGA009261 fibrohexamerin
comp7638_c0_seq1 1509 1118.49 0.75 -10.727 1.15E-26 8.65E-24 BGIBMGA009261 1.6506 fibrohexamerin
comp8883_c0_seq1 954 24544.0 15.54 -10.819 6.09E-28 1.46E-24 BGIBMGA009261 fibrohexamerin
doi:10.1371/journal.pone.0122837.t002
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 6/14
TranscriptomeAnalysisoftheWildSilkworm
Table3. SummaryofSNPsandIndelsforBombyxmandarina.
Length(bp) SNP(num) SNP_density(bp/num) Indel(num) Indel_density(bp/num)
Unigene 25,730,955 32,297 797 361 71,277
ORF 11,584,893 19,783/16,074/3,709* 586 14 827,492
Three_prime_UTR 4,051,613 4,974 815 179 22,635
Five_prime_UTR 1,640,684 2,024 811 54 30,383
*Totalnumber/synonymoussitenumber/nonsynonymousnumber
doi:10.1371/journal.pone.0122837.t003
ReadsmappingandSNPcalling
Toevaluateheterozygosisinthewildsilkwormgenome,cleanedreadsweremappedtoassem-
bledunigenes.Atotalof32,297SNPsand361INDELswereidentifiedinassembledunigenes
withatotallengthof25,730,995bp(Table3)andwithadensityof797bpperSNP;19,783
SNPswerelocatedinORFs,4974in3'-UTRs,and2024in5'-UTRswithadensityperSNPof
586bpforORFs,815bpfor3'UTRs,and811bpfor5'-UTRs.Only14INDELswereidentified
intheORFs,whichwaslessthanthe179in3'-UTRsand54INDELsin5'-UTRs.Moreover,
16,074(81.3%)SNPsinORFinducedsynonymousaminoacidsubstitutions(Table3).
Identificationofputativeorthologsandanalysisofpositiveselection
Areciprocalbest-hitmethodwasusedtoidentifyputativeorthologsbetweenthetwospecies
[21].Atotalof5295pairsofputativeorthologswereidentifiedbycomparingthewildsilkworm
unigenesandthesilkwormCDSsequences.TheKa/KsratiowascalculatedbyKaKs_Calcula-
torsoftwareshowing3855pairswithKaandKsforcalculatingKa/Ksand1440pairswithonly
KaorKs.The3855pairshadmeanvaluesof0.0624forKa,0.4845forKsand0.6409forKa/
Ks.Among3855pairsoforthologs,1049pairshadKs>1andwereconsideredpotentialpara-
logs.Afterremovingpotentialparalogs,2,806pairsbelongingtoorthologswithmeanvaluesof
0.0109forKa,0.0433forKs,and0.7316forKa/Ksratio(S4Table).
Atotalof83(2.96%)pairsoforthologshadaKa/Ks>1showinggenesthatmighthaveex-
periencedstrongpositiveselection(Fig3).AKa/Ksratioof0.5wasconsideredasausefulcut-
offtoidentifygenesunderpositiveselection[27],for317(11.30%)pairsoforthologswitha
Ka/Ksratiobetween0.5and1(Fig3).Thus,400(Ka/Ks>0.5)pairswereconsideredlikelyto
haveexperiencedorbeexperiencingpositiveselection.Comparingwith354candidatedomes-
ticgenesingenomicregionsofselectivesignals(GROSS)fromresequencingdata[1],13genes
werecommonbetweenKa/KsandGROSSanalysis(S5Table).Inaddition,twounigenes(Ka/
Ks=0.68,1.31)encodingcytochromeb-likeproteinswereconsideredtobeunderpositivese-
lection.Thecytochromebgeneshowedastrongsignalofpositiveselectioninthedomesticated
clade[2].
Functionalanalysisforpositiveselection
Among400unigeneswithKa/Ks>0.5consideredaspositivelyselected,168wereannotated
withKOassignments(42%).KEGGpathwayanalysiswasperformedwithKAASfor126
KEGGpathwaysidentifiedwithunigenes(S6Table).Theunigenesweredividedtofivebiologi-
calclassesusingpathwayenrichmentanalysis.Fig4Ashowsdistributionoffiveclassesbetween
positiveselectionandalltheunigenesinpathways.Onlytheclassmetabolismhadasignifi-
cantlyhigherpercentageamongpositivelyselectedgenesthanamongallunigenes(35.71%vs.
23.43%).Themetabolismclasswasapproximatelydividedintosixclassesbetweenpositively
selectedandallunigenesinpathways:aminoacidmetabolism(16.25%vs.18.07%),
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 7/14
TranscriptomeAnalysisoftheWildSilkworm
Fig3.DistributionofKaandKsvalues.Abovetheblackline,orthologouspairswithKa/Ksratio>1;
betweenblackandgraylines,pairswithKa/Ksratio0.5–1.
doi:10.1371/journal.pone.0122837.g003
carbohydratemetabolism(37.50%vs.35.63%),energymetabolism(8.75%vs.8.09%),lipidme-
tabolism(12.50%vs.14.15%),nucleotidemetabolism(12.50%vs.15.86%),andothers(12.50%
vs.8.21%)(Fig4B).Pathwayenrichmentanalysisidentifiedthreeenrichedpathways(p(cid:1)0.01):
glycosaminoglycanbiosynthesis-chondroitinsulfate/dermatansulfate(ko00532),retinoicacid-
induciblegeneI(RIG-I)-likereceptorsignalingpathway(ko04622),andcircadianrhythm-fly
(ko04711)(Table4).
Discussion
Inthisstudy,weperformeddenovotranscriptomesequencingforB.mandarinausingIllumina
HiSeq2000sequencingplatform.Weobtainedmorethan100millioncleanreadsand33,759
assembledunigeneswithN50length1437bpandaveragelength762.20bp.Atotalof883and
425unigeneswereup-regulatedexpressedintheMSGand425inPSG,respectively.Weidenti-
fiedsixunigenesencodingmajorsilkcomponents:threefibroinproteins(fibroin-Hchain,
fibroin-Lchain,andfibroin/P25)andthreesericinproteins(sericin1,sericin2,andsericin3).
TheseunigeneswerespecificallyexpressedintheMSGorPSGwithsimilarexpressionprofiles
oftheirhomologousgenesinthedomesticatedsilkworm[28].Fibroin-H(FPKM=1581;
logFC=6.04),fibroin-L(FPKM=137,786;logFC=6.53),andfibroin/P25(FPKM=31,455;
logFC=7.23)werespecificallyexpressedinthePSG,whereassericin1(FPKM=8340;
logFC=-10.53)sericin2(FPKM=7606;logFC=-10.37)andsericin3(FPKM=5175;
logFC=-9.88)werespecificallyexpressedintheMSG.Animmunosorbentassaywithspecific
antibodyagainsteachfibroinindicatedthatamolarratiooffibroin-Hchain,fibroin-Lchain,
andP25of6:6:1[29].However,theratioinwildsilkwormswas0.05:4.38:1.Weproposethat
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 8/14
TranscriptomeAnalysisoftheWildSilkworm
Fig4.DistributionofKEGGbiologicalcategoriesofpositivelyselectedandallunigenesinpathways.(A)Distributionoffiveclassesofpositively
selectedandallunigenesinpathways;(B)Metabolismshowsdistributionofsixclasses.(a)Distributionofpositivelyselectedunigenesinpathways,(b)
Distributionofallunigenesinpathways.
doi:10.1371/journal.pone.0122837.g004
lowexpressionoftheFibroin-Hgenemightbeduetoconstraintsofsequencetechnologyfor
highGC-contentandsimple-repeat-sequence-containinggenes.Wealsofoundunigenesin-
volvedinproteinbiosynthesisandsecretion.Fibroinhashighproportionsofglycine(44.6%),
alanine(29.4%)andserine(12.1%)[30].Forexample,twounigenesencodingglycine-tRNAli-
gase(logFC=2.82)andalanine-tRNAligase(logFC=2.45)weremorehighlyexpressedinthe
PSGthanintheMSG,suggestingthatthetwogenesparticipateinsericinandfibroinbiosyn-
thesis.Fourunigenesencodingfibrohexamerin-likeproteinswerespecificallyexpressedwith
Table4. KEGGpathwayenrichmentanalysisofpositiveselection.
KEGGpathway Unigene Bm_orthologs Ka/Ks Putativefunction
Glycosaminoglycanbiosynthesis—chondroitinsulfate/ comp12639_c0_seq1 BGIBMGA007765 0.5734 beta-1,3-galactosyltransferase
dermatansulfate(ko00532)
comp1835_c0_seq1 BGIBMGA003101 0.5113 carbohydratesulfotransferase
comp35605_c0_seq1 BGIBMGA012390 0.5474 chondroitinsulfatesynthase
Circadianrhythm—fly(ko04711) comp15415_c0_seq2 BGIBMGA007304 0.9945 double-timeprotein
comp18252_c0_seq1 BGIBMGA000486 0.5633 periodcircadianprotein(PER)
comp6926_c0_seq1 BGIBMGA000498 0.7571 circadianlocomoteroutputcycleskaput
protein(COLCK)
RIG-I-likereceptorsignalingpathway(ko04622) comp11847_c0_seq1 BGIBMGA003954 0.5546 autophagy-relatedproteinAtg12
comp12369_c0_seq1 BGIBMGA000198 0.5607 TANK-bindingkinase
comp6686_c1_seq1 BGIBMGA000679 0.9135 rotamasePin1
doi:10.1371/journal.pone.0122837.t004
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 9/14
TranscriptomeAnalysisoftheWildSilkworm
FRPMvalues(24,544,1118,530and3.18)intheMSG.Thesilkwormgenomehaseightfibro-
hexamerin-likegeneswithBmfhxh4expressedspecificallyintheMSGandinvolvedinprotein
assemblyandsecretion[31].
Thesilkglandisoneofthemostimportanttissuesthathasundergoneartificialselectionfor
moreandbetterqualitycocoons.Manyphenotypeshaveundergonechangesbetweenwildand
domesticatedsilkwormssuchassilkglandsizeandsilkoutput,andcocoonqualityandcolor.
BycomparingwithB.morigenes,weidentified5295pairsoforthologousgenesbyKa/Ksanal-
ysis,with,inamongwhich400genesmighthaveexperiencedortobeexperiencingpositivese-
lectionbyKa/Ksanalysis.Forinstance,threegenesencodingsericin2,sericin3,and
fibrohexamerin-likeexperiencedstrongpositiveselectionwithKa/Ksratiosof1.21forsericin
2,1.86forsericin3and1.65forthefibrohexamerin-likegene.Fibroin-Lwasfoundtobeunder
positiveselection(Ka/Ks=0.65).Inadditiontothesestructuralproteins,byKEGGanalysis,
severalsignalingpathwayswereinvolvedinorganortissuedevelopmentandnutritionalsig-
nals.Aunigeneencodinginsulin(Ka/Ks=1.12)hasexperiencedstrongpositiveselection.
OnlytwoSNPsweredetectedbetweenwildanddomesticatedsilkworms;bothledtononsy-
nonymoussubstitutions(I42-V42andS110-N110).Insulinisinvolvedinmanysignalingpath-
wayssuchastheinsulinsignalingpathway(ko04910),themTORsignalingpathway
(ko04150),andtheFOXOsignalingpathway(ko04068).Inthesesignalingpathways,eight
unigenesencodinganinsulingrowthfactor1receptor(Ka/Ks=0.79),aphosphorylase
kinase(Ka/Ks=0.62),aras-relatedGTPbindingprotein(Ka/Ks=0.53),aG-proteinbeta
subunit(Ka/Ks=0.86),aubiquitin-likeproteinATG12(Ka/Ks=0.55),acaseinkinase1
(Ka/Ks=0.99),aCREB-bindingprotein(Ka/Ks=0.73),andaserine/threonine-proteinphos-
phatasePP1catalyticsubunit(Ka/Ks=0.58),havebeenexperiencingpositiveselection.Thein-
sulinandmTORsignalingpathwaysareevolutionarilyconservedinmosteukaryotesandare
crucialforsystemictransductionofnutritionalsignalsregulatingcellgrowthandmetabolism
[32,33].Thesegenesandpathwayshaveexperiencedpositiveselection,suggestingthatthey
shouldbemoleculartargetsofartificialselectioninthedomesticatedsilkworm.
ThreepathwayswereenrichedbyKEGGpathwayenrichmentanalysis,includingglycos-
aminoglycanbiosynthesis-chondroitinsulfate/dermatansulfate(ko00532)(p=0.006),RIG-I-
likereceptorsignalingpathway(ko04622)(p=0.007),andcircadianrhythm-fly(ko04711)
(p=0.004).ThechondroitinsulfatemetabolicandRIG-I-likereceptorsignalingpathwaysare
relatedtohostimmuneresponse[34,35].Inthechondroitinsulfatemetabolicpathway,three
unigenesencodingabeta-1,3-galactosyltransferase(Ka/Ks=0.57),acarbohydratesulfotrans-
ferase(Ka/Ks=0.51),andachondroitinsulfatesynthase(Ka/Ks=0.55)havebeenexperienc-
ingpositiveselection.Chondroitinsulfateintheperitrophicmembranemightprotectthe
midgutepitheliumfromingestedpathogens[34].IntheRIG-I-likereceptorsignalingpathway,
threeunigenesencodedanautophagy-relatedproteinAtg12-likeprotein(Ka/Ks=0.55),a
TANK-bindingkinase(Ka/Ks=0.56),andarotamasePin1(Ka/Ks=0.91)andareunderposi-
tiveselection.RIG-I-likereceptorsrecognizeRNAvirusesandtriggerarobustinnateimmune
responseagainstRNAvirusinfection[35].Thesignificantenrichedtwopathwaysmightbea
resultofsilkwormdomesticationagainstpathogens.
Thedomesticationprocessalsointroducesdiversityinbehavioralphenotypesbetweenwild
anddomesticatedsilkworms.Forexample,thedomesticatedsilkwormmothhaslostflightabil-
itythroughartificialselection.Aunigene(Ka/Ks=1.14)encodingaleucine-richrepeatflight-
less-interactingprotein2(LRRFIP2)hasundergonestrongpositiveselection.LRRFIPinteracts
withaleucine-richrepeatflightlessproteinthatisamemberofthegelsolinfamilyandwasdis-
coveredasamutationaccompanyingdisorganizedflightmusclemyofibrilsinDrosophilamela-
nogasterthatleadstoflightlessness[36,37].ArtificialselectiontargetingLRRFIPmightbe
involvedwithlossofflightabilityofthedomesticatedsilkworm.Thecircadianrhythm
PLOSONE|DOI:10.1371/journal.pone.0122837 March25,2015 10/14