The test materials were selected in accordance with the International Cooperation on Harmonisation of Technical Requirements for Registration of Veterinary Medicinal Products guideline 34, which serves as a standard for assessing the safety, quality, and efficacy of veterinary drugs, particularly biological products (Table 1). To assess potential contamination, the guideline mandates the use of five specific organisms belonging to the class Mollicutes. Although these organisms belong to different families, they include one species of Acholeplasma (A. laidlawii) and four species of Mycoplasma (M. fermentans, M. orale, M. hyorhinis, and M. synoviae). These species were chosen based on several factors, including antibiotic susceptibility, culture requirements, potential for contamination, and pathogenicity (10). These organisms originate from a range of hosts, including mammals, birds, and humans. The five species used in this study were obtained from the American Type Culture Collection. In addition, E. rhusiopathiae was selected as a test material for the development of bioinformatics-based analytical methods due to its genetic similarity to Mycoplasma, which makes it difficult to distinguish using PCR. A commercial vaccine approved in the Republic of Korea containing CSFV and E. rhusiopathiae antigens was used in the study.

Initially, both E. rhusiopathiae and Mycoplasma species were diluted using phosphate-buffered saline. E. rhusiopathiae was adjusted to achieve a concentration equivalent to two vaccine doses. Equal volumes (1 mL each) of the E. rhusiopathiae vaccine and serial dilutions of Mycoplasma species were then mixed. This ensured that the final concentration of E. rhusiopathiae in the mixed sample corresponded to a single vaccine dose, while Mycoplasma concentrations varied. The spike assay confirmed that E. rhusiopathiae remained at 1.6 × 10⁹ CFU/mL, while the starting concentrations of the five Mycoplasma species before dilution were as follows: 1.35 × 10⁶ CFU/mL for A. laidlawii, 1.05 × 10⁶ CFU/mL for M. fermentans, 7.15 × 10⁵ CFU/mL for M. orale, 7.95 × 10³ CFU/mL for M. hyorhinis, and 2.5 × 10⁶ CFU/mL for M. synoviae. In addition to the spiked samples, a separate group containing only Mycoplasma, without E. rhusiopathiae, was prepared to assess the detection threshold by PCR. All experiments, including sample preparation and downstream analysis, were performed in triplicate to ensure reliability and reproducibility.

Along with the spiked samples, 31 negative vaccine samples—corresponding to half of the 62 vials assigned for quality control at APQA in 2024—were also analyzed. These samples consisted of commercial combination vaccines containing only CSFV and E. rhusiopathiae, and were presumed to be free of Mycoplasma contamination. The negative field samples were used to verify assay specificity and to serve as negative controls for downstream PCR and NGS-based analyses.

DNA extraction was performed using an automated nucleic acid platform (Maelstrom 4810, TANBead, Taiwan, China) and a magnetic bead-based protocol with the TANBead Nucleic Acid Extraction Kit (TANBead). Following the manufacturer’s instructions, 300 μL of sample and 10 μL of Proteinase K were used as input. The DNA was eluted in 80 μL of elution buffer.

The 16S rRNA molecule, approximately 1,500 base pairs in length, is a widely conserved sequence that contains essential structural regions along with variable regions that enable differentiation between bacterial species. This makes it a commonly used tool in microbiology for bacterial identification and phylogenetic analysis (7). In this study, primers were designed to amplify the 16S rRNA region (11) with the following sequences (Table 2): universal forward primer (5′-GGC GAA TGG GTG AGT AAC ACG-3′) and universal reverse primer (5′-CGG ATA ACG CTT GCG ACC TAT C-3′). These primers follow the National Regulatory Standards for Veterinary Biologicals established by the Republic of Korea (1). PCR was carried out using the Maxime™ PCR premix (i-Taq, iNtRON Biotechnology) in a total reaction volume of 20 μL comprising 1 μL (10 pmol) of forward primer, 1 μL (10 pmol) of reverse primer, 2 μL of DNA, and 16 μL of distilled water. Thermal cycling conditions included an initial denaturation at 94 °C for 5 min; 30 cycles of denaturation (1 min at 94 °C), annealing (1 min at 60 °C), and extension (1 min 30 s at 72 °C); followed by a final extension at 72 °C for 7 min. PCR products were analyzed by electrophoresis on a 1.5% agarose gel, and the presence of a specific band at 464 base pairs (bp) was confirmed using a UV transilluminator.

The PCR products were submitted to a commercial provider (BIONICS, Daejeon, Republic of Korea) for library preparation, quality control, and sequencing using the BITseq next-generation sequencing service. Sequencing was performed with a target output of 30,000 reads, and data were delivered in FASTQ format containing 150 bp paired-end reads. These reads were used as raw data for downstream analysis.

