Molecular characterization of an unauthorized genetically modified Bacillus subtilis production strain identified in a vitamin B2 feed additive.

Many food and feed additives result from fermentation of genetically modified (GM) microorganisms. For vitamin B2 (riboflavin), GM Bacillus subtilis production strains have been developed and are often used. The presence of neither the GM strain nor its recombinant DNA is allowed for fermentation products placed on the EU market as food or feed additive. A vitamin B2 product (80% feed grade) imported from China was analysed. Viable B. subtilis cells were identified and DNAs of two bacterial isolates (LHL and LGL) were subjected to three whole genome sequencing (WGS) runs with different devices (MiSeq, 454 or HiSeq system). WGS data revealed the integration of a chloramphenicol resistance gene, the deletion of the endogenous riboflavin (rib) operon and presence of four putative plasmids harbouring rib operons. Event- and construct-specific real-time PCR methods for detection of the GM strain and its putative plasmids in food and feed products have been developed.

Introduction

Riboflavin (7,8-dimethyl-10-[(2S,3S,4R)-2,3,4,5-tetrahydroxypentyl]benzo[g]pteridine-2,4-dione, vitamin B2) is a water-soluble vitamin naturally synthesized by many microorganisms and plants. Since not being produced by higher animals, it is an essential micronutrient in animal and human diets.

The product riboflavin is often used in food as an additive but also finds applications in small amounts as the colouring agent E101 or as a nutritional additive in animal feedstuffs (Abbas & Sibirny, 2011). In the past, riboflavin was mainly chemically synthesised for the production of very pure material. Biotechnological developments have resulted in microbiological processes that can compete with the chemical synthesis and nowadays commercial production of vitamin B2 is mostly done by fermentation.

In most cases, microbial synthesis of riboflavin involves genetically engineered selected strains of Escherichia (E.) coli, Bacillus (B.) subtilis, Ashbya (A.) gossypii, and Candida (C.) famata (Abbas & Sibirny, 2011). Among this plethora of genetically modified micro-organisms (GMMs), GM strains of E. coli and B. subtilis are probably the best known (see (Burgess, Smid, & van Sinderen, 2009)for a review), as in these bacteria the riboflavin biosynthetic pathway has been studied extensively (reviewed in (Bacher et al., 2001)). B. subtilis is known as an aerobic endospore-forming bacterium commonly found in nature and generally not considered to have a pathogenic or toxigenic potential. There is a history of safe use in large-scale fermentation production of speciality chemicals of enzymes used in food production processes, and of several traditional ways of food preparation.

Usually, non-sporulating derivatives of the B. subtilis strain 168, which often carry natural mutations inducing riboflavin overproduction, were genetically modified (GM). Introduction of different plasmids harbouring (i) both a (recombinant) B. subtilis riboflavin biosynthetic operon (rib operon, also known as ribDEAHT operon, i.e. including the ribD, ribE, ribA, ribH, ribT genes) under the control of a strong promoter and (ii) antibiotic resistance genes as selection markers (e.g. cat, tet, ermAM), resulted in GM B. subtilis strains with multiple copies of the rib operon. These strains are able to amplify the riboflavin expression by a magnitude of 10- to 25-fold (Mander and Liu, 2010, Perkins et al., 1999, Smolke, 2009).

According to EFSA guidelines for additives produced with GMM, it is necessary to show that, in the final product, neither the production strain nor its recombinant DNA can be detected (EFSA, 2011). In September 2014, it was notified in the European Rapid Alert System for Food and Feed (RASFF) that a German official enforcement laboratory in Hesse detected viable GM B. subtilis spores in a consignment of vitamin B2 feed additive (80% feed grade) imported from China (RASFF, 2014).

In April 2015, a report of a Belgian official control laboratory was published about the genome sequence of a GM B. subtilis strain (Barbau-Piednoir, De Keersmaecker, Wuyts, et al., 2015). This strain (isolate 2014-3557) was identified by a French competent authority in a lot of vitamin B2 (riboflavin, 80% feed grade) imported to France from China. On basis of next generation sequencing (NGS) data of a Belgian institution a TaqMan® qPCR method (named VitB2-UGM) for specific detection of this EU-unauthorized GM riboflavin-overproducing B. subtilis was developed (Barbau-Piednoir, De Keersmaecker, Delvoye, et al., 2015). However, the TaqMan® qPCR method targets a junction between the riboflavin biosynthesis genes and the vector backbone. It remains unsolved whether the targeted sequence is integrated into the bacterial genome or present on a plasmid. For the latter case, the detection might fail if the plasmid is lost and the corresponding target sequence is therefore missing. In addition, the authors have not reported the comprehensive molecular characterization of the GM strains’ genome and plasmids.

In the current study, microbiological and molecular analyses of the GM B. subtilis strain found in Germany in 2014 are presented. Whole genome sequencing (WGS) was performed with DNA extracted from two independent isolates to characterize in detail the genome of these riboflavin-overproducing GM B. subtilis strains, and to reconstruct the putative plasmids present. Subsequently, construct- and event-specific PCR-based methods for its detection in food and feed were developed and applied.

