Escherichia coli strain RR1 (= KCTC 2134 or ATCC 31343) for genome sequencing was purchased from the Korean Collection for Type Cultures (Jeongeup, Jeollabuk-do, South Korea). Cells were grown aerobically in LB medium at 37°C. Genomic DNA was isolated using the Wizard genomic DNA purification kit (Promega, Madison, WI, USA). Library construction using Illumina TruSeq DNA sample preparation kit v2 and 101 cycle paired-end sequencing using the Illumina HiSeq 2000 system were performed according to the manufacturer’s protocol at the National Instrumentation Center for Environment Management (Seoul, South Korea). For PCR-based validation of several mutations of interest, additional RR1 and HB101 cells were purchased from KCTC, Korean Culture Center of Microorganisms (Seoul, South Korea), and TaKaRa Bio (Kusatsu, Shiga, Japan): KCTC 1473 (= ATCC 31343; RR1), KCTC 1467 (= ATCC 33694; HB101), KCCM 70032 (= ATCC 33694; HB101), and HB101 competent cells (TaKaRa cat. No. 9051). Culture condition and DNA isolation method were all the same for the other strains.

Because it is well-known that E. coli HB101 and its descendant RR1 had a K-12:B hybrid genome structure, we first constructed a “backbone” genome sequence by in silico recombination of K-12 and B genome sequences, where parental regions were determined from the initial read mapping on K-12 MG1655 and BL21(DE3) sequences separately. Subsequent read mapping and sequence correction were carried out, and finally, residual regions that could not be resolved were validated using PCR amplification and Sanger sequencing. Pretreatment of reads (quality limit 0.01, maximum allowed ambiguous base 1, and minimum read length of 50), reference mapping, and subsequent sequence manipulation were carried out using the CLC Genomics Workbench version 6.5.1 (Aarhus, Denmark). Only paired reads passing the pretreatment step were mapped to the reference sequences of E. coli K-12 MG1655 (NC_000913.3) and E. coli B BL21(DE3) (NC_012971.2) separately, and quality-based variant calling was run. Using putative recombinational junctions inferred from the distribution of single nucleotide variants (SNVs), a hybrid genomic sequence containing part of the B strain genome (190.4 kb) was constructed in the K-12 genomic backbone. Based on this hybrid sequence as the starting reference, a series of sequence manipulation steps consisting of (i) read mapping, (ii) variant detection, and (iii) consensus sequence extraction were iterated until no further mutations were detected. Finally, 25 breakpoints (analyzed by CLC Genomics Workbench) selected on the basis of read number, p-values, and fraction non-perfectly mapped, as well as one K-12 island in the B genome background, were grouped in 11 clusters based on genomic location and then corrected using PCR and Sanger sequencing of the amplified products. The Phred/Phrap/Consed package¹ was used for the final sequence manipulation. Primer sequences and their information are shown in Table 1. Amplification targets for primer pairs starting with ‘P’ are shown in Table 2. For validation of the mrr-hsdRMS-mcrBC deletion, two primer pairs, L-outer:L-inner and R-inner:R-outer, were designed. Because the two inner primers were designed on the basis of the E. coli BL21(DE3) genome sequence that does not have deletion, these two primer pairs encompassing deletion junctions would produce amplification products from genomes that have an intact mrr-hsdR::IS1-hsdMS-mcrBC region (967 and 867 bp, respectively). In case of deletion, only the outer primer pair (L-outer:R-outer) would produce a 919-bp product. Primer pairs mlc_F:mcl_R and recA_F:recA_R were used to check mutations in mlc and recA genes, respectively. Primer pair rhsA-L:rhsA-R and the internal sequencing primer rhsA-I were used for the validation of the rhsA sequence. Final validation of the reconstructed RR1 genome sequence was carried out using CLC Genomics Workbench-based re-mapping of Illumina reads and breseq version 0.27.1 (Deatherage and Barrick, 2014). Genome annotation was carried out by the RAST server and NCBI Prokaryotic Genome Annotation Pipeline.