The raw sequencing data were analyzed using two distinct approaches: a reference-mapping pipeline and a metabarcoding workflow, each employing specialized software tailored to its respective analytical purpose (Figure 2). The reference-mapping method was developed to isolate specific strains from heterogeneous samples. BBMap (v39.01) was used to align sequencing reads to predefined reference genomes (12), while SPAdes (v3.13.1) served as the assembler for genome reconstruction (13, 14). Key features of this pipeline included high-accuracy alignment to reference genomes and targeted extraction of Mycoplasma-derived contigs.

Prior to analysis, reference sequences for E. rhusiopathiae and Mycoplasma spp. were compiled. For E. rhusiopathiae, the 16S rRNA partial sequence of strain ATCC 19414 was obtained from the NCBI reference sequence database. Concurrently, 16S rRNA gene reference sequences were obtained from the SILVA ribosomal RNA gene database (released ver. 138.2, SSURef NR99) (15). To construct a targeted database for mapping, sequences were filtered to retain only those assigned to the families Mycoplasmataceae, Acholeplasmataceae, Metamycoplasmataceae, and Mycoplasmoidaceae according to SILVA taxonomy (Supplementary Data Sheet 1). In parallel, a UniVec database containing adapter, linker, and primer sequences commonly introduced during cloning workflows was incorporated (16). Although the library preparation process did not explicitly confirm the presence of such sequences, vector screening was included as a precautionary measure to mitigate potential contamination and safeguard data integrity. Vector sequence removal and quality filtering are critical in reference mapping to prevent non-specific alignment of contaminant reads to the reference genome. This pre-processing step is essential for reducing false-positive alignments in downstream mapping.

The reference-mapping protocol consisted of four sequential stages, beginning with vector removal and quality filtering. The bbduk.sh script was used for both tasks. In the first stage, reads were aligned to the UniVec database and separated into vector-contaminated reads (outm1, outm2) and cleaned reads (outu1, outu2). The key parameters included k = 31 (k-mer length of 31) and hdist = 1 (maximum permitted Hamming distance of 1). In the second stage, adapter trimming and quality filtering were performed using the same script with the following parameters: ktrim = r (trimming from the right), k = 23 (k-mer length of 23), mink = 11 (minimum k-mer length of 11), hdist = 1 (maximum permitted Hamming distance of 1), qtrim = rl (quality trimming from both ends), trimq = 20 (quality score threshold of 20), and minlen = 50 (minimum read length of 50 bases).

The next two steps involved reference-based alignment using the bbmap.sh script. In step three, reads that did not align with the E. rhusiopathiae reference genome were retained using minid = 0.95 (minimum identity of 95%) and maxindel = 3 (maximum allowable insertion/deletion length of 3 bases). In step four, the remaining reads were aligned to the Mycoplasma reference database with stricter parameters: minid = 0.99 (minimum identity of 99%) and maxindel = 3. Genome assembly was then performed using the SPAdes assembler (spades.py) (14). The detection limit of Mycoplasma was established by evaluating colony-forming unit (CFU) counts and identifying contig formation at corresponding dilution levels.

To enable accurate comparisons across samples, sequencing depth was normalized using the Total Count (TC) method. This involved calculating the total coverage for each contig (contig length multiplied by average coverage) and dividing it by the total number of reads in the sample (17). This normalization approach minimizes technical bias introduced by varying sequencing depths and enables consistent interpretation across datasets. The formula used is as follows:

The metabarcoding analysis was conducted using Qiime2 (v2024.05) to examine sequencing data from diverse microbial communities present in the samples (18). The analysis included quality filtering, sequence clustering, and taxonomic classification at the order level to identify and categorize different species. MicrobiomeAnalyst was also used for downstream data visualization and interpretation (19).

Initial sequencing data were processed in Qiime2 (v2024.05) through a series of bioinformatic procedures (15, 18). Raw sequence reads were imported into the pipeline, and quality control steps were applied. This included preprocessing with DADA2 to remove noise, correct sequencing errors, and identify amplicon sequence variants (ASVs), providing a high-resolution representation of microbial diversity. Following denoising, a feature table was constructed to record the occurrence of each ASV across all samples. Taxonomic classification of ASVs was performed by aligning sequences with the SILVA 16S rRNA reference database (15), followed by filtering to remove low-confidence assignments and potential contaminants. The data were then normalized to account for variations in sequencing depth, using methods such as rarefaction to allow comparability between samples (18, 19). Finally, the processed data were exported for statistical analysis and visualized using various tools within MicrobiomeAnalyst (19). We assessed the prevalence of Mycoplasma in each sample based on visual outputs and evaluated detection limits accordingly.

