Re-sequencing and optical mapping reveals misassemblies and real inversions on Corynebacterium pseudotuberculosis genomes

The number of draft genomes deposited in Genbank from the National Center for Biotechnology Information (NCBI) is higher than the complete ones. Draft genomes are assemblies that contain fragments of misassembled regions (gaps). Such draft genomes present a hindrance to the complete understanding of the biology and evolution of the organism since they lack genomic information. To overcome this problem, strategies to improve the assembly process are developed continuously. Also, the greatest challenge to the assembly progress is the presence of repetitive DNA regions. This article highlights the use of optical mapping, to detect and correct assembly errors in Corynebacterium pseudotuberculosis. We also demonstrate that choosing a reference genome should be done with caution to avoid assembly errors and loss of genetic information.

In order to solve this assembly problem and to improve the generated data, we have the strategy of optical mapping, or Whole Genome Mapping (WGM), which is an approach that uses high-resolution restriction maps to generate the actual orientation of the organism's genome. It is the main method of large-scale genome analysis that provides complete visualization of the structural genome through a single image 12 . Optical mapping is based on the distance of the restriction sites for high precision map construction. This is a strategy in which data are obtained with greater precision since it is a physical result of the genome evaluation. This method, combined with the application of de novo assembly methodology, assists in the orientation of contigs 13 .
The technique of optical mapping was first developed by Schwartz and collaborators in 1993, with the purpose of studying the chromosomal gene ordering of Saccharomyces cerevisiae 13 . Samad et al., 1995, describe optical mapping as the novel approach for single-molecule DNA analysis using flowering microscopy to identify and estimate its size by the generated images 14 . Since then, several improvements have been added to the technique, especially in the images and algorithms for fragment size estimation 15 . Hence, WGM gained notoriety in several applications, such as in lineage typing in epidemic cases for clinical microbiology 16,17 ; ordering of contigs generated by de novo assembly 7 ; and in the study of inversions, insertions, deletions, duplications, and instability of bacterial genomes 18,19 . WGM has been successfully performed on very different types of organisms such as bacteria [20][21][22] , fungi 23 , plants [24][25][26] , and mammals 27 .
Regarding the genomic assembly strategy, WGM is an additional method that allows the ordering of the contigs and thus provides a size estimation of the gaps and their positions. This combination of methods is called a hybrid approach to scaffolding assembly, and it is feasible to acquire a complete genome of high quality and accuracy 7 .
The strategy is to make C. pseudotuberculosis the most used organism in genomic studies involving the Corynebacterium genus. Therefore, a total of 11 strains were selected from different hosts, isolation sites and distributed between ovis and equi biovars, so that complete genomes can be made available, well assembled, and updated by new sequencing. In this manner, this data can be explored with greater reliability by future comparative studies, intraspecific evolutionary relationship analyses.