Materials and methods

B. subtilis living cells were isolated independently by the Hessian State laboratory (LHL) and the Bavarian Health and Food Safety Authority (LGL) in Germany from a product lot of vitamin B2 feed additive (80%) powder imported from China and analysed in the framework of the RASFF notification reference number 2014.1249 (RASFF, 2014). Microbiological and molecular methods for the cultivation and identification of the microorganism and procedures for DNA extraction are described as Supplementary Material.

PCR analyses

Real-time PCR methods were applied to screen for the presence of DNA sequences from recombinant pUC plasmids (Table 1).

Table 1

Primers and probes used in this study.

Target	Name	Oligonucleotide sequence (5′-3′)	Reference
pUC-cloning vectors	401-F_PUC 18-F	TgT CgT gCC AgC TgC ATT A	(Mäde, Reiting, Strauch, Ketteritzsch, & Wicke, 2008)
	401-R_pUC 18-R	gAg CgA ggA AgC ggA AgA g
	401-Tex_Tm-Puc18	TexasRed–AAT Cgg CCA ACg CgC gg-BHQ2
	IPC-fw	TgT gAA ATA CCg CAC AgA Tg	(Messelhäusser et al., 2007)
	IPC-re	AgC Tgg CgT AAT AgC gAA G
	IPC-S	HEX-gAg AAA ATA CCg CAT CAg gC-TAMRA
Chloramphenicol acetyl transferase gene (cat)	356-F_Cat-Staph-F	CAg CTT TTA gAA CTg gTT ACA ATA gCg	This work
	356-R_Cat-Staph-R	gCA TgA TAA CCA TCA CAA ACA gAA T
Plasmid pAM-beta1 (ermAM adenine methylase gene -)	357-F2	CgT CTA TTg AAT TAg ACA gTC ATC TAT TCA	This work
	357-R2	Tgg AAC ATC TgT ggT ATg gCg
Integration site of cat in B. subtilis isolate e871	558-F	CgA gCT TTT gCg CgT ATA	This work
	558-R	gCC ATT CCA ATA CAA AAC CAC ATA
	558-Tm	FAM-Cgg ATC TAA CgC ATg CTC CgC A-BBQ
Plasmid pGMBsub01 (junction of pUB110 to pUC19)	690-F	gAT gAA TTA TAT CAA CAT ATT AAg CCT TTg g	This work
	690-R	gCT Atg ACC ATg ATT Acg CCA Ag
	690-Yak	Yak-AAg ATC Cgg ggA ATT gCT gCA gg-BBQ
Plasmid pGMBsub01 (junction of B. amyloliquefaciens rib-operon to pUC19)	691-F	CgA TTA AgT Tgg gTA ACg CCA	This work
	691-R	TTC TCT AAA gAA AAC TgC TCg TAC g
	691-Tm	FAM-ACg gCC AgT gAA TTC gCA AgA Cg-BBQ
Plasmid pGMBsub02 (junction of deleted B. amyloliquefaciens rib-operon)	804-F1	AgA CCg CgT TTA Cag TCA gCA T	This work
	804-R1	CTCg AAT TCT TTT TTC gTT CCA A
	804-Tm	FAM-ACC ACA AgC TgA CCg AAT ATg Cgg AT-BBQ
Plasmid pGMBsub03 (junction of 3′-rib-fragment to pUC19)	693-F	TCgTgCACAgCTTgAAATCTAgA	This work
	693-R	ggA AAC AgC TAT gAC CAT gAT TAC g
	Cy5-693	Cy5-CCT CTA gAg TCg ACC TgC Agg CAT gC-BBQ
Plasmid pGMBsub04 (junction of rib-operon-fragment to pSM19035)	694-F	CAT TCg ATT gTg CgA gCg	This work
	694-R	Tgg TAT TTT TTg TAT TCA gCg TAA Cag ACA TAA T
	694-Tm	FAM-Cag gCg AAT TCC AgT TAA ATT CCg TgT Agg--BBQ

Primers for screening and detection of an erythromycin resistance gene (ermAM) and of the chloramphenicol acetyl transferase gene (cat) were designed using the Primer Express 3.0 software (Life Technologies Inc.) on the basis of the Streptococcus faecalis plasmid pAM-beta1 adenine methylase gene (GenBank:Y00116) and the sequence information for plasmid pC194 of Staphylococcus (S.) aureus (GenBank:K01998.1), respectively (Table 1).

WGS data were used to develop a GMM event-specific real-time PCR assay targeting the integration site of the chloramphenicol resistance gene (cat) in the genome of the B. subtilis isolate. Further real-time PCR assays were designed to detect the putative recombinant extra-chromosomal plasmids (see Supplementary Material).

Conventional PCR was done in 25 µL using a 10× PCR buffer (Qiagen Inc.) with 15 mM MgCl2, 0.5 µM of each primer, 0.625 U Taq polymerase (HotStar, Qiagen Inc.) and 5 µL of template DNA corresponding up to 500 ng DNA. For thermal cycling, an initial denaturation step of 15 min at 95 °C was followed by 45 cycles of 30 s at 95 °C, 45 s at 60 °C and 45 s at 72 °C with a final elongation step of 7 min at 72 °C.