The genealogy from wild-type E. coli K-12 (F⁺ λ ⁺) to HB101, a hybrid strain of E. coli K-12 and B, was traced by literature search (Figure 1). The pedigrees from wild-type E. coli K-12 to Y10 (F^- thr-1 leu-6 thi-1 supE44) (Lederberg and Tatum, 1946; Lederberg, 1947, 1952; Tatum and Lederberg, 1947), and from Y10 to W2961 [= AB266; F^- araC14 leuB6(Am) Δ(gpt-proA)62 lacY1 glnX44(AS) galK2(Oc) λ - Rac-0 rfbC1 mgl-51 rpsL20(Str^R) xylA5 mtl-1 thiE1] (Sypherd, 1965) are well-documented in Bachmann’s work (Bachmann, 1972, 1996). The earlier auxotrophic mutations up to Y10 were introduced by X-ray irradiation, whereas the mutations in latter generations leading to W2961 were introduced by repetitive UV irradiation and subsequent use of various selective media. The generations from W2961 via HB101 (RecA^-) to RR1 (RecA⁺) were all produced using F′- or Hfr-mediated conjugation. However, HB101 has been incorrectly described in a literature (Singer and Berg, 1991) as being derived from RR1 by mutation of the recA gene.

HB101 is one of the restriction-modification deficient mutants (r^-_B m^-_B) produced by the pioneering work of Boyer and Roulland-Dussoix (1969) to investigate the genetic basis of the restriction-modification system. HB101 was the first mutant to indicate the presence of the third cistron (ramC, currently known as hsdS) that comprises the restriction-modification system. In their study, r⁺_B m⁺_B alleles and thr⁺ from the B strain were simultaneously introduced into the genome of AB266 to make HB16. The recA13 mutation was introduced into the HB16 strain by conjugation with the HB82 Hfr strain (= AB3045) to make HB100. Subsequently, different r^-_B m^-_B mutants were segregated from a cross between HB100 and HB77.

The entire evolutionary process from AB266 to HB101, including details of all participant strains, was further elaborated by Rothen (1997). For example, while the B donor strain for r⁺_B m⁺_B alleles was not specified in the original report (Boyer and Roulland-Dussoix, 1969), Rothen reported that AC2517 (B/r F^-) (Boyer, 1964) was converted to a conjugal donor F′ strain (AC2516) (Boyer, 1964) after crossing with AC2515 (F′-lac⁺) (Johnson et al., 1964), and that HB16 was produced by the cross between AC2516 and AB266, not by P1 production. By contrast, Boyer and Roulland-Dussoix stated that r⁺_B m⁺_B alleles were “co-transduced” with thr⁺ genes. The B genomic fragment (190.4 kb; including the ∼20 kb region that was later deleted by IS1-mediated recombination; see below) harbored in the genome of RR1 is too large to be introduced by a single phage P1 transduction event. Therefore, the term “co-transduced” seems to be used as a generic term describing the introduction of foreign DNA.

Although the RR1 strain is frequently mentioned in the scientific literature, and even in online lists of laboratory E. coli strains, we could not discover for what purpose RR1 was constructed from HB101. Raymond L. Rodriguez explained via personal communication as follows: RR1 was constructed by conjugal mating of HB101 (recA13) with an F⁺ RecA⁺ strain with a view to using it as a host for colicin-resistant plasmids that work best in a RecA⁺ background. We found that the genomic position (2,772,443–2,806,787) including recA (SR35_13775) in the RR1 strain was identical to the intact form found in K-12 MG1655 (RecA⁺). Sanger sequencing of PCR products amplified from recA loci of several strains confirmed that only HB101 strains have a mutant recA13 allele (Leu⁵²→Phe).

We produced a total of 48,255,170 paired reads (4.97 Gb) from the library with an average insert size of 359 bp. Of these reads, 44,074,332 (4.24 Gb; 92.9% paired) passed quality trimming and filtering, and were used for mapping. SNV distribution shows that the finalized RR1 genome is a hybrid of K-12 (∼95.9%) and B (∼4.1%) (Figure 2); these proportions were later revised taking indels into account. A complete list of SNVs identified from the first round of read mapping was given in Supplementary Table S1. Given that B genomic DNA was introduced by a single event, two recombinational junctions where crossover occurred were defined between the last SNV identified in one genome reference and the first SNV in the other reference. In fact, the two junctions are short stretches of nucleotides in which K-12 MG1655 and BL21(DE3) are identical with each other, such that homologous recombination can occur. When K-12 MG1655 was used as the reference, 270 SNVs were densely distributed in a narrow genome range (∼27 kb; 302–27,367), where the last nucleotide position corresponds to 31,439 in the BL21(DE3) genome. On the other hand, when BL21(DE3) was used as the reference, the first SNV occurred at 31,505, corresponding to 27,433 in the K-12 MG1655 genome. When the same approach was applied to the other recombinational junction, we could construct the first hybrid sequence (4,651,433 bp) consisting of 1–31,504 from BL21(DE3), 27,433–4,488,508 from K-12 MG1655, and 4,400,100–4,558,953 from BL21(DE3).