Results
Sequencing and assembly. The strains deposited by our research group were either re-sequenced (i.e., Cp1002 (1002B), CpI19, Cp31, Cp162, Cp258, CpCIP52.97), or were first sequenced using the Ion Torrent PGM TM platform (i.e., Cp29156, CpMB302, CpMB11, CpT1) ( Table 1). Different software packages were used for de novo assembly (Table 1). optical mapping analysis: biovar ovis strains. The strains Cp1002 (CP001809.2) and CpI19 (CP002251.1) (Fig. 1A,B) showed an inversion of approximately 1.6 Mb and 1.22 Mb, respectively. It is observed in the regions flanking the first and third clusters of Ribosomal RNA in the CpI19 strain; while in Cp1002, it occurs between the first and fourth clusters. Figure 1A,B show the starting and ending points of the inversion, labeled as R1 and R2, respectively. The central block in Fig. 1 corresponds to the physical restriction map, while the upper and lower blocks represent the in silico restriction map generated by MapSolver TM software. The red regions of the central block (Fig. 1A,B) indicate that the same region exists in both the first version (Upper Block) and the updated version (Lower Block). Thus, they show that there was no significant loss between the compared versions in that region. Strain CpFRC41 (CP002097.1) (Fig. 1D) showed a correct alignment of the restriction sites in the whole genome and thus, with no probable errors of assembly and ordering. Strains Cp29156 and CpT1 (Fig. 1C,E) were first sequenced in this work. optical mapping analysis: biovar equi strains. Figure 2A shows the alignment between the first (upper map) and the last (lower map) versions of the CpCp31 strain. The R1 region highlights the absence of corresponding restriction sites at the beginning of the genome. R2 region, in its turn, shows the absence of a chromosome region. This difference probably occurred due to errors in the assembly and gap closure process.
Worse problems were found in the previous version of Cp162 (Fig. 2B). The R1-labeled region, starting in 5′ end of the dnaA gene, shows no similarity between the restriction site patterns. In the R2 region, there is no linear alignment between the sites, probably due to the absence of genes. According to the optical map, an error may have occurred in the ordering of contigs in the R3 region, where the segment should be in another region of the genome. R4 region shows a ~0.85 Mb inversion in the middle of the genome, located explicitly between two clusters of ribosomal RNA.
The CpMB11 strain presented an error of choice of chromosome initiation site, in which a segment close to 5′ end of the dnaA gene should be situated at the end (Fig. 2C, regions R1 and R2). The Cp258 strain did not present chromosomal inversions in the deposited genome. However, a misalignment of the restriction sites (Fig. 2D, region R1) is shown next to the origin of replication at the 5′ end of the genome. Also, a region containing the genes moaE, molB, molA, narI, narJ, narH, narG, narK, narT, moeY, moaC is absent in the deposited chromosome ( Fig. 2D, R2). The same situation occurred with the strain CpCIP5297 (Fig. 2E, regions R1, and R2). Finally, as for the Cp29156 and CpT1 strains, it is the first sequencing of the CpMB302 strain (Fig. 2F). www.nature.com/scientificreports www.nature.com/scientificreports/ content and genomes plasticity. This analysis showed a reduction of 136 bp in the total genome of the CpI19 strain, in addition to an increase of 34 CDSs and a reduction of 12 pseudogenes ( Table 2). The updated version of strain Cp1002 28 showed a reduction of 6 bases in the total genome and a reduction in the number of annotated CDSs (n = 24) when compared with the older one. In this case, updating the annotation of genes identified as hypothetical protein might be the explanation. This comparison can be visualized in the map generated by BRIG, in which the last version is compared with the most updated one before the optical mapping. Strains The equi strains showed more changes among the genomes. In the Cp258 genome, obtained from the Ion Torrent PGM sequencing, an insert of ~55 kb was added. A total of 41 new CDSs and a reduction of 12 pseudogenes were included (Table 2). This had not been represented within the first genome, which was obtained by SOLiD v3 sequencing. Differences can also be visualized on the genomic map (Fig. 4A), where the highlighted red genes are important and essential genes for the strain classification e.g., the operon nar, with the moaE, moaD, molB, molA, narI, narL, narG, narK, narT, moEY, mobA, moAc, moeA1 genes. Presence of genes coding for proteins such as Beta-lactamase, Vitamin K-dependent gamma-carboxylase, Heavy-metal-binding protein, Transposases, Type I restriction-modification system, and N-6 DNA Methylase is also important. Moreover, errors related to positioning and presence of CRISPR associated proteins were found among the assemblies.
The strain Cp162 also presented an increase of genomic content (~72 kb). Ninety-seven CDS and a reduction of 44 pseudogenes were found ( Table 2). Figure 4B shows the absent regions of genes such as the complete cluster of the operon nar. Genes coding for Fe 3+ dicitrate transport, ATP-binding protein FecE, Beta-lactamase, Vitamin K-dependent, gamma-carboxylase, Heavy-metal-binding protein, Phytoene dehydrogenase, CAAX protease self-immunity, Restriction endonuclease or methylase, Collagen-binding surface protein Cna-like, B-type domain protein, Membrane protein, ATP-dependent exonuclease, and several hypothetical proteins were also absent. A possible assembly error in the cluster of genes coding for CRISPR-associated proteins was found.
The new sequencing by Ion Torrent of the strain Cp31 resulted in the most significant increase in the gene content (~106 kb) among the selected strains. An increase of 110 CDSs and a reduction of 42 pseudogenes (Table 2) were detected. In Fig. 4C, we can highlight the inclusion of the corynephage with the tox gene of diphtheria toxin.  The CpCIP52.97 strain had an increase of ~49 kb, which represented a gain of 127 CDSs and the reduction of 13 pseudogenes (Table 2). In Fig. 4D, we can highlight the absence of essential genes, such as genes associated with CRISP (cas2, cas1, cas3, cas4, cas5, cas6, cas7). Once more, the operon nar was absent, as well as genes coding for Beta-lactamase, Phytoene dehydrogenase, integrins, transposases, and several hypothetical proteins.
With the complete and finalized genomes, a multiple alignment analysis of the genome was done using Mauve. Figure 5 shows the genomes of 5 ovis strains (i.e., Cp1002B, Cp29156, CpFRC41, CpI19, and CpT1). Blocks with the same color represent conserved regions, in which they share a high similarity. Inversions and rearrangement events were established as changes in the Cp1002B reference synteny (first genome of the figure). By analyzing the ends of the inverted blocks, the inversions are flanked by clusters of rRNAs. Only the CpFRC41 and CpT1 strains showed the same gene order (synteny) in their genome. Using the same strategy to compare the equi strains (e.g., Cp31, Cp258, Cp162, CpCIP52.97, CpMB302, and MB11), the strain Cp162 was the only one that showed an inversion and rearrangements (Fig. 6). However, these blocks are not flanked by the rRNA clusters.