Real-time PCR was performed in an ABI PRISM 7500 (Applied Biosystems) and 25 µL PCR buffer (QuantiTect Multiplex PCR Mix, Qiagen Inc.) with 0.4 µM of each primer, 0.1 µM probe and 5 µL of template DNA corresponding up to 500 ng DNA. Thermal cycling conditions were a denaturation step of 15 min at 95 °C followed by 45 cycles of 1 min at 94 °C and 1 min at 60 °C. Template DNAs tested by the different PCR methods were either extracted from the different isolates of LHL and LGL or from isolate 2014-3557 of the French competent authority (kind gift of the Scientific Institute of Public Health, Brussels – WIV-ISP).

Whole genome sequencing (WGS)

Three WGS experiments using NGS were performed starting from B. subtilis DNA isolated by the two German laboratories (LHL and LGL) and by the Joint Research Centre (JRC) of the European Commission (Italy).

The analysis of the DNA sample isolated by the LHL were performed by a NGS service provider (StarSeq Inc., Germany) using a MiSeq apparatus (Illumina Inc.). For library preparation, 1 ng of extracted DNA was used for application in the Nextera XT DNA library preparation kit (Illumina Inc.). The generated genomic library was sequenced using the MiSeq Reagent Nano Kit v2 300 cycles (Illumina Inc.) and the pair-end option of 2 × 150 bp of the MiSeq sequencing system. Sequencing was monitored using the ‘Sequencing Analysis Viewer’ program (Illumina Inc.).

The NGS analysis of the DNA sample isolated by LGL was also carried out using Illumina Nextera XT library preparation of 1 ng purified DNA. After quality control, the library was sequenced using an Illumina HiSeq 1500 device using the pair-end flowcell v4 and the HiSeq SBS kit v4, 2 × 50 bp chemistry.

The third NGS analysis was performed by the JRC using another DNA sample from LHL. Here, NGS was done using a GS Junior System (GS Junior System, 454 Life Sciences, Roche Applied Sciences). Rapid libraries (medium length 400–600 bp) were prepared using Rapid library preparation kit (Roche) and all the steps were conducted in accordance with the manufacturers’ instructions.

Bioinformatics

Quality of NGS reads

After adapter trimming, the quality of the NGS Illumina reads was analysed with the Fast QC program (version 1.2.10, default setting) (Andrews, 2010).

For the NGS Roche 454 reads, the quality was inspected by using Roche software (gsRunBrowser, version 2.9).

Assembly and mapping of NGS reads

The following assemblers and mapper tools were used to analyse the produced NGS reads:

Burrows-Wheeler Aligners program (MiSeq Reporter adapted version, default setting) for mapping of Illumina reads. The genome of B. subtilis subsp. subtilis, strain AG1839 (NCBI No. CP008698) was used as a reference;

Spades (version 3.6.2, default setting) (Bankevich et al., 2012) and Velvet (version 1.2.10, K-mer length 49, low coverage 0, high coverage 200, maximum contig length 200, expected genome coverage 1) (Zerbino & Birney, 2008) for de novo assembly of Illumina reads;

runAssembly (Newbler, version 2.9, default setting) for de novo assembly of 454 reads;

runMapping (Newbler, version 2.9, default setting) for mapping of 454 reads on B. subtilis subsp. subtilis strains and for detection of SNPs, indels and inversions.

All cited tools were run with parameters specific for bacterial genomes, as suggested by their developers.

At LGL, data analysis was carried out on a local Galaxy instance and genome assembly was done using Velvet combined with the VelvetOptimiser tool (version 1.0.0; threads = 1, start hash length = 23, end hash length = 45, Kmer optimisation metric = N50, coverage optimisation metric = >1 kb, read type = short paired reads) (Zerbino, 2010). Manual inspections were carried out using the Artemis version 11 (Carver, Harris, Berriman, Parkhill, & McQuillan, 2012) and Mauve version 2.4.0 (Edwards & Holt, 2013) programs with the default settings.

At the JRC, similarity searches and comparisons with other sequences (refseq, nt/nr and pat NCBI databases), were performed by running a local installation of the BLAST (version 2.3.0+) suite (Altschul, Gish, Miller, Myers, & Lipman, 1990). The summary of the similarity searches and the names of the reference strains used are listed in Table 2.

Table 2

Results of similarity searches against known reference strains and their genome sequences.

Species (subsp.)	Strain (substrain)	Number of large deletions	Number of large insertions	Number of large inversions	Number of SNPs	GenBank accession No.
B. subtilis (subtilis)	168	4	1	0	541	NC_000964
B. subtilis (subtilis)	6051-HGW	4	1	0	638	CP003329
B. subtilis (subtilis)	AG1839	5	2	1	887	CP008698
B. subtilis (subtilis)	OH 131.1	≥5	≥2	≥1	≥1000	CP007409
B. subtilis (subtilis)	JH642 (AG174)					CP007800
B. subtilis (subtilis)	BSP1					CP003695
B. subtilis (subtilis)	RO-NN-1					CP002906
B. subtilis (subtilis)	BAB-1					CP004405
B. subtilis	PY79					CP006881
B. subtilis	BEST7003					AP012496
B. subtilis	BEST7613					AP012495
B. subtilis	QB928					CP003783

