Phylogenomic networks reveal limited phylogenetic range of lateral gene transfer by transduction

12 Pages·2017·1.87 MB·English

by  Ovidiu Popa

TheISMEJournal(2017)11,543–554
OPEN
©2017InternationalSocietyforMicrobialEcology Allrightsreserved 1751-7362/17
www.nature.com/ismej
ORIGINAL ARTICLE
Phylogenomic networks reveal limited phylogenetic
range of lateral gene transfer by transduction
Ovidiu Popa1, Giddy Landan and Tal Dagan
Institute of General Microbiology, Christian-Albrechts University of Kiel, Kiel, Germany
Bacteriophages are recognized DNA vectors and transduction is considered as a common
mechanism of lateral gene transfer (LGT) during microbial evolution. Anecdotal events of phage-
mediatedgenetransferwerestudiedextensively,however,acoherentevolutionaryviewpointofLGT
by transduction, its extent and characteristics, is still lacking. Here we report a large-scale
evolutionary reconstruction of transduction events in 3982 genomes. We inferred 17158 recent
transductioneventslinkingdonors,phagesandrecipientsintoaphylogenomictransductionnetwork
view.WefindthatLGTbytransductionismostlyrestrictedtocloselyrelateddonorsandrecipients.
Furthermore, a substantial number of the transduction events (9%) are best described as gene
duplications that are mediated by mobile DNA vectors. We propose to distinguish this type of
paralogybythetermautology.Acomparisonofdonorandrecipientgenomesrevealedthatgenome
similarity is a superior predictor of species connectivity in the network in comparison to common
habitat. This indicates that genetic similarity, rather than ecological opportunity, is a driver of
successful transduction during microbial evolution. A striking difference in the connectivity pattern
ofdonorsandrecipientsshowsthatwhilelysogenicinteractionsarehighlyspecies-specific,thehost
range for lytic phage infections can be much wider, serving to connect dense clusters of closely
related species. Our results thus demonstrate that DNA transfer via transduction occurs within the
contextofphage–hostspecificity,butthatthistightconstraintcanbebreached,onrareoccasions,to
producelong-range LGTsofprofoundevolutionary consequences.
TheISME Journal (2017) 11,543–554; doi:10.1038/ismej.2016.116; published online20 September2016
Introduction is induced.Inthelyticcyclenewphagesare produced
using the host metabolism and are released during the
DNA transfer is an important mechanism for natural
host cell lysis (Campbell, 2003). The excision of phage
variationintheprokaryoticdomains(Ochmanetal.,
DNA from the host genome and the production of
2000). Recombination at the species level has a role
phages may be accompanied by packing of host DNA
in selective sweeps through the population (Shapiro
into the phages, which can then transfer it to the next
etal.,2012)whilelateralgenetransferacrossspecies
host in a process that has been termed transduction
boundaries has important implications to microbial
(ZinderandLederberg,1952).Specializedtransduction
adaptation and evolutionary transitions (for exam-
occurswhenthephageintegrasescleave,inadditionto
ple, Nelson-Sathi et al., 2015). Viruses that infect
the prophage, bacterial genes that are encoded at the
bacteria—termed phages—are known vectors of
prophage flanking regions. These are packed with the
DNA transfer between microbial cells (Zinder
phage DNA into the phages. Generalized transduction
andLederberg,1952;Halaryetal.,2010).Temperate
occurs when random bacterial DNA is packed into the
(or lysogenic) phages multiply via the lysogenic
phages (Weinbauer, 2004). A recent analysis of enter-
cycle, which is established by an integration of the
obacterialgenomesrevealedanextensivedomestication
phagegenomeintothehostchromosomes,creatinga
of genes encoded in prophages. The prophage domes-
prophage within the host genome. The phage
tication process comprises of rapid prophage inactiva-
typically remains dormant within the host and is
tion followed by a gradual genetic degradation that is
replicated with the host genome until the lytic cycle
marked by a strong purifying selectionontheacquired
gene sequence followed by their vertical inheritance
within the lineage (Bobay et al., 2014).
Correspondence: T Dagan, Institute of General Microbiology,
The frequency of transduction in nature may vary
Christian-Albrechts University of Kiel, Am Botanischen Garten
11,24118Kiel,Germany. between environments. In marine environment the
E-mail:[email protected] frequency of transduction has been estimated to
1Presentaddress:InstituteofQuantitativeandTheoreticalBiology, range between 1.33×10−7 and 5.33×10−9 transduc-
Heinrich-Heine-University of Düsseldorf, Universitätsstraße 1,
tants/plaque forming unit (Jiang and Paul, 1998).
D-40225Düsseldorf,Germany
In a freshwater environment, a higher transduction
Received10December2015;revised24June2016;accepted8July
2016;publishedonline20September2016 frequency ranging between 0.3×10−3 and 8×10−3
Phylogeneticrangeoftransduction
OPopaetal
544
transductants/plaque forming unit was observed, (version 10/2012). Annotation of 14920 prophages
where 20% of the gene recipients retained their encodedin8540genomicsequenceswasdownloaded
viability (Kenzaka et al., 2010). Phage lethality, as from PHAST database (version 10/2012) (Zhou et al.,
measuredbytheratioofphageinfectiontoadsorption, 2011). PHAST entries not found in GenBank were
and host specificity may however differ between discarded. This resulted in 9468 annotated prophages
various phage taxa. For example, cyanophages of encoded in 1330 complete and 1281 draft genomes.
diverse taxa are highly host-specific and their interac- Coding sequences (CDSs) in PHAST database are
tion is characterized by 100% lethality while hetero- classified into viral or bacterial according to their
siphoviruses have been shown to be adsorbed by a homology to known viral or bacterial genes (Zhou
wide range of Pseudoalteromonas strains and their et al., 2011). Prophages encoding only phage genes
lethality ranges between 10 and 40% (Deng et al., were excluded. The remaining 9201 (97%) prophages
2012). The realized host range in the wild is however encode 281616 CDSs, of which 89234 (32%) are
determinednotonlybythehostpermissibilitybutalso classified as bacterial genes. Prophages in the prokar-
by the phage–bacteria co-occurrence in the same yoticgenomeswereannotatedinadditionwithPhiSpy
geographic habitat (Flores et al., 2013). (Akhteretal.,2012)andVirSorter(Rouxetal.,2015a).
Bacteriaandtheirparasiticphagesareco-evolvingin Prophages were clustered into orthologous proph-
a constant arms race, yet their interaction may include age clusters (Bobay et al., 2013) based on their gene
also mutualistic aspects. The beneficial contribution of content. The first step included an all-against-all
phage-mediated gene transfer to the host fitness has BLAST of prophage protein sequences. Reciprocal
been documented in diverse environments (Canchaya best BLAST hits (Tatusov et al., 1997; Wolf and
etal.,2003).Forexample,genomesofphagesthatinfect Koonin, 2012) with E-value o1×10−10 were aligned
marine cyanobacteria have been found to encode globallyusingneedle(Riceetal.,2000).Pairshaving
components of both photosystem I (Sharon et al., o95% amino-acid identity were excluded. The
2009) and photosystem II (Lindell et al., 2005). The remaining CDSs were clustered into orthologous
elevated dose effect of these gene products within the protein families using MCL (Enright et al., 2002)
hostisassumedtoincreasethephotosystemsrecycling with default parameters. Pairwise prophage gene
efficiencyandbythatcompensatefortheenergeticcost contentsimilaritywascalculatedfromfrequenciesof
of phage proliferation (Lindell et al., 2005). Recently shared protein families using the Jaccard index.
sequenced metagenomic samples from hydrothermal Prophagesimilarityasestimatedusingthehypergeo-
ventsrevealedahighabundanceofphagesthatencode metricsimilaritymeasure(FuxmanBassetal.,2013)
components of the dissimilatory sulfite reductase gene was significantly correlated with the Jaccard mea-
(rdsr)(Anantharamanetal.,2014).Thisgeneisessential sure (r=0.96, Po10−15; excluding zero values) and
for sulfur oxidation and may confer an energetic we chose to proceed with the former. In agreement
advantage to chemolythoautotrophic bacteria that typi- with previous reports (Bobay et al., 2014), we found
cally inhabit such environments. In addition to the that high prophage gene content similarity thresh-
transfer of metabolic functions between closely related olds are too stringent, whereas low values lead to
hosts, phages have been found to mediate intergeneric clusters of distantly related prophages. Here we
gene transfer across species boundaries as exemplified apply a threshold of 470% shared families for the
in the transfer of toxin genes between Staphylococcus classificationoforthologousprophages(Bobayetal.,
aureus and Listeria monocytogenes in raw milk (Chen 2014). Our pipeline clustered 6494 (43.5%) pro-
andNovick, 2009). phages into 2397 orthologous prophages those are
Bacterialgenomesthatincludeaprophagemaybe considered as a phage entity in our evolutionary
considered as recipients in gene transfer events. reconstruction. The remaining 8426 unclustered
Bacterial genes in prophages are the result of gene prophages are designated as singleton phages.
acquisition by transduction and their origin can be
identified by homology and phylogenetic analysis.
Here we study the extent of phage-mediated gene Donor inference
transfer during microbial evolution using network The donor inference procedure operates within the
approach. The networks are composed of donors, framework of orthologous protein families and is
phages and recipients that are connected by recent assisted by a phylogenetic tree. Bacterial genes
transduction events reconstructed from genomic encoded within prophages were clustered into
data. Structural properties of the network supply a orthologous protein families. In the first step, we
large-scaleviewofbarriersfortransductionandgene searched for homologs to the prophage genes in the
transferability by phages in nature. GenBank genomes using BLAST. Using an E-value
o1×10−10threshold, weidentified 3908830 homo-
logous sequences to 75172 of the query genes,
Materials and methods whereas no homologs were detected for the remain-
ing 14197 (15.89%) queries. Protein pairs were
Data aligned globally with needle (Rice et al., 2000).
Genomes of 2103 complete and 1879 draft prokar- Protein sequences were considered as homologs if
yotic strains were downloaded from GenBank theyhadatleast90%globalamino-acidsimilarityto
TheISMEJournal
Phylogeneticrangeoftransduction
OPopaetal
545
the query sequence. This results in a data set of Genome similarity measures
42760 (57%) prophage bacterial genes and 252159 Genome sequence similarity (S ) between a recipi-
gs
homologs. The protein sequences were clustered ent and a donor was calculated as the Jaccard
intoorthologousproteinfamiliesusingMCL(Enright coefficient based on the proportion of 20bp
et al., 2002) to yield 20904 protein families with at segments common to the two genomes (using
leasttwoproteins.Proteinclusterscontainingatleast MUMmer (Kurtz et al., 2004)). Proteome similarity
three protein sequences (12611) were aligned using (S ) between bacterial species was similarly calcu-
pr
MAFFT (Katoh and Standley, 2013) and maximum lated as the Jaccard index of identical segments,
likelihood trees were reconstructed using PhyML3 but restricted to segments that have an overlap at
(Guindonetal.,2010)withtheLGmodel.Threetrees least of 10% within a protein-coding region. GC
were rooted using the midpoint criteria. In protein content similarity (S ) was calculated as follows:
GC
familieswithmultiplerecipientgenes,weexamined 100−|%GC −%GC | (Popa et al., 2011).
recipient donor
the monophyly of recipient genes, and when these The genome codon usage distance (D ) was calcu-
CU
wereparaphyletic(2205trees,17.48%),wetestedfor qlatffiPeffiffiffidffiffiffiffiffiaffiffiffisffiffiffiffiffitffihffiffiffieffiffiffiffiffiEffiffiffiuffiffifficffiffiffilffiiffiffidffiffiffieffian distance DCU ðdon; recÞ¼
the likelihood of an alternative tree with recipients n ðdon (cid:2)recÞ2 between the vectors of relative
consolidated into one clade. For this we used i¼1 i i
CONSEL (Shimodaira and Hasegawa, 2001) with codon frequencies per amino acid within the donor
the approximately unbiased test and the multi-scale and recipient genomes.
bootstrap technique. Of the reconfigured trees, 829
(37.6%) were not significantly less likely than the
original tree (approximately unbiased test, P⩾0.05) Synonymous and non-synonymous substitution rates
andwereretainedfordownstreamanalyses.Foreach The number of non-synonymous (dN) and synon-
gene acquired by transduction, we identified the ymous (dS) substitutions, and their ratio ω, was
most likely gene donor as the genome bearing a calculated using the branch model implemented in
homologous gene, which is the sister taxon of the PAML (Yang, 2007). We used ‘model 2’ of PAML,
acquiredgenein thephylogenetic tree or theunique allowing ω (dN/dS) to vary among the donor
homolog in the case of clusters with two members. branches, recipient branches and the remaining
Treeswherethedonorsorrecipientsformedaclade, (background) branches. For the special case of gene
rather than consisting of a single taxon, were duplication (self-donor recipient loops), we esti-
excluded from further analysis. An overview of the mated dN, dS and ω using the software package
various steps in the steps in the analysis is given in PAL2NAL(Suyamaetal.,2006).PAL2NALcreatesa
Supplementary Table S1. codon alignment from a pair of protein and their
corresponding DNA sequences and calculated the
d and d values using PAML (Yang, 2007). Codon
N S
Network construction adaptation index (Sharp and Li, 1987) for each gene
Donor–recipient relations were coded into the was calculated using the EMBOSS package (Rice
directed lateral gene transfer (dLGT) unipartite net- et al., 2000).
work,inwhichnodesrepresentbacterialspeciesand
edges lateral gene transfers mediated by transduc-
tion. Bacteria–phage relations were coded into a Function, habitat and ICTV classification
bipartite directed network, where nodes represent Functionalclassification ofeachcluster wasderived
eitherphages(5064nodes)orbacterialspecies(3982 from the Clusters of Orthologous Groups database
nodes).Thisenablespartitioningofthenetworkinto (Tatusov et al., 1997) by a majority vote of cluster
two subsets: the recipient subset (R-dLGT) that members. The habitat classification of donor and
consists of directed edges from prophages to their recipient nodes was extracted from the GOLD
hostbacterialspecies(thatis,therecipients);andthe database version March 2014. We defined 11 main
donor subset (D-dLGT) that consists of directed habitat classes using the combination of isolation
edges from donor bacterial nodes to phages (that is, place and ecology annotation. Putative International
the transfer vector). As in the dLGT, edge weights Committee on Taxonomy of Viruses (ICTV) annota-
correspond to the number of transferred genes. tion was derived from the majority vote of ICTV
labelsofprophagegenes.Laterallytransferredtoxin–
antitoxin operons were surveyed using PanDaTox
Network randomization (Amitai and Sorek, 2012) as a query. An additional
RandomizationofthedLGTnetworkwascarriedout survey for laterally transferred genes for antibiotics
using the switching methodology (Artzy-Randrup resistance was performed using the genes in CARD
and Stone, 2005), which rewires the weighted edges database (McArthur et al., 2013) as queries.
whilepreservingthein-andout-degreeofeachnode.
The method was implemented in an in-house
MatLab(versionR2015b;TheMathworksInc.)script Statistics and visualization
and used to generate 1000 randomly connected All statistical calculations were done using the
networks. Statistics toolbox in MatLab. The network layout
TheISMEJournal
Phylogeneticrangeoftransduction
OPopaetal
546
wascalculatedwithCytoscape(Shannonetal.,2003) orthologous prophages that are encoded within the
using the force-directed graph-drawing module. genomes of two chloroflexi strains: Chloroflexus
aurantiacus J-10-fl and Chloroflexus sp. Y400-fl.
The two prophages have a 100% match of their
Results
proteincontentwhenapplyingasequencesimilarity
threshold of 95% identical amino acids, thus they
The transduction network
areconsidered as orthologous prophages.One of the
To study the general properties of LGT by transduc-
tion we combine individual donor–recipient infer- eleven bacterial genes encoded in this prophage is
annotated as a threonine synthase and was puta-
ences into a network representation. Transduction
tively acquired by the phage from Chloroflexus
events are characterized by two distinct phases: the
aggregans(Figure1c).Thissmallclusterexemplifies
uptake of a gene from a donor into a phage and
how phage-mediated laterally transferred genes can
theacquisitionofageneaspartofaprophagebythe
enter the lineage.
recipient. We constructed a dLGT network that
includes two types of entities: bacteria and phages.
Adirectededgefromaphagenodetoabacterianode
Phage-mediated gene duplication—autology
designates a gene acquisition following transduction
Although most of the genes acquired via prophages
as inferred from the prophage annotation, where the
are xenologs, the network reveals a substantial
bacterianodeistherecipient.Adirectededgefroma
numberofgenes,wheretherecipientgenomeisalso
bacteria node to a phage node specifies the acquired
the donor (for example, Figure 1d). Thus, per
gene origin as inferred from the phylogenetic
definition these genes are paralogous rather than
analysis, where the bacteria node is the donor.
xenologous genes. We suggest terming such genes
We analyzed a total of 2103 finished and 1879 draft
autologs. According to our definition, an autologous
microbial genomes, including a total of 9468 anno-
geneistheresultofgeneduplicationthatismediated
tated complete and partial prophages annotated in
by a mobile DNA vector, where the donor is also
PHAST (Zhou et al., 2011). Alternative prophage
the recipient. Our analysis revealed 1550 (9%)
annotations calculated with PhiSpy (Akhter et al.,
autologous genes that are distributed over 543
2012)andVirSorter(Rouxetal.,2015a)resultedina
(21%) microbial genomes. About half of the gene
larger collection of putative prophages, but with
duplications in the network are of a single gene and
lesser consistency among the three methods than
up to a maximum of 48 genes in Magnetococcus sp.
PHAST (Supplementary Figure S1). Applying con-
MC-1. Of the self-donor recipients, 72% are con-
servative sequence similarity thresholds in the
nected to a single phage (Supplementary Figure S2).
different inference stages, our approach identified
A maximum of nine phages are connected by a self-
17158 transduction events, where donor and reci-
edge to Methylobacterium nodulans ORS_2060.
pient are specified. Constructing the network from
Of the 1550 autologs, 697 have no nucleotide
those events, where a single most likely donor was
substitutions at all, whereas the remaining 54%
identified yielded a dLGT network (Popa et al.,
autologs show the hallmarks of gene duplications.
2011), comprising 2573 bacteria and 4650 phage
They contain significantly more synonymous than
nodesthatareconnectedby15298edgessummariz- non-synonymous substitutions (Po10−15, using
ing all 17158 transduction events (Figure 1a; paired-Wilcoxon test). The median d /d ratio
N S
Supplementary Table S2). Edge weight in the net- (ω) is 0.12, which is significantly larger than the
work is calculated as the total number of genes observed for the bona fide gene acquisitions
that were transferred between the bacteria and (ω=0.09, P=3×10−13, using Wilcoxon test). More-
phage nodes. over, the codon adaptation index is significantly
The dLGT network comprises a large component smallerfortheprophagegenethanthegenomiccopy
of 4982 nodes, including 1538 bacteria and 3444 (P=0.035, using paired-Wilcoxon test). These obser-
phages.Theremainingnodesinthenetworkfallinto vationsareconsistentwiththeobservedrelaxationof
326 smaller clusters including, on average, three purifying selection in Entrobacteriales prophages
bacteriaandfourphages.Forexample,theNatrialba (Bobay et al., 2014).
magadiiϕCH1virushasatemperateinteractionwith
the chemoorganotrophic euryarchaeon N. magadii
isolated from Magadi lake in Kenya (Klein et al., Donor and recipient components
2002).Thevirusencodesatotalof24bacterialgenes. Structural properties of the dLGT network are the
Our inference algorithm yielded putative donors for result of two different phage–bacteria interaction
two of those genes. One gene, annotated as a modes. Phages connected to recipients represent a
hypothetical protein, was putatively acquired from lysogenic interactionthatinvolvestemperatephages
Halobiforma lacisalsi, an extreme halophylic and their hosts. Links between donors and phages
archaeon. The second gene, annotated as a gas are evidence for a lytic phage infection, where
vesicle protein, was putatively acquired from Natro- donors connected to the same phage designate the
nobacterium gregoryi, a haloalkaliphilic euryarch- putative hosts of that phage. The network thus
aeon (Figure 1b). In another cluster we identified combines two components: edges that link donors
TheISMEJournal
Phylogeneticrangeoftransduction
OPopaetal
547
bacteria
phage
Number of genes
10 20 30 40 50 60
Crenarchaeota Thaumarchaeota Bacteroidetes Chrysiogenetes Elusimicrobia Negativicutes Lentisphaerae Alphaproteobacteria
Euryarchaeota Deinococci Chlamydiae Cyanobacteria Fibrobacteres Fusobacteria Nitrospirae Betaproteobacteria
Korarchaeota Acidobacteria Chloroflexi Deferribacteres Bacilli Synergistetes Planctomycetes Deltaproteobacteria
Nanoarchaeota Actinobacteria Chlorobi Dictyoglomi Clostridia Ignavibacteria Poribacteria Epsilonproteobacteria
Tenericutes Thermodesulfobacteria Thermotogae Spirochaetes Verrucomicrobia Gemmatimonadetes Zetaproteobacteria Gammaproteobacteria
Natronobacterium gregoryi SP2
Halobiforma lacisalsi AJ5 Natrialba magadii ATCC 43099 / (289579696_1)
1
2
gas vesicle protein GvpA 1 2 hypothetical protein
Natrialba magadii ATCC 43099
Chloroflexus aggregans DSM 9485 Chloroflexus aurantiacus J 10 fl / (163845603_1)
1 bacterial gene
1 threonine synthase with inferred donor
Chloroflexus Y 400 fl bacterial gene
without inferred donor
Chloroflexus aurantiacus J 10 fl Chloroflexus Y 400 fl / (222523302_1) phage gene
Propionibacterium freudenreichii shermanii CIRM BIA1 / (297625198-3)
Propionibacterium 1
freudenreichii
shermanii CIRM BIA1 pyridine nucleotide-disulfide oxidoreductase 1
Figure 1 The directed transduction network. (a) Graphic representation of the directed, bi-partite, transduction network (dLGT)
reconstructedfromtransductionevents,whereasinglemostlikelydonorcouldbeidentified.Thenodescorrespondtobacterialgenomes
(triangles)andphages(hexagons).Bacterialnodesarecoloredbytheirtaxonomicgroup.Directededgescorrespondtogenesthatwere
transferredbetweenbacterialgenomesviaaphage.Bacteriatophageedgesdescribeuptakeofgenesduringlyticinfection(D-dLGT)while
phagetobacteriaedgesdesignatetheacquisitionofgenesinprophages(lysogenicinteractions,R-dLGT).Edgecolorcorrespondstothe
numberoftransferredgenes(seecolorbar).(b–d)DetailedexamplesshowingenlargedpartsthedLGT,includingbacterialspeciesnames
and the prophages genomic maps. Circled numbers identify specific genes in both views. Putative ICTV assignments are listed in
SupplementaryTableS3.
to phages correspond to gene uptake by phages, two components reveal the differential contribution
whereas edges that connect phages to recipients oflyticandtemperatephage–bacteriainteractionsto
correspond to gene acquisition by bacterial recipi- transduction dynamics during microbial evolution.
ents. Large-scale structural differences between the Thedonorandrecipientpartsofthenetwork,termed
TheISMEJournal
Phylogeneticrangeoftransduction
OPopaetal
548
here D-dLGT and R-dLGT, respectively, comprise a D-dLGT are significantly lower in comparison with
similar number of bacteria and phage nodes theR-dLGT(Po10−15,usingKolmogorov–Smirnovtest)
(Supplementary Table S2). Yet, the node connectiv- with medians of single gene per donor edge and two
ity degree is significantly larger in the D-dLGT in genes per recipient edge (Supplementary Table S2).
comparison with the R-dLGT for both phage and The R-dLGT comprises 230 (4.6%) edges with an edge
bacterianodes(D-dLGT,Po10−15;R-dLGT.,Po10−15, weight ⩾10. In the D-dLGT network, for comparison,
using Kolmogorov–Smirnov test; Figure 2a). we observe only 79 (0.77%) edges having an edge
Consequently the D-dLGT nodes are more densely weight ⩾10 (Supplementary Figure S4).
interconnected in comparison with nodes in the The different structural properties of the donor and
R-dLGT network. Most (93%) of the phages in the recipientnetworkcomponentssuggestthatgenetransfer
R-dLGTareconnectedtoasinglerecipientnodeand intohostsduringthelysogeniccycleusuallycomprises
at most to eight recipients (Figure 2a). Only 46% several genes, while gene uptake from hosts into the
(2131) of phages in the D-dLGT network are phagegenomeduringlyticinfectiontypicallyincludea
connected to a single donor node, while 25% single gene. Yet, the high connectivity of the D-LGT
(1146) phages are connected to two donors and the network indicates that the lytic interactions serve as
remaining29%(1373)phagesareconnectedtothree linkers among clusters of highly connected recipients,
donors or more (Figure 2a). thus they constitute an important contribution to the
It is noteworthy that highly connected phages global consolidation of the transduction network.
includegenesthathavethepotentialtobebeneficial
for the recipient. The most connected phage in the
D-dLGTnetworkencodes29genesofbacterialorigin
1
for which we identified 20 Enterobacteriales donors
(SupplementaryFigureS3;SupplementaryTableS4; 0.9
PhageID: 10223). The phage encodes the MazE/F
0.8 bacteria (R-dLGT)
toxin–antitoxinsystemthatcanmediatescellgrowth bacteria (D-dLGT)
arrestandwasshowntoincreasethepersistenceand cy 0.7 phage (R-dLGT)
n
survival of Escherichia coli under antibiotic stress e phage (D-dLGT)
qu 0.6
(Zhang et al., 2005). Our analysis further uncovered e
the transfer of 73 toxin–antitoxin genes mediated by e Fr 0.5
v
3d2ucptihoangeesve(nStuspspulegmgeesntttahraytTpahbalgeesS5m).ayTheensceodtreanfosr- ulati 0.4
m
addiction mechanisms similarly to plasmids. Cu 0.3
Another highly connected phage in the D-dLGT
networkisconnectedto19Enterobacterialesspecies 0.2
(SupplementaryFigureS3;SupplementaryTableS4;
0.1
PhageID: 11150). The phage encodes for mdtH, a
multidrug resistance gene that confers resistance to 0
100 101 102
norfloxacin and enoxacin (Nishino and Yamaguchi,
Node Connectivity Degree
2001). Transferred genes in the network include
additional 46 genes coding for a broad range of
antibiotic resistance (Supplementary Table S6)
23
demonstratingaputativeroleofphagesinthespread kingdom 29
of antibiotics resistance. The most connected phage 2 D-dLGT
dLGT
in the R-dLGT is connected to eight Bacillus 39 R-dLGT
recipients (for details see Supplementary Figure S3; phylum 50
1
Supplementary Table S4; PhageID:5009). The phage
containseightgenesofbacterialorigin.Oneofthose, 26
class 37
bclA, encodes for a spore surface glycoprotein in 1
Bacillus anthracis (Sylvestre et al., 2002).
320
The different bacteria and phage connectivity order 497
11
pattern of the R-dLGT and D-dLGT is evident also
in their global structure. The D-dLGT contains 782
genus 1372
significantly less connected components in compar-
57
ison to the R-dLGT. Furthermore, nodes in the
1329
recipient network are clustered into significantly species 2664
smaller components in comparison with the donor 261
network (P=7×10−10, using Kolmogorov–Smirnov 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
test) and the number of nodes in the D-dLGT largest Phage frequency
component is 25-fold larger in comparison with that
Figure 2 Connectivity and taxonomic distribution. (a) Cumula-
of the R-dLGT largest component (Supplementary
tive distribution of node connectivity degree. (b) Taxonomic
Table S2). In consequence, edge weights in the distributionofdonorsandrecipients.
TheISMEJournal
Phylogeneticrangeoftransduction
OPopaetal
549
Host range in the transduction network infections. Many of these examples are observed in
Phages that are linked to more than one donor or microbial genomes sequenced as part of the human
recipient in the network supply an insight into the microbiome project, thus it is possible that the high
phage host range. In the D-dLGT component, about sampling density of that habitat facilitated that
half of the phages (2519, 54.17%) are connected to recovery of those rare interactions.
multipledonors.Mostofthesephagesareconnected
to donors of the same species (1329, 53%) or genera Barriers for gene transfer by transduction
(782, 31%) revealing a very narrow taxonomic host The majority of phages (2664, 64%) connect donors
range at the donor side (Figure 2b). Only 22 phages and recipient from different strains of the same
in the D-dLGT network are connected to two donors species (Figure 2b). These phage-mediated DNA
that are members of different phyla, 20 of which are transfers are best viewed as genetic recombination
connected to Firmicutes strains (Supplementary rather than lateral gene transfer events. A Siphovir-
Table S7). A single phage is connected to three idae phage connected to multiple Vibrio cholera
donors from different phyla including Bacteroides strains illustrated this phenomenon (Supplementary
sp.3_1_33FAAandClostridiumsp.M62/1thatwere Figure S3; Supplementary Table S4; PhageID: 8390).
isolated from the human gastrointestinal tract and The phage encodes the nqr operon that has an
Cardiobacterium hominis ATCC 15826 (Gammapro- important function in the bioenergetics and home-
teobacteria) that was isolated from the human ostasis of V. cholerae (Barquera et al., 2002).
cardiovascular system (Human Microbiome Project ThefrequencyofobservedLGTsdecreasesmarkedly
Consortium, 2012) (Supplementary Table S7; when the donor–recipient taxonomic separation
PhageID: 9283). increases(Figure2b).Attheinter-domainlevel,only
Phages connected to more than a single recipient a single phage was observed, connecting Methano-
intheR-dLGTnetwork(333;7%)showevenstronger brevibacter smithii DSM 2374 as the recipient with
species-specificity, with most phages (261, 78%) Bacillus cereus Rock3–28 as the donor
connected to recipients classified into the same (Supplementary Figure S3; Supplementary Table
species (Figure 2b). A total of 57 (17%) phages S4;PhageID:9888).Therecipientstrainwasisolated
are linked to recipients from different species from human feces (Human Microbiome Project
within the same genus (Supplementary Figure S3; Consortium, 2012), whereas the donor was isolated
Supplementary Table S4; PhageIDs: 5548 and 5273). from the soil (Zwick et al., 2012). The prophage
Only 11 (3.3%) phages are found in recipients of includes a gene encoding for tetracycline resistance
different genera within the same taxonomic order. that has 100% identical amino acids to the gene
The rare inter-generic transduction events include a encoded in B. cereus. To our knowledge, this is the
phage connected to two Clostridiales recipients: first genomic evidence for transduction of an
Blautia hansenii DSM 20583 and Ruminococcus archaebacterium by a eubacterial bacteriophage;
gnavus ATCC 29149 (Supplementary Figure S3; hence,this putative inter-domain transfer represents
Supplementary Table S4; PhageID: 5915). Both a very exceptional event.
strains were isolated from the human digestive Barriers for transduction may be related to the
system (Human Microbiome Project Consortium, genetic requirements for a successful gene acquisi-
2012), hence they probably share a common habitat. tion and the ecological co-occurrence of the
A single phage links two recipients from different connected partners. In contrast to transduction, in
classes within the Firmicutes phylum, Clostridium transformation and conjugation the integration of
M62-1 and the Lactobacillus ruminis ATCC 25644, acquiredDNAintotherecipientgenomeismediated
both isolated from the human gastrointestinal tract by homologous recombination and therefore
(Human Microbiome Project Consortium, 2012; depends on sequence similarity between the donor
Supplementary Figure S3; Supplementary Table and recipient (Thomas and Nielsen, 2005). During
S4; PhageID: 6299). Only two phages link to transduction, however, the acquired DNA is inte-
recipients from different phyla (Supplementary grated into the recipient genome using the phage
Figure S3; Supplementary Table S4; PhageIDs: mechanism (Thomas and Nielsen, 2005), hence no
5805 and 6260). One of those connects the Bifido- such dependency is expected. To test for genetic
bacteriumpseudocatenulatumDSM20438(phylum: barriers to DNA transfer by transduction we
Actinobacteria) and Parvimonas micra ATCC 33270 calculated the genome similarity between donors
(phylum: Firmicutes), both isolated from the human and recipients using four measures. Genome simi-
gastrointestinal tract (Human Microbiome Project larity (S ) is calculated as the Jaccard index of
GS
Consortium, 2012). identical ⩽20bp sequences between the donor and
The narrow taxonomic range of multiple donors recipient genomes. CDS similarity (S ) is calcu-
CDS
and recipients observed in the network components lated similarly but is restricted to protein CDSs.
is in agreement with experimental observations of Codon usage distance (D ) is calculated as the
CU
phagespecies-specificity(HymanandAbedon,2010) Euclidean distance between the relative codon
and is expected from the tight phage–host frequencieswithinthedonorandrecipientgenomes.
co-evolutionary dynamics. Our results reveal how- GC content similarity (S ) is calculated from the
GC
ever several genomic footprints of rare cross-species genomic content of guanine and cytosine in the
TheISMEJournal
Phylogeneticrangeoftransduction
OPopaetal
550
donor and recipientgenomes. The distribution of all 0.87; FPR 0.08). Restricting the analysis to a subset
similarity measures was compared between the of prophages that were detected also by PhiSpy
dLGT network and a set of 1000 networks, where (Akhteretal.,2012) orVirSorter (Roux etal.,2015a)
the edges have been randomly shuffled. prophage annotation tools reveals an even
We find that donors and recipients connected in sharper deviation from the expected by chance
the dLGT network are significantly more similar to (Supplementary Figure S5). Our results demonstrate
each other than expected by chance using all that low donor–recipient genome similarity is an
similarity measures (Figure 3). The four similarity important barrier that constrains the extent of LGT
measures are correlated—closely related genomes via transduction.
will score high on each measure, yet it is of interest Another possible barrier for transduction is the
to grade their importance as barriers to LGT. To this need for ecological co-occurrence of donors, phages
end, we consider each pairwise similarity measure and recipients. This barrier may be partially brea-
as a predictor of the connectedness state of the pair ched by phage mobility that is thought to enable the
of species, and conduct a receiver-operating char- transfer of genetic material between donors and
acteristics analysis (for example, Fawcett, 2006). We recipients across a larger spatial separation com-
find that genome similarity is the best predictor for paredwithotherLGTmechanismthataredependent
dLGT connectedness, with an area under the of physical proximity (Majewski, 2001). Donor–
receiver-operating characteristics curve (AUC) of recipient pairs share the same habitat in 3330
0.99, and an optimal discrimination of 0.97 true- (44%) cases, of which the largest group (1383,
positiverate(TPR)and0.03false-positiverate(FPR). 41%) are members of the ‘human-associated’ habitat
Thenextbestmeasureiscodonusagedistance(AUC group. Inthe remaining 4187 (56%) donor–recipient
0.98; TPR 0.93; FPR 0.03), followed closely by CDS pairs classified in different habitat groups (cross-
similarity (AUC 0.96; TPR 0.94; FPR 0.04). habitat transfer events), we observed the majority
GC content similarity is an inferior predictor in (858, 20.49%) of links between the donor group
comparisonwiththeothermeasures(AUC0.95;TPR ‘host’ and the recipient group ‘human-associated’
1 1
0.9 0.9
0.8 0.8
y
c
en 0.7 0.7
u
eq 0.6 0.6
Fr
e 0.5 0.5
v
ati 0.4 0.4
ul
m 0.3 0.3
u
C
0.2 0.2
0.1 0.1
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Genome Similarity Coding Sequence Similarity
1 1
0.9 0.9
0.8 0.8
y
c
en 0.7 0.7
u
eq 0.6 0.6
Fr
e 0.5 0.5
v
ati 0.4 0.4
ul
m 0.3 0.3
u
C
0.2 0.2
0.1 0.1
0 0
0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6
%GC Similarity Codon Usage Distance
Figure 3 Donor–recipient genome similarity. Cumulative distribution of donor–recipient genome similarity measures in the dLGT
network(green)and1000randomizednetworks(black).(a)Genomesequencesimilarity,(b)CDSsimilarity,(c)GCcontentsimilarityand
(d)codonusagedistance.
TheISMEJournal
Phylogeneticrangeoftransduction
OPopaetal
551
bacteria (Supplementary Figure S6). To evaluate Thenucleotidesubstitutionpatternofgenesinthe
whether these values are different from the expecta- transductionnetworkindicatesthattheiracquisition
tion given that habitat sampling is heavily skewed was very recent or that they evolve under extremely
toward certain habitats, we estimated the expected strong purifying selection. Half (52%) of the donor–
within- and cross-habitat frequencies from 1000 rando- recipient pairs have no nucleotide substitutions at
mized dLGT networks. Links between donors and all. Comparing the nucleotide substitution rate
recipients from the same habitat are significantly over- between the donor and recipient lineages for the
represented in the dLGT network, with the corollary remaining 48%, we observe a very slight and not
that most cross-habitat links are occurring at a lower significant increase in the recipient lineage rates
frequency than expected. However, some habitats do (Supplementary Table S8). The ω ratio is also not
show a higher than expected cross-habitat LGT significantly different between the two lineages and
frequencies (Supplementary Figure S6). For example, in 95% of the genes is below 0.5 in both lineages.
we found 51 (expected 17) links between ‘soil and Togethertheseobservationssuggestthatthestrength
sediment’ and the ‘plant’ group and 41 (expected 19) of purifying selection in recipient lineages remains
transfers between the habitat ‘plant’ and ‘soil and similartothatinthedonorlineageswithnoapparent
sediment’ group. Forty-eight percent of these transfers relaxation of selective constraints or nonfunctionali-
are intra-specific and 95% are intra-generic. Indeed, zation. The great majority (95%) of bacterial genes
habitat sharing is only a weak predictor of species thatareencodedinprophagesaresingle-copygenes,
connectedness,withequivalentAUCofonly0.64(TPR that is, there is no pre-existing homologous gene in
0.44; FPR 0.17). Our analysis thus reveals that the the recipient genome. Taken together with the
barriers for gene transfer via transductionare primarily evidence for gene functionality, this suggests that
genetic while ecological barriers have a smaller role. mosttransductioneventsresultinanacquisitionofa
newfunction.Furthermore,itcouldindicatethatthe
accessibility of the host to the new function is
Functional classification and evolutionary constraints maintained as long as the lysogenic interaction with
The functional composition of dLGT genes is the phage is maintained.
significantly different than that of the analyzed
bacterial genomes (Po10−15, using χ2-test). Informa-
tion processing functions are overrepresented in the
Discussion
network,whereascellularprocessesandmetabolism
functions are depleted (Supplementary Figure S7). Here we study the contribution of phage-mediated
Of the genes that could be classified into putative gene transfer to microbial genome evolution. The
functions (2274, 13%), 42% perform metabolism transduction network reconstruction revealed a
functions, whereas 35% were involved in informa- substantial frequency of autologs. Autologs may be
tion processing; most of those are annotated as the result of recurrent infections, where both donor
transcription genes. Another 23% of the genes were andrecipientaremembersofonelineage.Transduc-
classified into cellular processes, with a majority of tion within the lineage may thus contribute to
cell wall and membrane biogenesis function protein family expansion in bacteria, which was
(Supplementary Figure S7). Interestingly, informa- proposed to be mediated more often by LGT than by
tion genes are transferred between less similar gene duplication (Hooper and Berg, 2003; Treangen
donors and recipients than the other functions, and Rocha, 2011). Indeed, we have to assume that
while metabolism and cellular processes genes are low sampling density may obscure a gene donor
transferred between equally similar donors and amongcloselyrelatedstrains,whosegenomehasnot
recipients (α=0.05, using Tukey test). This observa- yetbeensequenced.Nevertheless,thehighsequence
tion may be attributed to the universality of similarity and lack of alternative homologs besides
information processing genes.In addition,we exam- the recipient genomic copy indicate that autologs
ined whether the functional categories distribute originate from within the pan-genome.
differently when considering the habitat of the host Thetopologicaldifferencesbetweenthedonorand
or the taxonomic classification, with the difference recipient network components suggest that host-
that the analysis was restricted to the three main specificity is much more prevalent in lysogenic
functional classes due to limited sample size (using interactions and that phages have a broader host
χ2-testandfalsediscoveryrate(FDR)of5%).Noneof range for lytic infection. Because lysogenic phages
the habitats was found to be significantly different are highly dependent on the host cellular processes
from the others in terms of the functional classifica- (for example, Tal et al., 2014) it is likely that closely
tion of the dLGT genes. Only two taxonomic groups related strains having a similar genetic background
appeared to deviate from the common functional can have a lysogenic interaction with the
distribution. In Actinobacteria we observed an same phage.
excess of metabolic genes and paucity of cellular Previous studies of LGT dynamics estimated that
processes genes; in Epsilonproteobacteria we most LGT events involve very few genes while bulk
observed an excess of cellular processes genes and transfers are relatively rare (Kunin et al., 2005; Popa
paucity of information genes. et al., 2011). Transduction, however, is known to
TheISMEJournal
Phylogeneticrangeoftransduction
OPopaetal
552
D1 D1
D1
D2 D2 D2
R R Phage R
D3 D3 D3
D4
D4 D4
D1 D2 D3 D4 R
Figure4 Observation,networkandinterpretation.(a)Therawobservationsaregenesencodedinprophageswithinrecipientgenomes(R)
andtheirhomologsinmicrobialgenomesthatareidentifiedasthedonors(D1–D4).(b)Unipartitegenome-onlynetworkcodingofthe
donor–recipientrelationships.(c)Bipartitecoding,includingthephageasthemobileelementfacilitatingthetransfer.Notethebundlingof
multiple edges into a single transduction event. (d) One possible interpretation of the network in c as a historical scenario, depicting
consecutivegeneuptakeduringlytic(D1,D2andD4)orlysogenic(D4andR;penetratingarrow)cycles,followedbytransductionintothe
recipient.Theprophageintherecipientisacurrentsnapshotofanongoingprocess.
involve several genes simultaneously. This is appar- recipient; hence, the translational barrier for their
ent in our network approach and its distinction adaptation is expected to be rather low. Several
between donors, mobile elements and recipients. temperate phages have been reported to encode
This inclusion of the mobile element allows us to genes that are transcribed independently from the
groupgenesthatarrivedintotherecipientinasingle prophage excision mechanism (Cumby et al., 2012).
transduction event, as opposed to genome-only Thus, the transcriptional regulation of genes
networks that would depict the same history as acquired via transduction is likely to be promoted
multiplicity of edges into the recipient (Figure 4). by prophage-encoded promoters so that genes
On the donor side, on the other hand, we still acquired by transduction are functional upon
observe multiple donors for the same transduction acquisition.
event. Hence, transduction is characterized by bulk In summary, our results demonstrate that LGT via
acquisition of mosaics of genes from multiple transduction occurs within the tight constraints of
donors. phage–host specificity. Consequently, transduction
The dLGT network reveals the existence of strong is probably more important in the evolutionary
taxonomic and genetic barriers for phage-mediated context of genetic recombination within the species,
lateral gene transfer. Previous studies advanced the and selfing in the case of autologs, than in the
viewthatgenetransferduringmicrobialevolutionis evolutionary context of long-range gene transfers
largely determined by ecological rather than phylo- between distinct lineages. LGT is commonly viewed
geneticfactors(Smillieetal.,2011).Althoughwedo asasourceforreticulatedeventsthatreducethetree
find an overrepresentation of transfers within habi- signal during prokaryotic evolution (Martin, 1999).
tats, habitat sharing is only a weak predictor of Our current results show that the reticulated events
species connectedness and is inferior to all introduced by transduction affect mostly clades of
sequence-derived similarity measures. The signifi- closely related species and very rarely do they
cant high codon usage similarity of donors and traverse the tree and disrupt its global topology.
recipientsisconsistentwithpreviousobservationsof
non-randomcodonusageinphagegenomes,leading
to the suggestion that phage codon usage is adapted Conflict of Interest
to that of the host (Sharp et al., 1985; Roux et al.,
2015b). Previous studies have emphasized the The authors declare no conflict of interest.
importance of codon usage similarity for LGT,
suggesting that high codon usage similarity between
Acknowledgements
acquired genes and the recipient genome will
increase the xenolog retention prospects (Medrano-
We thank Anne Kupczok, Nils Hülter and Elie Jami for
Soto et al., 2004; Tuller et al., 2011). Our results their helpful comments on the manuscript and acknowl-
reveal that genes acquired via transduction originate edge support from the European Research Council (Grant
in genomes with a similar codon usage to that of the No.281357).
TheISMEJournal