Discussion
Although the optical map initially focused on typing and identification of strains without the need of sequencing, it is an efficient approach to order contigs generated by de novo assembly, allowing error detection and an accurate contig ordering in assemblies 30 . By definition, the optical mapping technique is a single molecule barcode using restriction enzymes, and the distances among these restriction sites, which are the basis for the alignments between optical maps and the in-silico maps (by contigs or genomes). These features are excellent for the scaffolding process with de novo assembled contigs 31 . Onmus-Leone and collaborators in 2013 applied this technique to successfully order and correct contigs generated in the de novo assembly step using pyrosequencing data of the bacteria Providencia stuartii. The mentioned authors established the strategy of using the optical map to order and correct contigs generated from short reads 32 . www.nature.com/scientificreports www.nature.com/scientificreports/ Latreille and collaborators in 2007 highlighted the efficiency of the optical map suggesting it as a routine procedure in the assembly finishing process using de novo assembly 33 . According to these authors, it is possible to detect errors in the order and construction of the contigs by using this technique even when the organism presents several repeated regions. In the mentioned work, it was possible to finish the assembly by using optical map data even when cosmid libraries and overlapping restriction maps of BACs 33 have already been applied without success. Analyzing these results, it was possible to conclude that the optical maps are an excellent option for bacterial assembly finishing because the restriction sites cover these repetitive regions in most cases.
Repetitive regions are considered a major difficulty for de novo assembly algorithms, mostly in transposons and ribosomal RNA cluster regions 33 . The strain Cp1002 was sequenced by using 454 Genome Sequencer FLX (Roche), Sanger e PacBio technologies, but the error of inversion in rRNA clusters continued. Only when  www.nature.com/scientificreports www.nature.com/scientificreports/ sequencing with Ion Torrent PGM TM and ordering the contigs by using optical mapping data were done, the genome was correctly finished. The strain CpI19 was, at first, sequenced by using SOLiD v2 technology with 25 bp mate-paired reads and a coverage depth of 321-fold. The mate-paired technology may have contributed to the correct construction of the contigs, but due to its small size of reads sequenced, the inversion of the ribosomal RNA cluster regions occurred.
Trost and collaborators sequenced CpFRC41 strain in 2010 34 by using 454 Genome Sequencer FLX (Roche) sequencer; it was the first C. pseudotuberculosis strain deposited in GenBank by NCBI in 2010. The assembly was performed using Corynebacterium diphtheriae NCTC 13129 (BX248353.1) as the reference genome, and the gaps were closed by using Polymerase Chain Reaction (PCR) and the software r2cat 35 . Gap filling by using PCR probably contributed to no inversions found in the final genome (Fig. 1D). Schröder and collaborators in 2011 successfully used this approach in Corynebacterium variabile DSM 44702 36 .
In CpMB11 strain (Fig. 1E) a standardization problem regarding the region of the 5′ end of the dnaA gene was found. Several works have shown that the most common origin of replication in bacteria is oriC and the first gene is the chromosomal replication initiator protein DnaA (dnaA) 37 . Thus, it is used as a standard pattern to linearize bacterial genomes.
Cp31 have been sequenced using several platforms: Solid v2 (CP003421.1), Ion Torrent PGM TM (CP003421.3) and PacBio. This fact confers more reliability to the assembly of this strain, which leads us to consider it as the reference strain in equi biovar. In its last version, published by Viana and collaborators in 2017, it was assembled using optical mapping technology. The strains Cp258, Cp162, Cp31 e CpCIP52.97 were initially sequenced in the SOLiD platform; these strains also belong to equi biovar and were characterized by using biochemical tests. Those were re-sequenced using Ion Torrent PGM TM platform and novel genomic regions were added to the genome sequence. Essential genes in regions, such as nar cluster and clusters associated with CRISPR proteins, which are only present in equi biovar, were added. The missing regions may be caused by an error propagation due to reference contig ordering, because of the first equi genome available, the Cp258 38 strain, was assembled using an ovis biovar strain as the reference, which does not completely present the referred clusters. The same issue happened to Cp31 39 isolate.
Similarly, the Cp162 40 strain presented the same issue because its contigs have been ordered according to C. pseudotuberculosis 316 (Cp316) (CP003077.1) 41 . Cp316 strain belongs to equi biovar and does not present nar operon in its former genomic sequence; it was assembled using the ovis strain CpFRC41 as a reference genome for contigs ordering. Presumably, the same problem happened in CpCIP52.97 42 , but the genome used as a reference in the assembly process is not described in the article.
Using a complete genome deposited in public databases as a reference to assemble novel genomes is a risky strategy because even if it generates complete genomes more efficiently, it may disseminate assembly errors from one strain to other 8 . In this article, the optical map is used to validate contigs, and it is shown that extension and gap filling using read mapping or de novo assembly may generate assembly errors. We highlighted this technique because it does not present sequencing bias. Assembly statistics such as N50, coverage and depth coverage may generate false positive answers. Another strategy suggested by Lehri and collaborators is to use non-paired reads together with paired or long reads. A genomic region presenting assembly errors caused by insertions, deletions, inversions or rearrangements may hide significant biological variations or produce false interpretations, mostly in genomic analysis 43 . Even before the NGS platforms boomed, Schmutz and collaborators (2004) were concerned about the quality of the human genome 43 , mostly because of the possible assembly errors.
The inversions caused by RNA ribosomal clusters in ovis biovar strains may occur due to the high similarity of these clusters. This kind of inversion has already been shown in literature in bacteria as Salmonella paratyphi A using restriction enzymes and pulsed-field electrophoresis gel (PFGE) 44 . The inversions among strains of the same species may be comprehended as homologous recombination events 45 . It can be highlighted that these inversions do not occur in strains belonging to equi biovar, except for Cp162. The inferred data concerning the genomic order of C. pseudotuberculosis strains were only achieved because optical mapping technology provides an accurate in vitro evidence.