16S and 23S rDNA analyses

The contigs obtained from assembly were identified by using a local and installed copy of RNAmmer (version 1.2, default setting).

Results

Microbe identification

Immediately after the RASFF notification for a lot of vitamin B2 feed additive (80% feed grade) imported from China, the LHL laboratory and subsequently the LGL laboratory were able to independently isolate living spore-forming bacteria from this batch of vitamin B2 product. Cultures of isolated cell colonies were producing an intense yellow colour secreted to agar plates or broth suggesting the synthesis of riboflavin. Further microbiological and molecular characterization led to the assignment of the isolated bacteria to the species B. subtilis (see Supplementary Material):

Microbiological analyses of five single bacterial colonies showed mass spectrometry (MALDI-TOF) score values above 2.0 and a 99% confidence value in a biochemically based identification system.

Sequence analysis of amplified parts of the 16S–rRNA-gene using DNA extracted from the vitamin B2 product as a template showed a 100% agreement with the corresponding B. subtilis sequence.

These results immediately highlighted that the product was not compliant with the European food and feed safety regulatory requirements for feed additives containing GMMs (EU, 2003a, EU, 2003b).

Identification of the genetically engineered elements

Screening PCR analyses of the vitamin B2 product

The first real-time PCR screening tests (with DNA extracted from the B. subtilis cultures isolated from the vitamin B2 product) targeted DNA sequences of recombinant pUC plasmids (Messelhäusser et al., 2007, Mäde et al., 2008). The two PCR methods target the left and right DNA junctions of the lac-operon from E. coli to sequences of the plasmid pBR322, a DNA construct which is present in all pUC plasmids and derivatives thereof. Clear amplification profiles with Cq-values in the range of 12–14 were obtained indicating the presence of these cloning vector sequences (data not shown).

The review of scientific literature on GMM strains in the context of riboflavin production revealed that antibiotic resistance genes are used as selection markers and thus potentially present in production strains. The most likely candidates were the erythromycin resistance gene (ermAM) and the chloramphenicol acetyl transferase gene (cat). The presence of these antibiotic resistance genes was tested on the B. subtilis strain isolated by LHL using the developed PCR systems (Table 1). The PCR tests resulted in PCR products of 385 bp and of 396 bp corresponding in size to parts of ermAM gene and cat gene, respectively.

Next generation sequencing (NGS)

At the time the RASFF notification was published (end of 2014), three independent NGS experiments were performed. The DNA of the LHL isolate was analysed in two independent NGS runs (using a MiSeq or a GS Junior System), while DNA of the LGL isolate was sequenced in a single run (using a HiSeq system). The obtained sequence reads were independently assembled by the institutions. The aim was to test the ability of different NGS devices and different throughputs to identify the B. subtilis strain and to detect artificially introduced genetic elements. The outcome of these experiments is summarised in Table 3.

Table 3

Summary of the NGS experimental analyses results.

Sample	Institution/NGS platform (specifications)	Read Length	Obtained throughput	BC within Q30	BC error rate	Average assembled contigs depth
LHL	JRC/Roche 454 (WGS modus)	400 bp	83 Mb	99.7%	0.10%	18×
LHL	StarSeq Inc./Illumina MiSeq (paired-end modus)	2 × 150 bp	250 Mb	94.3%	0.18%	58×
LGL	LGL/Illumina HiSeq 1500 (paired-end modus)	2 × 50 bp	12 Gb	94.2%	0.18%	500×

To identify the closest related B. subtilis strain, contigs corresponding to the 16S and 23S rRNA genes were used as a query in BLAST analyses (Altschul et al., 1990) on all known complete B. subtilis genomic sequences available in the GenBank database at the time of the analyses (end of 2014). The LGL, LHL and JRC complete 16S rDNA sequence assemblies are 100% identical to the sequence of the B. subtilis strain 168 and also to 16S rDNA sequence of the French isolate 2014-3557 (Supplementary Material, Fig. S2). Genomes with the best alignment scores (minimum 99% identity over the full length of the locus) were selected and each one was used as a reference genome to map the reads produced (Table 2). It is expected that the better the reads map to a reference genome, the less differences (i.e. indels, inversions and SNPs) are present between the sequenced genome and the reference genome. B. subtilis subsp. subtilis strains 168, 6051-HGW and AG1839 showed to have genomes with the highest similarity to the B. subtilis strain present in the vitamin B2 lot imported from China. These genomes are also fully covered by the reads, with the exception of the detected indels (Table 2). The genomic sequences of B. subtilis strains AG1839, 6051-HGW and 168 show high degree of identity (Kabisch et al., 2013, Smith et al., 2014) and the sequence reads generated here presented some of the strain-specific mutations, but also other variations not reported up to now (data not shown). Therefore, the sequenced isolates were considered as a B. subtilis subsp. subtilis closely related to strain 168, but with specific differences. It is noted, that strain 168 has a long history of laboratory use because of its competence for DNA uptake and transformation (Barbe et al., 2009).

