Rapid transcriptome sequencing of an invasive pest, the brown

22 Pages·2014·1.66 MB·English

Ioannidisetal.BMCGenomics2014,15:738
http://www.biomedcentral.com/1471-2164/15/738
RESEARCH ARTICLE Open Access
Rapid transcriptome sequencing of an invasive
pest, the brown marmorated stink bug
Halyomorpha halys
Panagiotis Ioannidis1,4, Yong Lu2, Nikhil Kumar1, Todd Creasy1,5, Sean Daugherty1, Marcus C Chibucos1,3,
Joshua Orvis1, Amol Shetty1, Sandra Ott1, Melissa Flowers1, Naomi Sengamalay1, Luke J Tallon1, Leslie Pick2
and Julie C Dunning Hotopp1,3*
Abstract
Background: Halyomorphahalys(Stål)(Insecta:Hemiptera;Pentatomidae),commonlyknownastheBrownMarmorated
StinkBug(BMSB),isaninvasivepestofthemid-AtlanticregionoftheUnitedStates,causingeconomicallyimportant
damagetoawiderangeofcrops.NativetoAsia,BMSBwasfirstobservedinAllentown,PA,USA,in1996,andthispest
isnowwell-establishedthroughouttheUSmid-Atlanticregionandbeyond.InadditiontotheseriousthreatBMSB
posestoagriculture,BMSBhasbecomeanuisancetohomeowners,invadinghomegardensandcongregatinginlarge
numbersinhuman-madestructures,includinghomes,tooverwinter.Despiteitssignificanceasanagriculturalpestwith
limitedcontroloptions,only100bpofBMSBsequencedatawasavailableinpublicdatabaseswhenthisprojectbegan.
Results:Transcriptomesequencingwasundertakentoprovideamolecularresourcetotheresearchcommunityto
informthedevelopmentofpestcontrolstrategiesandtoprovidemoleculardataforpopulationgeneticsstudiesof
BMSB.Usingnormalized,strand-specificlibraries,wesequencedpoolsofallBMSBlifestagesontheIlluminaHiSeq.
Trinitywasusedtoassemble200,000putativetranscriptsin>100,000components.Anovelbioinformaticmethodthat
analyzedthestrand-specificityofthedatareducedthisto53,071putativetranscriptsfrom18,573components.By
integratingmultipleotherdatatypes,wenarrowedthisfurtherto13,211representativetranscripts.
Conclusions:Bacterialendosymbiontgeneswereidentifiedinthisdataset,someofwhichhaveacopynumber
consistentwithbeinglateralgenetransfersbetweenendosymbiontgenomesandHemiptera,includingankyrin-repeat
relatedproteins,lysozyme,andmannanase.SuchgenesandendosymbiontsmayprovidenoveltargetsforBMSB-specific
biocontrol.Thisstudydemonstratestheutilityofstrand-specificsequencingingeneratingshotguntranscriptomesand
thatrapidsequencingshotguntranscriptomesispossiblewithouttheneedforextensiveinbreedingtogenerate
homozygouslines.Suchsequencingcanprovidearapidresponsetopestinvasionssimilartothatalreadydescribedfor
diseaseepidemiology.
Keywords:Brownmarmoratedstinkbug,Halyomorphahalys,Transcriptome,Hemiptera,Invasivespecies,Lateralgene
transfer,Horizontalgenetransfer,Mannanase,Lysozyme
*Correspondence:[email protected]
1InstituteforGenomeSciences,UniversityofMarylandSchoolofMedicine,
Baltimore,MD,USA
3DepartmentofMicrobiologyandImmunology,UniversityofMaryland
SchoolofMedicine,Baltimore,MD,USA
Fulllistofauthorinformationisavailableattheendofthearticle
©2014Ioannidisetal.;licenseeBioMedCentralLtd.ThisisanOpenAccessarticledistributedunderthetermsoftheCreative
CommonsAttributionLicense(http://creativecommons.org/licenses/by/4.0),whichpermitsunrestricteduse,distribution,and
reproductioninanymedium,providedtheoriginalworkisproperlycredited.TheCreativeCommonsPublicDomain
Dedicationwaiver(http://creativecommons.org/publicdomain/zero/1.0/)appliestothedatamadeavailableinthisarticle,
unlessotherwisestated.
Ioannidisetal.BMCGenomics2014,15:738 Page2of22
http://www.biomedcentral.com/1471-2164/15/738
Background spread, making it a particularly difficult invasive pest to
Halyomorpha halys (Stål) (Insecta: Hemiptera: Pentato- manage. Over sixty researchers funded by the USDA and
midae), otherwise known as the Brown Marmorated commodity organizations are conducting experiments
Stink Bug (BMSB), is an invasive pest that has ravaged aimed at better understanding the biology and ecology of
farms and distressed homeowners in the mid-Atlantic thispestandtofindmanagementsolutions.Whenthispro-
region of the US in recent years. It is a polyphagous in- ject began, only 100 bp of sequence data was available in
sect pest that causes economically important damage to the public database. Transcriptome sequencing wasunder-
many crops including vegetables, tree fruit, field crops, takentoprovideamolecularresourcetotheresearchcom-
andornamentals[1]. munity for basic and applied research purposes. The data
Native to Asia, BMSB was first observed in North can be used to examine the population genetics of BMSB,
AmericainAllentown,PA,USAin1996[1].BMSBisnow including gaining a better understanding of the original
well-establishedinthemid-AtlanticregionsoftheUSand invasion in Allentown, as well as subsequent invasions
hasspreadto41different statesandDC[2].Thereisalso throughouttheUSandtheworld.Understandingthesein-
evidence for established populations in Canada [3] and vasionswillincreaseourknowledge,preventingfutureinva-
Switzerland [4-6]. It is a pest of tree fruit; grapes; other sions. Our study demonstrates that the transcriptomes of
small fruit; row crops including vegetables, legumes and invasive species can be rapidly sequenced, providing a
cotton; shade trees; ornamentals; and nursery crops [7-9]. resource to the research community without extensive
BMSB crop damage reached economically significant breedingtocreatehomozygouslines(e.g.generating8gen-
levels during the 2010 growing season. Growers consid- erationsofanisofemaleline).Suchsequencingcanprovide
ered BMSB to be the single most important pest in the arapidresponsetopestinvasionssimilartothatalreadyde-
mid-Atlantic region, leading to ‘desperate measures’ in- scribedforepidemiologyofinfectiousdisease(e.g.[15])and
cluding increased use of broad-spectrum pesticides to jumpstartmolecularbiologyandgenetic-basedstudies.
controlBMSB,detractingfromsustainableIntegratedPest
Management practices. In keeping with this, the EPA Results
granted emergency approval for the use of several pesti- Sequencingresults
cides, including one approved for use in organic farming, Whole transcriptome sequencing was undertaken in order
to prevent further economic loss in the mid-Atlantic re- to obtain gene sequences and jumpstart molecular biology
gion [10]. In addition to being an agricultural pest, many studies focused on BMSB. RNA was collected for sequen-
mid-Atlantic homeowners are troubled by BMSB. Unlike cing from 10 different stages and conditions including: (a)
nativestinkbugs,BMSBaggregatesinhuman-madestruc- eggs, (b) 1st instar nymphs, (c) 2nd instar nymphs, (d) 3rd
tures, including houses, to overwinter [11]. One home instarnymphs,(e)4thinstarnymphs,(f)5thinstarnymphs,
owner reported removing 26,000 BMSB from his resi- (g) an active adult male, (h) an active adult female, (i) an
denceinthefirsthalfof2011[12]. adult male in diapause, and (j) an adult female in diapause
BMSB has a long history of hitchhiking to new areas. (Figure 1). Two pools were created with equimolar
The successful introduction in Allentown, PA, USA, was amountsofeachRNAsamplewithpoolAcontainingallof
preceded by the successful interception of BMSB by the thepre-adultstagesandpoolBcontainingtheadults.
USDA preventing such invasions including an intercep- Each pool was used to generate a strand-specific library
tion of BMSB on a 1983 flight from Japan, an intercep- that was subsequently normalized with DSN (Figure 2).
tion in baggage from Korea in 1984, and eight further Strand-specific sequencing allows for assignment of the
interceptions from China, Korea, and Japan from 1989– strand of transcription while normalization decreases the
1998 [8]. BMSB has been unintentionally shipped from sequencing of the most abundant transcripts allowing for
Japan to New Zealand in a used vehicle [13]. In 2005, sequencing of more rare variants. Sequencing of these li-
over a dozen BMSB were recovered from a storage unit braries on the Illumina HiSeq2000 resulted in 196,233,912
in Vallejo, Solano County, CA, that was rented by an in- readstotaling>19billionbasesforPoolAand170,455,294
dividual who had relocated from Pennsylvania [14]. Eco- reads totaling >17 billion bases for pool B. All of the reads
logical niche modelling and climate data predict that were pooled and assembled with Trinity [16], which uses
areas between 30–50 degrees latitude are at risk of inva- Inchworm to produce contigs via greedy k-mer extension,
sion including Northern Europe, north-eastern North Chrysalis to produce components containing related con-
America, the northern portions of the North American tigs, and Butterfly to generate transcripts using the various
west coast, southern Australia, the North Island of New de Bruijn graph paths in the components. We will refer to
Zealand as well as Angola in Africa and Uruguay in theassemblageofsequencesgeneratedbyButterflyasputa-
South America [8]. tive transcripts, given that some of the sequences that
As summarized above, BMSB are significant agricultural Butterflyassemblesarenotactuallytranscripts,asdescribed
pests with limited treatment options and an ability to below.
Ioannidisetal.BMCGenomics2014,15:738 Page3of22
http://www.biomedcentral.com/1471-2164/15/738
mappingsthatallowedformultiplematches,whichwillbe
referred to as the SSLR for strand-specific log ratio. We
wereabletoshowthatoverallthelibraryisstrand-specific
with the plus strand preferred to the minus strand 3.40:1
(Figure4A).WhiletheSSLRdemonstratestheplusstrand
was preferred, there was a significant amount of leakage
that results in two putative transcripts that are reverse
complements to one another. Leakage transcripts were
identified as those that (a) matched another putative
transcriptinthealternateorientationusingBLASTN,and
(b) had a positive SSLR as described above. Using this
approach,wediscarded32,308putativetranscriptsorigin-
ating from 24,507 components. Of those putative tran-
scripts, there were 20,795 that had an ORF predicted on
only one of the strands, and thus in only one direction.
Notunexpectedly,20,350hadtheirORFspredictedonthe
minus strand compared to only 445 that had their ORFs
in the plus strand. Overall, those with a SSLR<−5 had
>10-fold increase in ORFs on the plus strand relative to
Figure1LifestagesofBMSB.ThelifestagesofBMSBareshown the minus strand and those with a SSLR >5 had >10-fold
startingwitheggsfollowedby1stinstarnymphs,2ndinstarnymphs,3rd
increase in ORFs on the minus strand relative to the plus
instarnymphs,4thinstarnymphs,5thinstarnymphs,andanadultina
strand (Additional file 1: Figure S1). This indicates that
counter-clockwisespiraloutwardsandfromlargesttosmallest.Thebar
inthelowleftrepresents1cm. usingaSSLRcancorrectlypredicttherightstrand.
Upon mapping the distribution of putative transcripts
by strand specificity, one observes a bimodal distribution.
The resulting assembly contained 194,729 putative The bimodal distribution becomes more symmetrical
transcripts with 50,599 that were >1,000 bp and 20,076 upon applying a coverage threshold. In particular, there
having a top BLASTX match (e≤10−10 and over >80% arenumerousputativetranscriptswithastrandspecificity
of the reference) to 8,354 unique Uniref90 protein se- slightlylessthanzerothatdisappearedcompletelyasmore
quences. For comparison, Trinity produced 196,000 pu- coverage was required to support the underlying putative
tative transcripts with 14,522 that were >1,000 bp with transcript (Figure 4, Additional file 2: Figure S2). One ex-
4,323 having a top BLASTX match (e≤10−10) to 2,880 planation for this would be that these putative transcripts
unique Uniref90 proteinsequences over >80% ofthe ref- arise from contaminating genomic DNA that lacks strand
erence protein length for an outbred Bemisia tabaci specificity.Thiswouldmanifestasshortcontigssupported
white fly population (Insecta: Hemiptera: Aleyrodidae) by few reads and little strand specificity. Consistent with
[17]. The longest assembled putative transcript from genomic DNA being present in the samples, applying a
BMSB is 30 kbp in length; most components contained coverage cutoff of only two removed >10,000 putative
a single putative transcript; and the distribution of the % transcripts that are <300 bp. In other projects, we have
GCpeaked around 37%(Figure 3). observed contaminating E. coli DNA that we suspect ar-
rives in the sample through molecular biology reagents
Strand-specificityandcoverage that may allow for an examination of our ability todetect
One problem with assigning a function to transcripts is genomic DNA contamination. In BMSB, bacterial DNA
determining which open reading frame is the correct may also come from the microbiome, for example from
one when there are open reading frames on both strands an endosymbiont. However, regardless of the source, E.
that lack evidence of function. This can be particularly coli DNA could serve as a proxy for identifying genomic
true for organisms where good reference annotation is DNA in this sample. A coverage cutoff of two removed
lacking, like insects. Sequencing a strand-specific library 3,501 putative transcripts that matched E. coli. If we plot
should assist in identifying the transcribed strand the mean length, number of putative transcripts of size
(Figure 2), reducing the number of possible ORFs by <300bp,andnumberofreadswithamatchtoE.coli,the
half. However, strand-specific libraries can have leakage, greatestimprovementsareobtainedwithacoveragecutoff
where experimentally the wrong strand is sequenced, of 20 (Figure 5). Three putative transcripts still persist
thatneedstobeaddressed. with homology to E. coli. One of these (comp3244_c38_-
To measure strand specificity, we calculated the log ra- seq1) has little strand specificity (SSLR=−0.45) and
tio of sense reads to antisense reads using Bowtie2 matchesPhiX,whichisaddedtoIlluminasequencingruns
Ioannidisetal.BMCGenomics2014,15:738 Page4of22
http://www.biomedcentral.com/1471-2164/15/738
Random primed
first strand cDNA synthesis Second strand synthesis
& removal of unincorporated nucleotides with dUTP mix & RNase H
Sense cDNA
Sense mRNA with dUTP
Sense mRNA 5’ 3’ 5’ 3’
555’’’ 3’
3’ 5’ 3’ 5’
End repair, addition of A-tail,
adapter ligation with index
Digestion of 2nd strand
PCR amplification DNA with UNG 5’ Sense cDNA 3’
Sense DNA with dUTP
5’ 3’ orange adaptor contains
the index sequence
3’ 5’
3’ 5’ 3’ 5’
DSN normalization &
an additional PCR
amplification
3’ 3’ 3’ 3’
Sequence read 1 Sequence read 2
A
A A will be antisense A Flipping and will be sense
N First strand N N second strand
e D sequencing e D Indexing e D sequencing
s s s
n n n
e e e
S S S
Index read
5’ 5’ 5’ 5’
Flow Cell Flow Cell Flow Cell Flow Cell
Figure2Strand-specificsequencingwithDSNnormalization.Instrand-specificsequencing,thefinalsequencingreadscanbeassignedtoa
particularstrandofDNA.Thisisaccomplishedbyasecond-strandsynthesiswithdUTP,ligationofadaptorstodouble-strandedDNA,and
degradationofthesecondstrandwithUNG.Theresultisafragmenttobesequencedthatisdifferentiallylabeledonthetwoendsinamanner
thatdictateshowtheyareloadedonthesequencerandthustheorderinwhichtheendsaresequenced.
Figure3Featuresofthetrinityassembly.(A)Thesizedistributionoftheputativetranscriptsonalogarithmicscalerevealsthatthevast
majorityofassembledtranscriptsare<1kbpbuttranscripts>25kbpcanbeassembled.(B)Forthevastmajorityofgenesthereisonlyasingle
transcriptassembled.(C)The%GCcontentoftheputativetranscriptspeaksaround35%.
Ioannidisetal.BMCGenomics2014,15:738 Page5of22
http://www.biomedcentral.com/1471-2164/15/738
and is not always completely removed through standard
25000
AA >=1 B >=20 filtering. The remaining two, comp1098_c0_seq1 and
20000 comp1100_c0_seq1, are strand-specific, unique, and have
n:p = 3.85 n:p = 2.26
BLASTNmatchestothe16Sand23SrRNA,respectively,
15000
with best matches to the genome of an endosymbiont of
10000 the BMSB (GenBank: AP012554). This highlights that
combining the results of the SSLR, a coverage threshold,
5000
and the BLASTN search significantly improves the tran-
0 scriptome assembly and aids in identifying transcripts
15 10 −5 0 5 10 1515 10 −5 0 5 10 15 from both the host and its endosymbiont(s). The SSLR
− − − −
would be an easy and useful addition to implement into
Figure4Assessmentofthestrand-specificityoftheRNA-Seq
library.Thestrandspecificityofthesequencinglibrarywasbased Trinityandotherassemblersthatwouldprovidethegreat-
oncomparingthenumberoffirst-in-pairreadsmappedtotheplus estgains.
strand,tothecorrespondingnumberofsecond-in-pairreads.Ina
strand-specificlibrary,thestrand-specificlogratio(SSLR)shouldbe
Filteringtoasinglegenepercomponent
awayfromzero.Intheprotocolusedinthisstudy,wherethe
WhileTrinity can generate almost200,000putativetran-
second-in-pairreadcomesfromthesensestrand,anegativeSSLR
isexpected.(A)ThedistributionofSSLRsforputativetranscripts scripts in >100,000 components, there are thought to be
havingatleastoneproperreadpairmappingtothemdemonstrates only 20,000-40,000 genes in Metazoan genomes. Apply-
apreferenceforthosewithanegativeSSLR.(B)Whentheplotis ing the thresholds above reduces this to 53,071 putative
limitedtothosewithtwentyreadpairsmappingtothem,the
transcripts from 18,573 components. However, analyses
bimodaldistributionbecomesmoredistinct,likelyfromtheremoval
related to function, particularly those that result in stat-
ofcontigsarisingfromcontaminatinggenomicDNA.
istical analyses, need the data to be reduced to a repre-
sentative number of genes, thereby removing splice
variants. Using an integrated approach, we were able to
12000 40000
35000
10000
30000
8000
25000
6000 20000
15000
4000
10000
2000
5000
0 0
0 10 20 30 40 50
Coverage
Figure5Transcriptpropertiesrelativetotranscriptcoverage.Themeanlength(orange),numberofputativetranscripts<300bp(gray),and
homologytoanE.colisequence(yellow)wereplottedontheleftaxiswhilethenumberofputativetranscripts(blue)wereplottedontheright
axisafterapplyingacoveragethresholdtothetranscript(x-axis).Acoveragecutoffof20maximizesthemeanlengthandminimizesthenumber
oftranscriptswithhomologytoE.coli.
Ioannidisetal.BMCGenomics2014,15:738 Page6of22
http://www.biomedcentral.com/1471-2164/15/738
score each putative transcript based on two types of evi- substitution mutation. These positions were identified
dence including similarity to known sequences and the in 11,462 different putative transcripts, which is 86.8%
intrinsic features of each contig. We calculated a score of the total number of putative transcripts. The genetic
for each of the remaining 53,071 putative transcripts of heterogeneity observed may be significantly less than
theTrinity assemblybased on(a) similarityto knownse- that found in the native Asian population, since the
quences, (b) features of any predicted ORFs, and (c) se- founder population in the US was likely small, repre-
quencing coverage along the putative transcript. The senting a significant bottleneck, and an additional
exactmetricsusedforcalculatingthescorearedescribed bottleneck occurred when founding the laboratory
in Table 1 and Methods. This filtering resulted in 13,211 population. However, these polymorphisms may be a
putative transcripts representing a 14-fold reduction in useful resource for population studies examining the
the number of putative transcripts and a 9-fold reduc- spreadofBMSBintheUS.
tion in the number of components (Table 2). This num-
ber of putative transcripts is consistent with sequenced Differentialexpressionbetweenadultsandnon-adults
hemipteran genomes such as the 34,604 predicted genes In order to compare the expression of juveniles and
inA.pisum[18]. adults, the normalized log Ratio of the read count was
2
calculated between the two pools (Additional file 3:
COGClassification Table S1). A total of 22,707 putative transcripts had
The NCBI Cluster of Orthologous Groups (COG) data- more than a 2-fold change in ratio, and 6,590 had more
base was used to classify the predicted proteins in these than a 16-fold change in ratio. The most highly
13,211 transcripts. Only 5,315 functions were assigned expressed genesin the adult pool B that have functional
using NCBI COGs. Of those, almost a quarter (1,277) are annotation were vitellogenin genes, which is expected
functions that are not well categorized, namely general given that they encode the yolk protein produced in
functional prediction and unknown function (Figure 6). adult females. Similar levels of differential expression
The remaining categorizations are almost evenly divided were found in a large number of hypothetical proteins
betweentheotherthreecategoriesof(a)informationstor- and a protein with homology to an ankyrin protein
age and processing, (b) cellular processes and signaling, found in bacterial endosymbiont Wolbachia strains. In
and (c) metabolism. The top five sub-categories were (a) the juvenile stages, the most abundantly expressed
general function prediction only (22.2%); (b) replication, genes were the larval cuticle proteins, which is also
recombination and repair (8.2%); (c) posttranslational expected. This, combined with the qRT-PCR validation
modification, proteinturnover, and chaperones(7.9%);(d) of specific transcripts of interest described below, dem-
transcription(7.5%);and(e)aminoacidtransportandme- onstrates that DSN normalized libraries can be used for
tabolism(7.1%).Thethreecategoriescomprising(a)extra- a gene expression analysis that successfully identifies
cellular structures, (b) nuclear structure, and (c) cell genes to be targeted for further characterization with
motility were the least abundant COGs represented, with other methods. The dynamic range of the expression
0, 2, and 4 matches (0.0%, <0.1%, and 0.1%), respectively. changes is expected to be muted by the normalization
This is very similar to the distribution of COGs in other procedure, and this is demonstrated for specific tran-
transcriptomesofeukaryotes,includinginsects[19-23]. scriptsbelow.
Juvenile BMSB are more susceptible to insecticides
Polymorphisms than adults ([24,25], and G. Dively, personal communi-
The laboratory colony used to generate RNA for tran- cation). Therefore, we sought to identify transcripts
scriptome sequencing was established recently from with functions related to pesticide detoxification that
field-caught BMSB. While there was at least one bottle- wereup-regulatedinBMSB adults.Therewere25genes
neck in founding the colony, the assembled sequences encoding proteins related to cytochrome P450 and
were expected to have heterogeneity that reflects some glutathione S-transferase that were more abundant in
of the genetic heterogeneity of the US BMSB popula- adults (Table 3). There were nine such proteins that
tion. Of the 26.6 million positions in the 13,211 tran- were more abundant in the pre-adult RNA pool
scripts that were examined with MPILEUP, 23.7 million (Table3).
positions were supported by >20 reads when all reads
were aligned with Bowtie2. Of those positions, 212,422 Taxonomicprofile
had a substitution mutation that differed from the con- The putative transcripts were searched against NR using
sensus base call that was supported by >5% of the BLASTX and a lowest common ancestor (LCA) assign-
underlyingreads,whichrepresents~1%ofthepositions. ment [26] was made by aggregating all of the matches
If more stringent requirements are applied, 109,920 of with the best e-value. Prior to filtering the putative tran-
these positions had >20% of the underlying reads with a scriptstoremoveoneswesuspectedarosefromgenomic
Ioannidisetal.BMCGenomics2014,15:738 Page7of22
http://www.biomedcentral.com/1471-2164/15/738
Table1Variablesandtheirweightsusedtofilterputativetranscripts
Attribute Weight Variable Reason
A)Basedontheputativetranscript
sequence
1.Whatproportionofthedatabase 10 Proportioncovered Rewardtheputativetranscriptbasedontheproportionofthe
proteiniscoveredinthefirst databaseproteinthatiscoveredinthebestBLASTXhit
Uniref100hit?
2.Whatproportionoftheputative 8 Proportioncovered Rewardtheputativetranscriptbasedontheproportionofthe
transcriptiscoveredbythefirst queryputativetranscriptthatiscoveredinthebestBLASTXhit
Uniref100hit?
3.Whatisthelengthcoveredonthe 7 Databasehitlength/longest Rewardbasedontheabsolutedatabaseproteinlengthcoveredin
databaseproteininthefirstUniref100 databasehitlength thebestBLASTXhit,comparedtothelongesthitlengthinthe
hit? component
4.Whatisthelengthcoveredonthe 5 Putativetranscripthitlength/ Rewardbasedontheabsolutequeryputativetranscriptlength
putativetranscriptinthefirst longestputativetranscripthit coveredinthebestBLASTXhit,comparedtothelongesthit
Uniref100hit? length lengthinthecomponent
5.IsthestrandoftheUniref100match, 4 Matchstrand* RewardmatchesinplusstrandifSSLR<0ormatchesinminus
theexpectedone(basedonSSLR)? (−SSLR/max|SSLR|) strandifSSLR>0.Incontrast,penalizematchesinplusstrandif
SSLR>0ormatchesinminusstrandifSSLR<0
6.Whatproportionofthedatabase 9 Proportioncovered SameasthecorrespondingmetricforUniref100
proteiniscoveredinthefirstNRhit?
7.Whatproportionoftheputative 7 Proportioncovered SameasthecorrespondingmetricforUniref100
transcriptiscoveredbythebestNR
hit?
8.Whatistherelativelengthcovered 6 Databasehitlength/longest SameasthecorrespondingmetricforUniref100
onthedatabaseproteininthefirst databasehitlength
NRhit?
9.Whatistherelativelengthcovered 4 Putativetranscripthitlength/ SameasthecorrespondingmetricforUniref100
ontheTrinityputativetranscriptin longestputativetranscripthit
thefirstNRhit? length
10.IsthestrandoftheNRmatch,the 3 Matchstrand* SameasthecorrespondingmetricforUniref100
expectedone(basedonSSLR)? (−SSLR/max|SSLR|)
11.IstheSSLRnegative(i.e.the 7 -SSLR/max|SSLR| Rewardputativetranscriptswiththenormal,negativeSSLR
expected)?
12.Howlongistheputativetranscript 7 Putativetranscriptlength/ Rewardlongerputativetranscripts
comparedtothelongestinthe longestputativetranscript
component? length
B)BasedontheORFs
1.IsthebestmatchforeachORFthe 10 (1-Numberofbestmatches)/ PenalizeputativetranscriptshavingORFsthathavedifferenthits.
same? numberofbestmatches
2.ArethereORFsinbothstrandswith 10 -NumberofORFsinstrand"A"/ MaximumpenaltyifbothORFshaveaNRhit
bothhavinganNRhit? numberofORFsinstrand"B"
3.ArethereORFsinbothstrandswith 8 -NumberofORFsinstrand"A"/ IntermediatepenaltyifonlyoneoftheORFshasaNRhit
onlyonehavinganNRhit? numberofORFsinstrand"B"
4.ArethereORFsinbothstrandswith 3 -NumberofORFsinstrand"A"/ SmallpenaltyifnoneoftheORFshaveaNRhit
noneofthetwohavinganNRhit? numberofORFsinstrand"B"
5.HowmanyORFsarecalled? 8 (1-numberofORFs)/number Penalizeputativetranscriptshaving>1ORFs
ofORFs
6.AretheORFsfoundonlyinthe 8 ORFstrand* RewardputativetranscriptshavingORFscalledinonlythe
expectedstrand(SSLR)? (−SSLR/max|SSLR|) expectedstrand
C)Sequencingcoveragedips
1.Howmanysequencingcoverage 10 -Numberofdips/max Penalizeputativetranscriptswithsequencingcoveragedips
dips? numberofdipsinthe
component
Ioannidisetal.BMCGenomics2014,15:738 Page8of22
http://www.biomedcentral.com/1471-2164/15/738
Table2Summaryofassemblyandannotation
Characteristic Beforefiltering Afterfiltering
Numberofreads(bothpools) 366,689,206 N/A
Numberofputativetranscripts(bothpools) 194,729 13,211
Averagetranscriptlength(bp) 1,005 2,026
Standarddeviation(bp) 1,474 1,592
Mediantranscriptlength(bp) 439 1,649
Maximumtranscriptlength(bp) 27,655 24,046
Transcripts>1000bp 50,599 9,657
TranscriptswithaUniref100hit(e-value<1e-10) 80,536 11,513
TranscriptsmatchinguniqueUniref100proteins 37,160 9,993
TranscriptswithaNRhit(e-value<1e-10) 80,262 11,497
TranscriptsmatchinguniqueNRproteins 37,346 10,007
NumberofTrinitycomponents 123,175 13,211
NumberofORFs 89,684 13,210
NumberofORFs>450bp 61,569 11,141
NumberofORFswithafunctionassigned 57,197 9,811
DNA as described above, only 78,489 (89%) putative filtering, 11,359 (96%) putative transcripts had a
transcripts had aeukaryotic LCAof atotal of88,630 pu- eukaryotic LCA of a total of 11,854 putative transcripts
tative transcripts with BLAST results (Additional file 4: with a BLAST result (Additional file 6: Figure S5). Of
Figure S3). Of the 7,111 (8%) putative transcripts with a the 120 (1%) putative transcripts with a bacterial LCA,
bacterial LCA, 6,366 had an LCA to γ-Proteobacteria only 18 had an LCA to γ-Proteobacteria, and E. coli was
with Escherichia coli as the most common species level no longer the most common species level assignment
assignment (Additional file 5: Figure S4). Following (Additional file 7: Figure S6). This suggests that the
1400
1200
1000
Information storage and processing
800
Cellular processes and signalling
600 Metabolism
Poorly characterized
400
200
0
Translation, ribosomal structure and biogeRnNeAs ipsrocessing and modification TrRaensplcirciapttiioonn, recombination and rCehrpoairCmeallt icn ysctlre ucctounrtre oal, ncd elld ydinvaisimiocns, chromosome partitioning Nuclear structureDefense mechaniSisgmnsal transduction Cemlle cwhaalln/ismemsmbrane/envelope biogenesis Cell motility CytoskeleItnotnracEelxltrulaacr eltlrualaffir csktirnugc, tsuercPersoesttitoran, naslnatdi voensail culmaor dtirfiacnastipoonr,t protein turnover, chapEernoerngesy production andC acrobnovherysdiroante transport and metAabmiolnios amcid transport and metNaucbloeliostimde transport and mCeoteanbzolyismem transport and metabolisLimpid transport andI nmoeSrtegacabonionlci disaormny trmaentsapboorlit taesn db iomseytnatbhoelisiss, mtransport and catabolisGmeneral function prediction onlyFunction unknown
Figure6Transcriptfunctionalcategories.TheNCBIClusterofOrthologousGroups(COG)databasewasusedtoclassifythepredictedproteins
inthe13,211representativetranscripts.AssignmentofCOGcategoriesshowedthatalargenumberofORFsbelongedtocategoriesofproteins
whosefunctionsarepoorlycharacterized,namelythosethathavegeneralfunctionpredictiononlyandthosewithunknownfunction.
Ioannidisetal.BMCGenomics2014,15:738 Page9of22
http://www.biomedcentral.com/1471-2164/15/738
SSLR-based filtering described in the prior section was endosymbionts in this population of BMSB. This is con-
effective at removing contaminating genomic DNA in sistentwiththerecentidentificationofPantoeaagglomer-
thelibrary. ans as a BMSB endosymbiont [28]. The highest number
Intriguedbythepresenceofatranscriptionallyregulated of bacterial LCAs (53 putative transcripts or 44%) were
gene of Wolbachia ancestry in the BSMB transcriptome, from Amoebophilus asiaticus in the Bacteroidetes, which
we sought to investigate the ancestry of these bacterial isanendosymbiontoffree-livingamoebae[29].A.asiaticus
transcripts further. The remaining bacterial putative tran- is related to insect endosymbionts including Blattabacter-
scriptshavehomologytonumerousendosymbionts.LCAs iumspp.incockroaches[30],Sulciamuelleriincicada[31],
for the endosymbiont of the brown-winged green stink and Cardinium hertigii in parasitoid wasps [32]. While the
bug Plautia stali [27] and Pantoea spp. were detected, BLASTXreturnedanLCAofAmoebophilus,theBLASTN
suggesting relatives of these bacteria may be present as searches did not return any significant matches for any of
Table3Differentiallyexpresseddetoxificationgenes
Name Functionalannotation Stageup-regulated
Comp25785_c0_seq1 CytochromeP450gene Adults
Comp26191_c0_seq1 GlutathioneS-transferase Adults
Comp15122_c0_seq1 CytochromeP450gene Adults
Comp16607_c1_seq2 CytochromeP450gene Adults
Comp15049_c0_seq1 CytochromeP450gene Adults
Comp15912_c1_seq2 CytochromeP450gene Adults
Comp20672_c0_seq4 CytochromeP450gene Adults
Comp13685_c1_seq2 CytochromeP450gene Adults
Comp40610_c0_seq2 CytochromeP450gene Adults
Comp8954_c0_seq1 GlutathioneS-transferase Adults
Comp20241_c0_seq1 CytochromeP450gene Adults
Comp18070_c0_seq2 CytochromeP450(2ORFs) Adults
Comp11443_c0_seq1 CytochromeP450gene Adults
Comp18881_c0_seq1 CytochromeP450gene Adults
Comp7095_c0_seq2 CytochromeP450gene Adults
Comp14891_c0_seq2 CytochromeP450gene Adults
Comp8170_c0_seq1 CytochromeP450gene Adults
Comp4236_c0_seq1 ProbablecytochromeP450gene Adults
Comp2339_c0_seq1 CytochromeP450gene Adults
Comp6322_c0_seq2 CytochromeP450gene Adults
Comp17381_c0_seq6 CytochromeP450gene Adults
Comp23582_c1_seq1 CytochromeP450gene Adults
Comp4238_c0_seq1 CytochromeP450gene Adults
Comp3991_c0_seq3 Glutathioneperoxidase(3ORFs) Adults
Comp8540_c0_seq1 ProbablecytochromeP450gene Adults
Comp6146_c0_seq1 GlutathioneS-transferase Pre-adults
Comp11026_c2_seq1 CytochromeP450gene Pre-adults
Comp21713_c0_seq1 CytochromeP450gene Pre-adults
Comp3892_c0_seq1 CytochromeP450gene Pre-adults
Comp18921_c0_seq3 Catalase Pre-adults
Comp25932_c0_seq1 CytochromeP450gene Pre-adults
Comp10873_c0_seq1 ProbablecytochromeP450gene Pre-adults
Comp12303_c0_seq1 CytochromeP450gene Pre-adults
Comp8344_c0_seq1 ProbablecytochromeP450gene Pre-adults
Ioannidisetal.BMCGenomics2014,15:738 Page10of22
http://www.biomedcentral.com/1471-2164/15/738
these, indicating a high level of divergence for the nucleo- amino acids at the N-terminus have little homology for
tide sequences. Most of the proteins with a bacterial LCA comp549_c15_seq3 and comp2753_c5_seq1, respect-
had homology to ankyrin repeat containing proteins found ively. No changes in sequencing coverage around these
indiverseendosymbiontandinvertebrategenomes,includ- regions were identified in either transcript that would
ing all of the proteins with an Amoebophilus LCA. Add- suggest a misassembly. The N-terminus of the ORF
itionaltranscriptswithbacterialLCAsencodeproteinswith encoded in comp549_c15_seq3 contains a JNK_SAPK-
homologytomannanases,amylase,C4-dicarboxylatetrans- associatedprotein-1domain.
porter,andthreedifferenthypotheticalproteins. Examination of the phylogenetic relationship of
comp549_c15_seq3 and comp2753_c5_seq1 with hom-
Ankyrinrepeatcontainingproteinswithsignificant ologous sequences in GenBank (<e-15), reveals that
bacterialhomology while they have homology to Wolbachia endosymbionts
Therewere29LCAstoα-Proteobacteria,allofwhichwere (α-Proteobacteria), collectively they also have homology
in the Rickettsiales, including 27 that were to Wolbachia to coding sequences (CDSs) in Rickettsiella grylii (γ-
endosymbionts.SinceWolbachiaendosymbiontsarefound Proteobacteria) and Diplorickettsia massiliensis (α-Pro-
in many insects, we tested for the presence of Wolbachia teobacteria) (Additional files 8 and 9: Figures S7-S8).
using a robust set of 40 primers that collectively amplify The phylogenetic diversity of the bacterial taxa suggests
genes from a diverse array of Wolbachia endosymbionts as that there have been lateral gene transfers involving di-
wellasothermicroorganismsinthefamilyAnaplasmataceae verse arthropod-associated bacterial lineages and the
[33].WedidnotdetectanyamplificationproductsforWol- precise donors and recipients of these bacteria-bacteria
bachia endosymbionts using BMSB genomic DNA, while lateral gene transfers are not clear. While we suspect
the primers amplified products using control DNA. This that this may also be a lateral gene transfer involving
suggests that this BMSB colony is not colonized by a BMSB, we cannot rule out the presence of these genes
Wolbachiaendosymbiont. arisingfromanendosymbiontofBMSB.ByqPCR,these
Threeputativetranscriptswereidentifiedthatwerediffer- geneswereatthesamerelativeabundanceasBMSBnu-
entiallyexpressedandhadhomologytoWolbachiaankyrin clear genes, but we were not successful at placing these
proteinsincludingcomp549_c15_seq1,comp2753_c5_seq1, genes in the BMSB genome using inverse PCR. The on-
and comp18511_c0_seq1, which had SSLRs of −4.8, −5.4, going BMSB genome sequencing project will likely clar-
and −6.2, respectively. Differential expression analysis of ify this issue. Should this prove to be a lateral gene
thetranscriptomedatareveals28-foldand776-foldoverex- transfer, it is not clear whether the gene moved from
pression of comp549_c15_seq1 and comp2753_c5_seq1 in bacteriaintoinsects,or viceversa.
the adult pool while comp18511_c0_seq1 was 104-fold The third transcript (comp18511_c0_seq1) contains an
overexpressedinthejuvenilepool.Theseresultswere con- ORFwithfulllengthmatchestoWolbachiaproteinsthat
firmed by qRT-PCR using the original RNA samples, prior contain ANK repeats and a PRANC domain. PRANC
to pooling. For qRT-PCR we examined comp549_c0_seq3, domains were at one time described only in Pox viruses,
instead of comp549_c0_seq1, as it is slightly longer and where PRANC-containing proteins inhibit NF-κB activa-
contains part of the latrotoxin domain discussed below. tion by TNF-α, thereby altering the innate immune re-
These qRT-PCR experiments demonstrate a specific sponse of the mammalian host. More recently PRANC
20,000-fold and 17,000-fold increase in expression of domainswereidentifiedinthehumanpathogenOrientia
comp549_c15_seq3 and comp2753_c5_seq1 in active tsutsugamushi (Bacteria: Rickettsiales), multiple Wolba-
adultfemaleswhencomparedtoallotherstages,andaspe- chia endosymbionts (Bacteria: Rickettsiales), Wolbachia
cific170-foldincreaseinexpressionofcomp18511_c0_seq1 phage,andNasoniavitripennis(Eukaryota:Hymenoptera).
in5thinstarnymphs. Previous phylogenetic analysis suggests that the Nasonia
While all three putative transcripts have homology to lineage acquired one or more of these proteins from a
knownankyrin-repeatcontainingproteins,comp2753_c5_- Wolbachia endosymbiont via lateral gene transfer [34].
seq1 and comp549_c15_seq3 do not contain ankyrin ExaminationofthephylogeneticrelationshipoftheBMSB
(ANK)repeats.Bothhavehomologytoalargemultipledo- sequence with homologous sequences in GenBank (<e-
mainproteininWolbachiaendosymbiontsofCulexpipiens 15), reveals that while it has homology to Wolbachia en-
(e.g., WP_007302981.1) that has a latrotoxin domain at its dosymbionts (α-Proteobacteria), it is distinct from the
C-terminus and numerous ankyrin repeats in the N- Wolbachia proteins. It is not clear if the BMSB protein
terminal half. The homology for comp549_c15_seq3 spans represents an ANK-PRANC protein in another bacterial
a portion of the latrotoxin domain and the adjacent ANK lineage or arthropod lineage. As with the genes described
repeats, while the homology for comp2753_c5_seq3 spans above,qPCRresultsrevealedthesamerelativeabundance
only part of this region. This region of homology is at asBMSBnucleargenesandisconsistentwiththesegenes
the C-terminus of BMSB ORFs while the 560 and 223 being in the BMSB nuclear genome. However, in both