Regulation of genetic flux between bacteria by restriction–modification systems

Data.

Curated reference protein sequences of restriction endonucleases (REases) and methyltransferases (MTases) belonging to Types I, II, IIC, and III R-M systems and Type IV REases were downloaded from the dataset “gold standards” of REBASE (26) (last accessed in October 2014). For purposes of building a recognition target predictor of R-M systems, we excluded R-M proteins having more than one recognition motif or recognition motifs smaller than 4 bp, control proteins involved in the regulation of expression of certain REases, nicking enzymes, homing endonucleases, truncated proteins, and orphan MTases for which no R-M type could be assigned either by searching REBASE or the literature. We searched for R-M systems in 884 bacterial genomes from 79 bacterial species (Dataset S1). The sequences and corresponding annotations of these genomes were retrieved from GenBank Refseq (ftp://ftp.ncbi.nih.gov/genomes; last accessed in February 2014) (41). We excluded genes indicated in the GenBank files as partial genes, as well as those lacking a stop codon or having one stop codon within the reading frame. The clades analyzed were based on taxonomy, i.e., the genomes of a named species were put together. In some rare cases, like Escherichia coli and Bacillus cereus, this does not exactly define monophyletic clades (because of Shigella spp. and Bacillus anthracis, respectively).

Identification of Core Genomes and Pangenomes.

We built core genomes for species with at least four complete genomes available in GenBank RefSeq (Dataset S1). We used a methodology previously published (42). Briefly, a preliminary list of orthologs was identified as reciprocal best hits using end-gap free global alignment, between the proteome of a pivot (typically the first completely sequenced isolate) and each of the other strain’s proteomes. Hits with less than 80% similarity in amino acid sequence or more than 20% difference in protein length were discarded. This list of orthologs was then refined for every pairwise comparison using information on the conservation of gene neighborhood. Because (i) few genome rearrangements are observed at these short evolutionary distances, and (ii) HGT is frequent (28), genes outside conserved blocks are likely to be xenologs or paralogs. Thus, positional orthologs were defined as bidirectional best hits adjacent to at least four other pairs of bidirectional best hits within a neighborhood of 10 genes (5 upstream and 5 downstream). These parameters (four genes being less than one-half of the diameter of the neighborhood) allow retrieving orthologs on the edge of rearrangement breakpoints and therefore render the analysis robust to the presence of rearrangements. The core genome of each species was defined as the intersection of pairwise lists of positional orthologs.

Pangenomes are the full complement of genes in the species and were built by clustering homologous proteins into families for each of the 79 species. We determined the lists of putative homologs between pairs of genomes (including plasmids) with BLASTP v.2.2.28+ (default parameters) and used the e values (<10⁻⁴) to cluster them using SILIX (v1.2.8, lbbe.univ-lyon1.fr/SiLiX) (43). SILIX parameters were set such that a protein was homologous to another in a given family if the aligned part had at least 80% of identity and if it included more than 80% of the smallest protein.

Phylogenetic Analyses.

We built a tree to display the phylogenetic distribution of our dataset using the 16S rRNA sequences of the sequenced type strain of the 79 analyzed bacterial species (Fig. 1A and Dataset S2A). We made a multiple alignment of the 16S sequences with MAFFT, followed by manual correction with SEAVIEW. The tree was computed by maximum likelihood with PHYML v3.0 (44) under the general time reversible (GTR)+Γ(4)+I model. This tree is never used in the calculations; it is only used in Fig. 1A to display the relative position of each species in the phylogeny of bacteria.

We made core genome trees for each clade using a concatenate of the multiple alignments of the core genes (available upon request). Each species tree was computed with RAxML v8.00 (45) under the GTR model and a gamma correction (GAMMA) for variable evolutionary rates. All trees are shown in Dataset S2B. We performed 100 bootstrap experiments on the concatenated alignments of each clade to assess the robustness of the topology of the tree. The vast majority of nodes were supported with bootstrap values higher than 90% (Dataset S2B). We inferred the root of each phylogenetic tree using the midpoint-rooting approach of the R package “phangorn” v1.99.14 (46).

The phylogenetic depth was defined as the average root-to-tip distance, and was computed as the diagonal mean of the phylogenetic variance–covariance matrix of each tree, using the vcv.phylo function in the R package “ape.” We retrieved the patristic distances between taxa using the cophenetic.phylo function from the same R package.

Inference of HR.

We inferred HR on the multiple alignments of the core genes of each species using several programs. Multiple alignments were made on protein sequences using MUSCLE v3.8 with default parameters (47), and then backtranslated to DNA, to increase accuracy. They were not edited to avoid biasing the inference of HR rates. We computed for each multiple alignment the neighbor similarity score (NSS) (48), the maximum χ² (MaxChi) method (49), and the pairwise homoplasy index (PHI) test using PhiPack (50) (downloaded in February 2015) with 10,000 permutations. The Geneconv v1.81 program (51) was used with the options /w123 to initialize the program’s internal random number generator, /lp to compute the local pairwise P values, and the parameter gscale = 2 to allow for mismatches between fragments. Geneconv outputs “inner” fragments that have arisen by gene conversion events between ancestors of aligned sequences, and “outer” fragments that represent either conversion events involving an ancestor outside the alignment or events that have been eroded by subsequent mutations. Only inner fragments were considered in this study. Besides, only gene conversion tracts with P < 0.05 in the global analysis were used. To analyze the exchanges between specific taxa within clades, we used the fragments of Geneconv, because it is the only program producing this information. We also identified events of HR with ClonalFrameML (CFML) v10.7.5 (3) with a predefined tree, default priors R/θ = 10⁻¹, 1/δ = 10⁻³, and ν = 10⁻¹ and 100 pseudobootstrap replicates, as previously suggested (3). Mean branch lengths were computed with the R package “ape” v3.3 (52), and transition/transversion ratios were computed with the R package “PopGenome” v2.1.6 (53). The priors estimated by this mode were used as initialization values to rerun CFML under the “per-branch model” mode with a branch dispersion parameter of 0.1.

Reconstruction of the Evolution of Gene Families.

We assessed the dynamics of gene family repertoires using Count (25) (downloaded in April 2015). This program offers the most general phylogenetic birth-and-death models currently available in the literature. It models the gains of novel families, as well as expansion and contraction of existing families, while accommodating rate variations across phylogenetic lineages and across families. The analysis with Count starts with the estimation of the parameters of the model by maximum likelihood using the pangenome matrix of gene presence and absence. Count then uses these parameters to calculate the expected size of each family in every internal node of the species tree. It also computes the expected number of gain, loss, expansion, and contraction events along each branch.

Rates were computed with default parameters, assuming the Poisson family size distribution at the tree root, and uniform gain, loss, and duplication rates. One hundred rounds of rate optimization were computed with a convergence threshold of 10⁻³. After optimization of the branch-specific parameters of the model, we performed ancestral reconstructions by computing the branch-specific posterior probabilities of evolutionary events, and inferred the gains in the terminal branches of the tree.

The manual of Count provides no guidelines on which threshold probabilities to use in a branch. We consulted the literature and found very different values (ranging from 0.2 to 0.95) and methods to compute these (54, 55). To have an objective way of defining this threshold, we computed the gain/loss scenarii using the Wagner parsimony (same parameters, relative penalty of gain with respect to loss of 1). We then made a correlation analysis of the number of HGT events inferred by maximum likelihood and under Wagner's parsimony using posterior probability thresholds of 0.2, 0.3, 0.5, and 0.9. The correlations were always very high (Spearman’s ρ > 0.87), and maximal for 0.2 (Spearman’s ρ = 0.96, P < 10⁻⁴) (Fig. S6C). Hence, the posterior probability matrix was converted into a binary matrix of presence/absence of HGT events using a threshold of 0.2. At a given terminal branch, the expected total number of acquisitions was computed by summing all family-specific gene gains obtained from the posterior probability binary matrix. Transfer co-occurring in pairs of terminal branches in a tree, was computed by summing all common pairwise gene acquisitions.

To control for the effects of the choices made in the definition of our model, we made three complementary analyses. First, we verified that genetic flux was higher for pairs of genomes encoding cognate R-M systems than for the others using different thresholds for the branch-specific posterior probabilities of evolutionary events (0.5: P < 0.01, Mann–Whitney–Wilcoxon test; 0.9: P < 0.05, Mann–Whitney–Wilcoxon test). Second, we allowed for different gain–loss and duplication–loss rates within branches and found very similar distributions of gene transfer events (always Spearman's ρ > 0.90, P < 10⁻⁴). Third, we allowed rates to follow a gamma distribution with four categories and found very similar values of gene transfer events (Spearman's ρ > 0.88, P < 10⁻⁴).

Identification of R-M Systems.

Identification of R-M systems was performed as previously described (4). Briefly, all-against-all searches were performed for REase and MTase standard protein sequences retrieved from REBASE using BLASTP v2.2.28+ (default settings, e value < 10⁻³). The resulting e values were log transformed and used for clustering into protein families by Markov Clustering (MCL) v10-201 (56). Each protein family was aligned with MAFFT v7.205 (57) using the E-INS-i option, 1,000 cycles of iterative refinement, and offset 0. Alignments were visualized in SEAVIEW v4.5.3 (58) and manually trimmed to remove poorly aligned regions at the extremities. Hidden Markov model (HMM) profiles were then built from each multiple sequence alignment using the hmmbuild program from the HMMER v3.0 suite (59) (default parameters). Type II MTases were retrieved using the PFAM-A profiles PF01555.12, PF02086.9, PF00145.1, and PF07669.5 (last accessed in February 2013). Both Type II and Type IV REases diverge rapidly, resulting in sequences that produce poor multiple alignments (6) that cannot be used to build protein profiles. In these cases, BLASTP was used to scan the genomes for homologs (default settings, e value < 10⁻³, and minimum coverage alignment of 50%). Types I, II, and III R-M systems were identified by searching genes encoding the MTase and REase components at less than four genes apart. The output was subsequently curated to eliminate multiple occurrences of the same R-M system (this occurs when two REase or MTase genes are encoded in the same locus). R-M systems containing more than one Specificity (S) gene were considered as a single system. Situations involving ambiguous identification may also occur, for example, between REases of Types II and IV, or between Type IIC systems and other MTases or REases. In these cases, the R-M type was defined on the basis of the corresponding genomic context (presence or not of a linked REase or MTase) and on the output of the analysis of the system using REBASE. Type IIC R-M systems were defined as those including a gene encoding both an MTase and a REase function with similarity to Type IIC MTases and REases.

Robustness of the Recognition Target Motif Predictor.

We tested the robustness of our recognition target motif predictor in two ways.

Type I and III systems showed rapid change of sequence recognition in terms of sequence divergence. This precluded the inference of cognate systems in these types. We made an additional analysis to check whether these problems were due to the analysis of all of the sequences of the MTases and REases, which may evolve fast and recombine. We therefore focused our analysis on the evolution of the target recognition domains (TRDs), which are directly involved in the recognition of the site (Fig. S3A). Type I Specificity domain sequences were retrieved from the REBASE gold standard database (Table S2), whereas Type III TRDs were selected from our dataset of Type III mod genes. The former are typically flanked by the classical MTase specific motifs IV–VIII at the N terminus of the mod core and the I–III and X motifs at the C terminus (60).

Since the previous test confirmed our earlier results, we made a second analysis where we increased the size of the initial dataset by retrieving additional MTase and REase sequences from the REBASE PacBio database (rebase.neb.com/rebase/rebpbe.html) (Fig. S3 B–E). At the time it was last accessed (January 2016), REBASE PacBio had information on 642 organisms corresponding to a total of 2,446 PacBio records. More than 50% of the recognition sequences had not yet been assigned to known enzymes, more than 26% of the 642 sequencing projects were still not classified as complete (still at “shotgun” phase), and many MTases (and cognate REases) had not been thoroughly characterized from the biochemical point of view (i.e., were not classified by REBASE as gold standards). To minimize the risk of adding a large number of eventually poorly characterized protein sequences and/or recognition motifs, we retrieved only the sequences pertaining to the eight clades analyzed in Fig. 4 excluding data not yet assigned to known enzymes and genomes from unfinished sequencing projects. This resulted in 213 Type I, II, and III MTase and REase sequences from 54 genomes. Importantly, this allowed to increase the dataset of MTases in ∼38% and of REases in ∼70% (for Type I and Type III). Nevertheless, the results remained unchanged (Fig. S3).

Robustness of the Count Analysis to Phylogenetic Reconstruction.

The robustness of HGT analysis using Count depends on the quality of the phylogenetic reconstruction. We checked the robustness of the results of Count in the light of phylogenetic uncertainty in two ways.

• i) We used the 100 bootstrap trees mentioned above (SI Methods, Phylogenetic Analyses) of each of the eight clades analyzed in Fig. 4. We then analyzed the matrix of presence and absence of genes in the pangenome with Count for every single bootstrap tree. This analysis was done exactly as the one of Fig. 4. For each set of bootstrap experiments (eight bootstrap trees, one per clade), we computed the difference between the median HGT gains among R-M cognate genomes and R-M noncognate genomes (Δ_medianHGT). The analysis of this distribution showed that, in the 100 experiments, there was no single value smaller than zero (Fig. S7). This means that all sets of bootstrap experiments showed an overrepresentation of HGT between cognate genomes (relative to noncognate). Hence, incertitude in the phylogenetic reconstruction does not affect our conclusions.

• ii) The phylogenetic trees that we have used in our work have not been modified to account for distortions associated with recombination. Hence, we also tested the robustness of the method by running Count using the trees obtained from ClonalFrame (CF) v1.2 (61). This program provides rooted trees purged from recombination; hence it allows to test the effect of the use of sequences not purged from recombination and of midpoint rooting. For each core genome of our eight-clade dataset, we ran CF (200,000 iterations with the first half discarded as burn-in, and the second half sampled at each iteration). This tree was then used in the Count analysis. For each clade we ran at least three more independent runs of CF to check convergence and mixing between runs. This was then manually checked. We then reran Count using the CF rooted trees (using the same method as for Fig. 4). The results were qualitatively simple, because R-M cognate genomes had twice the frequency of HGT as the others (an average of 10 and 5 events, respectively; P < 0.001).