Using this sequence as the starting reference, three rounds of mapping, variant calling, and consensus extraction were performed. During the initial mapping step, five zero-coverage regions that could be due to large-scale deletions were identified from visual inspection of aligned reads (Supplementary Figure S1). With iterative mapping and reference correction, discrepancies between reads and the reference sequence were gradually mitigated and finally, confirmatory read mapping revealed that a correct genome sequence with regard to large deletions could be reconstructed without the need for PCR and Sanger sequencing.

Large-scale insertions, however, could not be accurately reconstructed by mapping reads and revising reference sequence only. Therefore, candidate regions identified by manual inspection of read depth and by breakpoint analysis were amplified by PCR and confirmed by Sanger sequencing. Large-scale insertions and deletions, mostly results of prophage- or insertion sequence (IS)-mediated events, are summarized in Table 2. Meanwhile, a moderate length deletion (123 bp) resulting in a truncation at the C-terminus of yghQ (SR35_15225) (data not shown) did not appear as a conspicuous zero-coverage region, leaving a stretch of 66 Ns after the final consensus sequence extraction. This region was corrected by manual in silico extension of reference sequence using the unaligned ends of partially mapped reads at the position followed by sequence joining.

The length of finalized genome sequence of RR1 is 4,587,291 bp with 50.8% G+C. If recombinational junctions are defined as a and b (see Figure 2B), we can say that 31,349 ≤a < 31,505 and 4,447,909 < b ≤ 4,451,871 in RR1 genome coordinates, which means the length of extant B genomic DNA introduced by transduction is 166.8–170.9 kb. A single nucleotide C to T variation at 4,448,417 [4,489,016 and 4,400,608 in K-12 MG1655 and BL21(DE3), respectively; within yjgR of the right recombinational junction] was identified whichever references were used, implying this mutation occurred after the introduction of the B genomic fragment.

During visual investigation of aligned reads on the B reference genome, an SNV-rich segment that was almost identical to K-12 MG1655 (99.97% identical; only 1-bp difference out of 3521 bp) sequence was found at position 4,577,076–4,580,596 and was verified using PCR and Sanger sequencing. The segment corresponded to 4,548,791–4,552,258 in B genome coordinates and contained slt-trpR-yjjX-ytjC (gpmB)-rob (SR35_22630-SR35_22655; Supplementary Figure S2). Out of 48 SNVs occurring in the B genome segment, 45 were concentrated in this narrow 3.5 kb region (Supplementary Table S1). They were ∼136 and ∼38 kb apart from each end of recombinational junctions, respectively. We confirmed the presence of this small K-12 island in all five strains chosen for this study (two RR1 strains and three HB101 strains) using PCR and Sanger sequencing. We also encountered an opposite situation in the K-12 genomic background of RR1 strain, where B-like 17 SNVs were concentrated in a short 2.1-kb region (yaaU-kefF-kefC; SR35_00255-SR35_00265; Supplementary Table S1). Because this region is only ∼18 kb apart from the last recombinational junction (rihC; SR35_00165), the integration at the secondary site might have occurred concomitantly with the primary integration event. Such dispersed integration was discussed in the previous report (Studier et al., 2009). The presence of flanking sequences (>3 kb both) where K-12 and B genome sequences are identical with each other (Supplementary Figure S3) might have facilitated the integration of B-like “islet” in the K-12 genome region.