The detection limits for five Mycoplasma species—A. laidlawii, M. fermentans, M. orale, M. hyorhinis, and M. synoviae—were evaluated (Table 3). Amplification for A. laidlawii was consistently observed at 3.13. log CFU/mL across all replicates, with no detection at lower concentrations. M. fermentans showed a consistent detection threshold of 3.02 log CFU/mL in all three replicates. Similarly, M. orale produced amplification at 3.85 log CFU/mL across all replicates. In the case of M. hyorhinis, the detection limit varied slightly, with amplification observed at 2.9, 1.9, and 1.9 log CFU/mL, respectively; however, 1.9 log CFU/mL was the most consistent threshold. M. synoviae demonstrated reproducible amplification at 3.4 log CFU/mL across all replicates. Collectively, these results indicate variation in PCR detection efficiency among Mycoplasma species, with amplification generally decreasing as target concentrations decline.

The reference-mapping method was applied to assess its performance in detecting multiple Mycoplasma spp. in spiked samples. The results present an overview of the read retention at each step of the reference-mapping process and demonstrate how mapped read counts vary according to bacterial concentration.

Sequencing yielded consistently high read counts across all species: 27,633 reads for A. laidlawii, 33,496 for M. fermentans, 33,210 for M. orale, 36,347 for M. hyorhinis, and 33,904 for M. synoviae, with 100% read retention at the initial stage. After vector sequence removal and quality trimming, nearly all reads were retained, indicating minimal data loss and high initial sequence quality.

Following the first step against E. rhusiopathiae, read counts decreased significantly across species: A. laidlawii retained 2.6% (726 reads), M. fermentans 3.4% (1,129 reads), M. orale 4.2% (1,410 reads), M. hyorhinis 13.0% (4,729 reads), and M. synoviae 7.1% (2,419 reads). A second mapping step against the Mycoplasma database further reduced reads, yielding 0.3% (84 reads) for A. laidlawii, 1.3% (453) for M. fermentans, 1.7% (577) for M. orale, 8.3% (3,018) for M. hyorhinis, and 2.5% (863) for M. synoviae (Table 4). These findings demonstrate the effectiveness of the two-step reference-mapping strategy in removing non-specific reads and emphasize the value of sequential mapping for accurate downstream analysis.

The metabarcoding method was employed to evaluate the relative abundance of Mycoplasma spp. across various dilution levels, offering insights into sample composition and detection limits, For comparison, a table presenting both relative coverage and relative abundance is included (Table 5).

Relative coverage, derived from the reference-mapping method, refers to the proportion of reads mapped to a reference genome out of the total reads in a sample. It is calculated as the sum of the average coverage values across all genome regions divided by the total read count. Relative abundance, in contrast, reflects the proportion of each species within the total microbial population, providing information on species diversity and distribution. Although the relative coverage data are presented for reference, the primary focus is on abundance patterns obtained from the metabarcoding method.

For A. laidlawii, the relative abundance at the highest concentration (α × 10⁵ CFU/mL) was 88.0%. As the dilution increased, abundance dropped sharply, reaching 0% at α × 10⁰ and α × 10⁻¹ CFU/mL, where no contigs were detected. M. fermentans showed a similar pattern, with 94.4% abundance at α × 10⁵ CFU/mL, decreasing to 0% at the lowest concentrations. M. orale followed this trend, starting at 82.2% and declining to 0% with increasing dilution. M. synoviae exhibited the highest starting abundance (92.7%) at α × 10⁵ CFU/mL, which also dropped to 0% at lower levels. In contrast, M. hyorhinis, which began testing at α × 10³ CFU/mL due to limited starting material, showed 89.3% abundance at that level. Higher concentration experiments were not conducted for this species, resulting in “Not Tested” values. Despite this, the metabarcoding successfully detected M. hyorhinis at lower concentrations, demonstrating its utility in identifying less abundant species in spiked samples. Non-specific reads from genera such as Bacillus and Staphylococcus were also detected.

The detection limits of Mycoplasma spp. were compared across three methods: conventional PCR, the reference-mapping method, and the metabarcoding method (Table 3). Across all three methods, standard deviation remained below one log unit across triplicate experiments. PCR exhibited the highest detection thresholds, while both NGS-based methods yielded lower detection limits across all species.

Initial reference-mapping method yielded mycoplasma-specific contigs in 5 of the 31 negative vaccine samples. The contigs were 447, 444, 439, 447, and 449 bp in length, respectively, and all exhibited sequencing depth below 10. BLAST analysis revealed that four contigs corresponded to M. synoviae, whereas the 449 bp contig was identified as M. fermentans. Concurrently, the sequenced data from these samples were reanalyzed using our metabarcoding method pipeline (Qiime2). This method did not detect any reads classifiable to Mycoplasma.