In 2015, Barbau-Piednoir, De Keersmaecker, Wuyts, et al. (2015) announced the sequencing of a GM B. subtilis isolate overproducing riboflavin. These authors did not provide a complete genome sequence but a total of 39 gap-closed scaffolds consisting of 143 contigs with a maximum gap-closed scaffold size of 1,018,461 bp and a minimum size of 370 bp. These contigs were compared to the sequences generated in this study. The comparison (made by BLAST) revealed an almost identical genome with more than 99% sequence identity over the whole length of all the 143 contigs.

In order to detect artificially introduced elements, chromosomal rearrangements were investigated by using the runMapping 454 software and the genome sequence of strain 168 as reference. When compared to this reference, the two sequenced isolates shared the same variations. In particular, four deletions and one insertion were identified:

Nucleotides 529,421 to 549,933. This deletion involves ICEBs1, an integrative and conjugative element (ICE) integrated in the trnS-leu2 gene in B. subtilis, which has been described earlier (Auchtung, Lee, Monson, Lehman, & Grossman, 2005).

Nucleotides 2,428,481 to 2,436,938. This deletion includes ribD, ribE and ribAB genes, in combination with the fswa flavin riboswitch, which is known to be involved in the fine regulation of the rib operon (Winkler, Cohen-Chalamish, & Breaker, 2002). As no other ribDEAHT operon was found on the main chromosome, this indicates that the strain is unable to produce riboflavin without additional plasmid-encoded rib operons. Reads spanning the deletion site are present, while reads that correspond to ribD, ribE and ribAB genes are not included by the software in this region, but are rearranged in a different contig, marked as not present in the reference.

Nucleotides 3,236,246 to 3,236,442. This deletion involves the 3′ end of the hypothetical yufK gene encoding a putative transmembrane protein.

Nucleotides 3,812,644 to 3,812,747. This deletion occurs in the rpoE gene. This gene is known to encode the DNA-directed RNA polymerase subunit delta involved in both the initiation and recycling phases of transcription; in the presence of the delta subunit, RNA polymerase displays an increased specificity of transcription, a decreased affinity for nucleic acids and an increased efficiency of RNA synthesis because of enhanced recycling (Hiratsu, Amemura, Nashimoto, Shinagawa, & Makino, 1995).

A unique identified insertion occurs on chromosomal nucleotide position 1,764,920, and it disrupts a gene (recA, alternative name recE) encoding a multifunctional protein involved in homologous recombination and DNA repair. This insertion is a fragment (1,290 bp) that contains the complete cat gene (European Nucleotide Archive (ENA) accession number LT622644). According to the sequence assembly, it has been cloned within a ClaI restriction site (AT|CGAT). The cat gene encodes a protein that confers chloramphenicol resistance. This finding supports the idea that the chloramphenicol resistance gene present in the chromosome of the analysed B. subtilis strain was intentionally introduced by a crossing over recombination event (Fig. 1).

Fig. 1

Graphic representation of the B. subtilis genomic region in which a transgenic cassette carrying the cat gene has been inserted within ClaI site of the genomic recA gene, thus disrupted. The location event specific detection method is also represented on the junction between the 3′ prime transgenic cassette and the genome.

Finally, when compared to the B. subtilis strain 168, more than 400 potential SNPs have been identified. Searches for specific mutations conferring resistance to antibiotics or peculiar bacterial features, like auxotrophies for amino acids or related compounds were performed. Some examples are described in the Supplementary Material.

GMM event-specific detection method

The integration site of the cat gene in the genome was identified by bioinformatics analyses on the assembled NGS data. The determined DNA sequence across the junction between cat and the recA gene was used to develop a real-time PCR-based event-specific method (Fig. 1).

Details for this PCR system (No. 558) and the specificity are given in the Supplementary Material (Table S1). The event-specific method showed equal detection capability for DNA extracts from the LGL and LHL isolates and for DNA of the French isolate 2014-3557 analysed by NGS (Barbau-Piednoir, De Keersmaecker, Delvoye, et al., 2015). The results of the PCR tests are given in the Supplementary Material (Table S2).

Recombinant extra-chromosomal plasmids

Apart from the described insertion, no other chromosomal integration site could be identified. On the contrary, a fraction of the obtained sequence reads did not map on the reference genome. They could be assembled into four distinct contigs that do not show any read joining them with the genomic sequence of B. subtilis strain 168. Moreover, reads overlapping both their ends were identified and thus marked by the assembler tools (such as newbler) as potentially circular molecules, suggesting a plasmid-like structure. Similarity searches with publicly available nucleic acid databases revealed correspondences with artificial plasmids. For these reasons, we assume that these contigs correspond to four different recombinant extra-chromosomal plasmids (Table 4):