Special care had to be taken when resolving the boundaries for invertible P-DNA segment (GenBank X01805.1) (Plasterk and van de Putte, 1985), a 1,797 bp long element flanked by two inverted repeat sequences (5′-TTGGTTTGGGAGAAGG-3′) within the cryptic prophage e14. New “junction sites,” identified by the structural variant detection function of CLC Genomics Workbench and breseq during the validation process of the final genome sequence, suggested that the sequencing library is mixed such that roughly half of the library molecules have inversions. PCR results of junction regions also suggested that the genomic DNA used for sequencing library construction had both orientations of the internal sequence element (data not shown). When reads were mapped with a high-stringency condition (match score 1, mismatch cost 5, length fraction 0.99, and similarity fraction 0.99) simultaneously to the two P-DNA segments in both directions having 1 kb flanking sequences beyond terminal repeats, average read coverage was 460.17 (normal direction in compliance with NC_000913.3) and 511.73 (inverted orientation). Even though the cell culture for genomic DNA preparation was inoculated from a single colony, we observed the dynamic inversion of the P-segment in the single strain population during a short-term period. A BLAST search against the nucleotide collection at NCBI showed that there are at least five genome sequences having a P-DNA segment with inverted orientation besides X01805.1. CLC Genomics Workbench and breseq also suggested a low coverage region at the 3′-end of rhsA gene, but deletion was not detected except for correcting two nucleotides using PCR and Sanger sequencing (data not shown).

Most of the genes in the genomic fragment derived from the B strain are syntenic with respect to their K-12 counterparts, but some of them have different amino acid sequences, indicating that they may have functional differences with respect to the parental K-12 strain. The presence or absence of genes specific to each strain confers phenotypical differences. For example, SR35_00100–SR35_00110 encoding fimbrial proteins and type III effector-like protein (between nhaR and IS1) is present only in B, a choline transporter downstream of the fec cluster (SR35_22010) is present only in B, the yjhIHGFU cluster is present only in K-12, the nanCMS cluster for the utilization of sialic acid is present only in K-12, fimB (SR35_22195) is interrupted by IS1 in B but is intact in K-12, a yjiV-hpaB (C-terminal) deletion encompassing the entire mrr-hsdRMS-mcrBC locus is present only in RR1, the hpa cluster for 4-hydroxyphenylacetic acid catabolite pathway is present only in B, and yjjJ is present only in K-12.

The immigration control region (ICR), consisting of the mrr-hsdRMS-mcrBC gene cluster, together with the O-antigen region, is the most divergent region between E. coli strains (Milkman et al., 2003). The hpa cluster is downstream of the ICR. Because the ICR and neighboring genes are completely absent in the RR1 genome (see below), the KpLE2 prophage-like region (SR35_21940–SR35_22160) in the B genomic fragment is the most divergent compared with its cognate K-12 region. The aforementioned choline transporter gene and yjhIHGFU cluster are all located in the KpLE2 prophage region. The presence of ISs was one of the main factors contributing to inter-strain sequence variation in this prophage region. In K-12 MG1655, there are single copies of IS2, IS4, IS911 (interrupted by the following IS30), IS30 and IS1, but RR1 has only IS911 (uninterrupted; SR35_22030) and IS1 (SR35_22080-SR35_22090) surrounding fecIRABCDE operon.

The complete genome sequence of RR1 can help us to elucidate all known and yet-to-be-discovered characteristics of RR1 (Table 3), but also tell us more about its history. We should be aware that genomic differences between RR1 and K-12 MG1655 are in fact the sum of mutations occurring in each descendant from the common ancestor, wild-type E. coli K-12 (Bachmann, 1972, 1996). Strain MG1655 (F^- λ ^- rph-1) was rapidly obtained from wild-type K-12 via W1485 (F⁺ λ ^- rph-1) after only one round of UV irradiation and acridine orange mutagenesis (Bachmann, 1996). Compared to genomic differences between E. coli B strains REL606 and BL21(DE3) (Jeong et al., 2009), most of which were caused by disparate integration of K-12 DNA in a narrow region of the genome (Studier et al., 2009), the 258 SNVs between RR1 and MG1655 constitute a much larger difference in the same K-12 genomic background. The extensive mutagenesis by X-ray and UV radiation of the RR1 lineage (Figure 1) might account for this observation.