Methods
Strain and DnA isolation and Genome sequencing. The methodology described below was applied to the novel sequencings of the strains C. pseudotuberculosis 1002(1002B), C. pseudotuberculosis 29156, C. pseudotuberculosis I19, C. pseudotuberculosis Cp162, C. pseudotuberculosis 258, C. pseudotuberculosis 31, C. pseudotuberculosis MB302, C. pseudotuberculosis MB11, C. pseudotuberculosis T1 and C. pseudotuberculosis CIP5297. The strains were cultivated in solid media with 1.5% of bacteriological agar. Subsequently, an isolated colony was used to grow in liquid media with brain-heart-infusion media (BHI-Hi Media Laboratories Pvt. Ltd, India) supplemented with 0.5% of Tween 80, at 37 °C for 20 hours under rotation. Genomic DNA was extracted following the protocol of Pacheco in 2006 46 . After the extraction step, the libraries were constructed with IonXpress ™ Plus DNA Fragment Library Preparation Kit. The DNA samples were fragmented using Ion Shear TM Plus for five minutes at 37°. Then, adaptors from Ion Xpress ™ Barcode Adapters kit were ligated for library quantification. Subsequently, the fragments were amplified using Ion PGM ™ 200 bp or 400 bp kits. These reactions were transferred to the semiconductor chip (ION 318 TM Chip v2), and it was put into Ion PGM ™ . During all the steps described above, all the manufacturer's instructions were strictly followed. No novel sequencing was performed for C. pseudotuberculosis FRC41.
Genome assembly and annotation. The analysis of the reads quality was performed by using FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). No trimming was performed on reads with Phred score above 20, which were the majority. For contigs construction, we applied the de novo assembly strategy (no reference used) by using Mira v. 3.9.18 46  www.nature.com/scientificreports www.nature.com/scientificreports/ software. Scaffolds construction was manually performed in CLC Genomics Workbench (CLC-gw) version 7.0 (Qiagen) software using the visualization of the contigs mapped and ordered according to the restriction sites of the strains in MapSolver TM (OpGen). Then, dnaA gene was fixed in the probable oriC position in the chromosome by using an in-house python script. In order to fill the gaps and finish the assembly, GapBlaster 49 and FGAP 50 software was used and subsequently, the contigs were mapped to the scaffold or a reference genome by using CLC Genomics (Qiagen). The annotation was performed by using in-house scripts for annotation transference from C. pseudotuberculosis strains, which were manually curated in the UniProt database (http://uniprot.org). Finally, pseudogenes were manually curated by using Artemis software 51 and CLC Genomics (Qiagen). optical mapping. The optical maps were acquired from Opgen, Inc. The MapSolver ™ (OpGen Inc.) software was used for the comparison of the physical restriction map and the restriction sites present in the assembled genome. Several pieces of information about metrics for the quality of each optical map are available in Table 3.
Genome plasticity and genetics content. This analysis was performed using genomic sequences before and after assembly assisted by optical mapping data. The maps comparing different versions of the studied strains were generated by using Blast Ring Image Generator (BRIG) v0.95 52 . For inversion, deletion and rearrangement analysis, the Mauve v. 2.3.1 53 software was used with progressiveMauve 53 option set.

conclusion
The results obtained from optical mapping data analysis pointed errors in the assemblies of C. pseudotuberculosis genomes deposited in Genbank. Thus, the optical map was efficient in the assembly error detection of the strains Cp1002, CpI19, Cp31, Cp162, CpMB11, Cp258, and CpCIP52.97. Regarding the novel genomes, such as Cp29156, CpT1, and MB302, the optical map data contributed in the contigs ordering step, which contributed to a more efficient assembly finishing considering there are no assembly errors in the final version of the genomes. Furthermore, the update of the genomic sequences done by re-sequencing the genomes with Ion Torrent PGM TM platform was essential to the relevant genomic content increase, which happened in the strains previously sequenced using the SOLiD platform. We also pointed out several inversions caused by ribosomal RNA gene clusters in strains of the ovis biovar. Thus, we can suggest that the genomes deposited after applying this strategy made these strains more reliable for novel studies.

Strains
Enzyme Length (bp)  Table 3. Information about the quality metrics of the optical maps used.