Table 4

Description of the four extra-chromosomal plasmids identified in this study.

Plasmid name	Size	Genes conferring resistance to	RibDEAHT operon	RibDEAHT origin	Similarity to known vectors	Found in sample	Presence described in Barbau-Piednoir et al. (2015)
pGMBSub01	11,378	ampicillin, kanamycin, bleomycin	Full	B. amyloliquefaciens	pUC19, pUB110	LGL, JRC, LHL	Yes(splitted in 5 contigs
pGMBSub02	7,580	ampicillin, kanamycin, bleomycin	Truncated	B. amyloliquefaciens	pUC19, pUB110	LGL, JRC, LHL	Yes(splitted in 5 contigs)
pGMBSub03	8,544	ampicillin, tetracycline	Truncated	B. subtilis	pUC19, pLS1	LGL	Yes(splitted in 2 contigs)
pGMBSub04	29,760	erythromycin	Full	B. subtilis	pSM19035	LGL	Yes(splitted in 5 contigs)

pGMBsub01 (Fig. 2, panel A) (submitted to the European nucleotide Archive – ENA - accession number LT622641). This putative plasmid contains the full B. amyloliquefaciens ribDEAHT operon (ribD, ribE, ribA, ribH, ribT genes, 11,378 bp), together with a B. amyloliquefaciens conserved hypothetical protein (upstream of the rib operon in the B. amyloliquefaciens) and (downstream) the B. amyloliquefaciens genes encoding the segregation and condensation proteins A and B (ScpA and ScpB). The scpA open reading frame (ORF) is present as full length, the scpB ORF is truncated and directly linked (not in-frame) to the S. aureus replication initiation protein B (RepB). The plasmid also carries an ampicillin resistance gene from the pUC19 vector and coding sequences for the kanamycin and bleomycin resistance genes from the vector pUB110 (GenBank M19465.1).

Fig. 2

Schematic representation of the putative recombinant extra-chromosomal plasmids. The different PCR detection methods developed and applied for experimental tests are indicated by the black bars.

pGMBsub02 (Fig. 2, panel A) (ENA accession number LT622641). This putative plasmid (7,581 bp) is a truncated version of pGMBsub01, as mapped reads overlapping the deleted part are present in the NGS data (see Supplementary Material, Fig. S2). In fact, only the ribD and ribE genes from the B. amyloliquefaciens ribDEAHT operon are present. Like pGMBsub01, it carries an ampicillin resistance gene from the pUC19 vector and coding sequences for the kanamycin and bleomycin resistance genes from the pUB110 vector.

pGMBsub03 (Fig. 2, panel B) (ENA accession number LT622642). Like pGMBsub01 and pGMBsub02, this putative plasmid (8,544 bp) carries an ampicillin resistance gene from the pUC19 vector (GenBank M77789.2). In addition, a full tetracycline resistance gene (tet) cassette is present, together with a (truncated) scpB gene upstream, both probably derived from Streptococcus agalactiae plasmid pLS1 (Lacks, Lopez, Greenberg and Espinosa, 1986; GenBank M29725.1). Differently from the previous plasmids, it includes part of the B. subtilis ribDEAHT operon (only ribA, ribH, ribT genes), together with (downstream) the B. subtilis ORFX6 gene.

pGMBsub04 (Fig. 2, panel C) (ENA accession number LT622643). This putative plasmid (29,760 bp) includes the full B. subtilis ribDEAHT operon (about 12 Kb) together with genes typical for Enterococci and Streptococci plasmids and conferring resistance to erythromycin. At the sequence level, the non-rib operon sequence is identical to that of plasmid pSM19035 (GenBank AY357120), a low-copy-number theta-replicating plasmid of the pathogenic bacterium Streptococcus pyogenes, stably maintained in a broad range of gram-positive bacteria (Lioy, Pratto, de la Hoz, Ayora, & Alonso, 2010). The B. subtilis ribDEAHT operon is surrounded by two EcoRI restriction sites. pGMBsub04 shares with pSM19035 large inverted repeats which are described in the latter as responsible for the erythromycin resistance (Soberon, Lioy, Pratto, Volante, & Alonso, 2011). Furthermore, the presence of the pSM19035 sequences explains the detection of the erythromycin resistance gene (ermAM) at the early stage of the characterization of the isolate.

Sequence reads for plasmids pGMBsub01 and pGMBsub02 were present in data obtained from the LHL and the LGL isolate, whereas for pGMBsub03 and pGMBsub04 reads could be found only in the LGL data (Table 4). Presumably these two plasmids were coincidently lost during the enrichments and cultivations at LHL, most probably because the respective antibiotics were not supplemented to the medium. It is noted, that NGS sequence data for all four putative plasmids identified for the LGL isolate are also present in the data of the French isolate 2014-3557.It is thus plausible that there is no incompatibility among the plasmids. In the course of the analyses of the B. subtilis isolate, five real-time PCR methods for detection of the extra-chromosomal recombinant plasmids were developed. Details on the specificity tests are included in the Supplementary Material (Table S2). In conformity with the PCR results obtained for the event-specific method, all plasmid specific methods reacted positive for the DNA extracted from isolate 2014-3557 (Barbau-Piednoir, De Keersmaecker, Delvoye, et al., 2015).

Discussion

In this study, NGS was performed for the molecular characterization of a GM B. subtilis strain overproducing vitamin B2. The aim was to use NGS coupled with bioinformatics to identify the nucleotide sequence of all genetic modifications inserted into the microorganism's genome and to characterize complementing extra-chromosomal recombinant plasmids. WGS by NGS may represent a promising approach for gaining sequence-based information for detailed molecular characterization of insertions in GM organisms (Pauwels et al., 2015, Wahler et al., 2013, Yang et al., 2013). However, as demonstrated by these studies, the approach requires appropriate NGS devices, specific expertise in bioinformatics and the respective analysis tools. Difficulties in the distinction between true variants and sequence alterations and the challenges in making error corrections (for repetitive regions, uncalled bases, ploidy etc.) are currently discussed (Akogwu, Wang, Zhang, & Gong, 2016). On the other hand, NGS may be applied to rapidly generate sequence data in order to identify unknown genomic insertions/deletions and genetic modifications with potential application as basis for developing specific detection methods (Arulandhu et al., 2016, Barbau-Piednoir, De Keersmaecker, Delvoye, et al., 2015). The present work is based on the use of two second generation sequencing technologies (Illumina and 454). Both sequencing chemistries deliver unparalleled accuracy, with a vast majority of bases scoring at Q30 and above (see Supplementary Material). This level of accuracy is ideal for a range of sequencing applications, including clinical research. By mapping the generated reads to the B. subtilis reference genomes, it was possible to identify in the recA gene an insertion site in the genome of the isolated strain where recombinant DNA is integrated. Presumably the recA gene was disrupted to avoid homologous recombination of the production strain. The complete nucleotide sequence of this unique insertion could be determined and identified as the complete chloramphenicol resistance gene (cat) sequence, which may serve as selection marker. The junction of this artificial and stably integrated cassette as target for the design of an event-specific detection method. The specificity was verified by the experiments with DNA samples derived from the LHL, LGL and the French isolate 2014-3557 (see Supplementary Material).

The genome of the isolates carries a deletion of the endogenous ribDEAHT operon, indicating that the strain is unable to produce riboflavin without the recombinant plasmids encoding the rib operon. The absence of the endogenous ribDEAHT operon is presumably complemented with four recombinant plasmids carrying full or truncated ribDEAHT operon from B. subtilis and B. amyloliquefaciens. Together, these operons carry one or more specific antibiotic resistance genes for selection purposes and stable riboflavin expression during fermentation. In the DNA samples sequenced by JRC and LHL only two of the four putative plasmids could be identified (pGMBsub01 and pGMBsub02). Probably during non-selective cultivation of bacteria without antibiotic supplements, the plasmids pGMBsub03 and pGMBsub04 have been lost. The construction of these recombinant plasmids has not been reported by (Barbau-Piednoir, De Keersmaecker, Delvoye, et al., 2015). This study provides the corresponding sequences as unrelated contigs not distinguished from the chromosomal contigs. Presumably, the local similarities, like the common pUC19 backbone regions, do not allow untangling of ambiguities during the assembling procedure. Only the pGMBsub04 plasmid shows a deletion of 200 bp if compared to the counterpart contigs published previously (Barbau-Piednoir, De Keersmaecker, Wuyts, et al., 2015). This deletion puts the neighbouring additional copies of tau and gamma genes in-frame. The presence of the reads across the point of deletion suggests that pGMBsub04 is presumably present in two forms, i.e. a long one found and reported previously (Barbau-Piednoir, De Keersmaecker, Wuyts, et al., 2015) and a short one sequenced within our study of the isolate, maybe the result of a deletion event naturally occurring.

We further investigated the origin of the GM B. subtilis strain. This strain seems not to be described in the literature, apart from the recent study (Barbau-Piednoir, De Keersmaecker, Wuyts, et al., 2015). However, a detailed search in patent databases revealed the presence of an application patent from Russia (Mironov et al., 2004), in which the authors describe a method for producing riboflavin. In this document, the applicants developed and claimed both GM B. subtilis strains and plasmids to over-produce riboflavin. No sequence information is provided, but the following detailed descriptions are reported:

A vector named pEK14 is described which is identical to pGMBsub01, i.e. consisting of the E. coli plasmid pUC19, a kanamycin nucleotidyl transferase gene (kan) from the pUB110 plasmid, and a desensitized (ribO2 mutation) riboflavin operon from B. amyloliquefaciens.

The development of a GM B. subtilis GM51 strain, that is made chloramphenicol resistant by inserting the chloramphenicol acetyl transferase gene (cat) from plasmid pBT69 into the Clal site of the recE gene.

The use of a B. subtilis GM51 strain, which is auxotrophic for riboflavin.

The latter two descriptions are in complete agreement with the identified insertion and rib operon deletion at the chromosomal level.

Furthermore, the authors describe the transformation of B. subtilis GM51 strain with plasmid pMX45 containing the full rib operon from B. subtilis. pMX45 is referenced as a patented recombinant plasmid described in an older French patent (FR Pat. 2546907) (Stepanov, Kukanova, Glazunov, & Zhdanov, 1977). Here the pMX45 plasmid is described as a 30 Kb DNA molecule carrying the full B. subtilis rib operon inserted into the pMX33 plasmid with an erythromycin resistance operon. A bibliographic search revealed that the patent authors published two book chapters (Rabinovich et al., 1985, Rabinovich et al., 1984) in which they describe in detail the development of the recombinant plasmid pMX33 by ligation of an EcoRI-cut B. subtilis DNA fragment containing the full 12 Kb rib operon with an EcoRI-cut vector pMX30. pMX30 in turn is defined as a deletion mutant of plasmid pSM19035 (Rabinovich et al., 1985) with the large inverted repeats conferring erythromycin resistance, like the sequenced pGMBsub04. For these reasons, we argue that pGMBsub04 corresponds to either the pMX45 plasmid or its parent pMX33.

At the time of the RASFF alert, a request for clarification was sent by the German diplomatic service to the Chinese company that marketed the vitamin B2 feed additive (80%). A brief description of a strain claimed to be used as the vitamin B2 production strain was provided, together with the sequence information of a recombinant plasmid that confers resistance to tetracycline. This plasmid sequence was found to be identical to the pGMBsub03 sequence. The company claimed that a plasmid named pMX45 was used for the transformation of the B. subtilis production strain, which confirms our hypothesis that plasmid pGMBsub04 correspond to pMX45. A comparison of the information provided by the Chinese company and our findings is shown in Table S2 of Supplementary Material. In particular, we did not find any pUC19 integration at the genomic level, but additional recombinant plasmids neither mentioned by the Chinese company nor expected to be present. The differences in the information provided by the Chinese company and the results of the molecular characterization strongly imply that the production strain must have been contaminated or switched before or during production.

In response to the detection of the presence of GMOs in rice products exported from China, to a case of Bt63 presence in feed additive to the EU and to the case described in our study which lead to a notification in the RASFF system, the European Commission’s Food and Veterinary Office (FVO) carried out an audit in China in order to evaluate their control systems for GMO. It is reported that the Chinese competent authorities carried out investigations at all vitamin B2 producers, but no irregularities were identified (FVO, 2016).

In order to evaluate the possible microbiological risks for human and animal health, the microbiological status of the production process requires consideration with regards to the safety and identity of the producer organism used in the fermentation process. Riboflavin production by fermentation of specific GMMs must therefore involve the approved safety assessment of the genetic modification and the respective strain. The continuous quality control testing of the production process and cultural purity is required to exclude contaminations and/or impurities. Degradation of all DNA needs to be verified by the producer in polymerase chain reaction (PCR) tests (Hermann & Schurter, 1995).

Food and feed additives, produced in closed production systems with GMMs, do not fall under the scope of Regulation (EU) No. 1829/2003 on GM food and feed (EU, 2003a). According to its recital 16, this Regulation covers food and feed produced “from” a GMO but not food and feed “with” a GMO. This implies, and it is laid down in EFSA guidelines, that for additives produced with GMMs it has to be shown that, in the final product, neither the production strain nor its recombinant DNA can be detected (EFSA, 2011). In these cases, an event-specific GMO identification method is not required. The product is not subject to any labelling requirement according to Regulations (EC) 1829/2003 on GM food and feed or No. 1831/2003 (EU, 2003a, EU, 2003b), assuming that the producing company has verified the purity of the final product. To our knowledge, the finding of a contamination of living GM B. subtilis and of recombinant DNA in a marketed product has induced an increase in official controls of additives and enzymes produced with GMMs, at least in Germany. Such measures are based on the precautionary principle and the risk-based approach of enforcement to ensure the high level of food and feed safety in the EU and its Member States (EU, 2002). Since the current case, no other report concerning a GMM has been notified in the RASFF. It is thus assumed that the finding of a GMM in this product group was an exceptional and singular case. However, the lack of a specific detection method for a respective production strain made the analytical work quite challenging.

Conclusions

We have described the value of the use of NGS approaches in the context of assisting food/feed competent enforcement authorities and official control laboratories in their investigations concerning the presence of GMMs in biotechnological products. The approach described in our study allowed to provide a detailed molecular characterization of an unknown GMM and thereby facilitated the risk assessment and as well as the development of specific PCR methods for its detection.

Nevertheless, the study showed that the bioinformatics analysis of the data produced by NGS is still a challenging task and would require the development of adapted bioinformatics tools in order to be implemented for routine data analysis and management in the frame of GMO characterization and detection. The investigation of the GM B. subtilis strain showed at the same time, that in-depth literature review is fundamental to interpret the sequence data and to fully understand and characterize the sample, especially if completely unknown.

Fundings

This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profit sectors.