As shown in Table 3, the characteristics of RR1 and HB101 strains reported in the literature or on websites show some discrepancies. For example, Bolivar et al. (1977) mistakenly described RR1 as r^-_K m^-_K, whereas its parental strain HB101 had been clearly described as r^-_B m^-_B. The most significant discrepancy concerns Δ(mcrC-mrr) or the mcrC mrr genotype. Our sequencing results demonstrated an IS-mediated deletion (∼20 kb) from the yjiV gene at the 3′-end of the hpaB gene (Supplementary Figure S1). The deleted region includes mrr-hsdRMS-mcrBC, whose deletion results in more permissive host strain than that of strains harboring point mutations (Woodcock et al., 1989; Grant et al., 1990; Doherty et al., 1991). We assumed that the mcrC-mrr deletion leading to the r^-_B m^-_B phenotype did not occur when HB101 was made from HB100 (r⁺_B m⁺_B). Reportedly, HB101 is trans dominant to the r⁺ phenotype, as shown by complementation analysis (Boyer and Roulland-Dussoix, 1969), which means that the original HB101 strain produces mutant HsdS proteins (caused by hsdS20 mutation) that interfere with wild-type restriction function (Arber and Linn, 1969). Therefore, our sequencing results showed that IS-mediated deletion of mrr-hsdRMS-mcrBC must have occurred in a later generation of the original HB101 cells harboring the hsdS20 mutation. We checked the prevalence of the mrr-hsdRMS-mcrBC deletion in the five HB101 strains using PCR and Sanger sequencing. All tested HB101 strains harbored this deletion, implying that it might have occurred immediately after the construction of HB101 and might not be confined to the RR1 lineage.

The genotype of wild-type K-12, which comprises most of the RR1 genome, differs slightly from source to source. Bachmann (1972, 1996) stated that wild-type K-12 is simply F⁺ λ ⁺ in the pedigree figures omitting rfb-50 or Δrfb-51 mutations common to K-12 wild-types, but mentioned these mutations in the main text of references. Hayashi et al. (2006) stated that, without mentioning rfb mutations, wild-type E. coli K-12 has rpoS(Am) and rph-1 (frameshift) mutations, while K-12 MG1655 has the pseudorevertant allele (Q33) at the rpoS locus. However, in Bachmann’s pedigree, the rph-1 mutation appears in W1485 for the first time, not in its parent wild-type K-12. The genotype of K-12 MG1655 (F^- λ ^- ilvG^- rfb-50 rph-1) is in good agreement with current genome sequence information of K-12 MG1655. The genotype inconsistency implies either that mutations found later have not yet been added to the list of wild-type characteristics, or that different culture stocks of the same strain have independent variations (Freddolino et al., 2012). Alternatively, it may arise from misinterpretation of experimental data or literature search results.

Notably, the rfbD1 mutation in the dTDP-4-dehydrorhamnose reductase gene is absent from currently available genotypes of HB101, RR1, Bachmann’s first pedigrees of K-12, and their mutant derivatives (Bachmann, 1972). However, later work (Bachmann, 1996) stated that rfbD1 is present in Y10 and its direct descendants. We found a frameshift mutation in the rfbD gene (SR35_10405) in the genome of the RR1 strain, which is identical to that found in recently sequenced E. coli ER1821R (Jobling et al., 2016), a K-12 derivative laboratory strain harboring the ancestral rfbD1 mutation. An IS5 insertion (SR35_10355), designated rfb-50, that is common in K-12 wild-type strains, was also found in the RR1 strain. Our complete genome sequence can explain the genotype or mutations of the RR1 strain (Table 4).

Whereas an amber mutation in rpoS was found in RR1 (SR35_13990), no mutation in the rph gene was observed. This implies that the rph-1 mutation might have occurred in the descendant of wild-type K-12, contrary to a previous report (Hayashi et al., 2006). Because Jensen (1993) observed the rph mutation only in MG1655, W3110, and their common ancestor W1485, the original wild-type strain may have a normal rph allele. Unexpectedly, we observed an amber mutation (C1105T) in the open reading frame of the mlc gene (also known as dgsA) in the RR1 strain, which encodes a global regulator of carbohydrate uptake including the glucose phosphotransferase system. Although it is not known how the RR1 strain obtained the amber mutation, DNA sequencing of the mlc gene of HB101 strains revealed the same mutation, indicating that disruption of the mlc gene occurred before the construction of the RR1 strain.