The first whole-genome sequence of a prospective novel sponge-associated Streptomyces strain from Indonesia with a long 1.5 Mbp terminal inverted repeat

Bacterial isolation and culture condition

Streptomyces sp. BTA 1-131 was previously isolated from marine sponge M. sarasinorum in December 2015. The sponge was collected in Tanjung Batu Angus, Lembeh Strait, Bitung, North Sulawesi, Indonesia (1.50572° N, 125.24541° E; depth 13 m) as described in the previous study.⁶ The strain BTA 1-131 was previously identified by sequencing the 16S rRNA gene and is closely related to soil-derived S. kunmingensis NBRC 14463.⁷ The nucleotide sequence of the 16S rRNA was deposited in the NCBI GenBank database under accession number MT280129, and the cultivable isolate was deposited in the Indonesian Culture Collection BRIN (InaCC) under accession number A1205.⁷ In this study, the strain BTA 1-131 was revived from a glycerol stock stored at –80 °C. Prior to genomic DNA extraction, the strain BTA 1-131 was grown on an International Streptomyces Project (ISP)-recommended medium, ISP2 agar medium, at 28 °C, for 14 days. The ISP2 medium contains yeast extract (BD Difco, USA; 4 g L⁻¹), malt extract (BD Difco, USA; 10 g L⁻¹), and dextrose (BD Difco, USA; 4 g L⁻¹). The medium was prepared in 1 L of seawater with an addition of bacteriological agar (BD Difco, USA; 20 g L⁻¹).

DNA extraction and sequencing

Genomic DNA was extracted from Streptomyces sp. BTA 1-131 colonies grown on ISP2 medium agar plate using the Quick-DNA^TM HMW MagBead kit according to the manufacturer’s instructions (Zymoresearch, USA). The integrity of the extracted genomic DNA was assessed by gel electrophoresis (0.8% tris borate-EDTA agarose w/v), while its quality and quantity were analysed by a NanoDrop spectrophotometer (Thermo Fischer Scientific, USA) and a Qubit 4 fluorometer (Thermo Fischer Scientific, USA) using the Qubit^TM dsDNA High Sensitivity assay kit (Thermo Fischer Scientific, USA), respectively. Moreover, total gDNA was used for the input of library preparation for both short read-based Illumina sequencing and long read-based ONT sequencing platforms.

To sequence the genome of the strain BTA 1-131 on the Illumina platform, the DNA library was prepared using the Illumina NextSeq500 system with the NextSeq 500/550 kit v2.5 (300 cycles) to generate 2 x 150 bp paired-end sequence reads following the manufacturer’s protocol (Illumina, USA). For the ONT platform, the DNA library was prepared for GridION ONT using the SQK-NBD114.24 ligation sequencing kit according to the manufacturer’s recommendation (Oxford Nanopore Technology, UK). In brief, total gDNA was repaired using an end prep enzyme mix, generating DNA with 5’ phosphorylated, 3’ dA-tailed ends. Barcodes were ligated with an ONT-compatible adapter. The library was quantified with a Qubit Fluorometer before loading to the FLO-MIN114 flow cell. Sequencing control was conducted in MINKNOW v23.04.6, and the output raw reads were set to POD5 for a separate basecalling process.

Data processing and hybrid assembly

Illumina read processing

Illumina-generated FASTQ files were initially processed using Fastp v0.24.0¹² to assess read quality and remove low-quality bases (Data file 1).¹³ Reads were filtered with the parameters -q 20 -l 75 -c, ensuring only high-quality reads were retained (Data file 2).¹⁴ Unpaired reads (those lacking a mate from either R1 or R2) were kept using the flags –unpaired1 sample_unpaired1.fastq.gz --unpaired2 sample_unpaired2.fastq.gz and were included alongside paired reads for subsequent hybrid assembly with the ONT data.

ONT read processing and hybrid assembly

Raw POD5 files were basecalled at super accuracy using Dorado v0.8.2, generating a simplex BAM file. This BAM file was converted to FASTQ format using Samtools v1.21.¹⁵ Read quality distributions were assessed with NanoPlot v1.44.0,¹⁶ (Data file 1).¹³

Hybrid assembly began by assembling ONT reads first using Flye v2.9.5,¹⁷ (Data file 3).¹⁸ Parameters included --nano-hq -g 8m --asm-coverage 50, reflecting an estimated 8 Mbp genome size and a coverage limit set to 50X. Assembly graphs were visualised with Bandage v0.9.0,¹⁹ (Data file 4).²⁰ Moreover, to refine the ONT-based assembly, Medaka v2.0.1 (Medaka consensus) was applied with recommended settings optimised for v14 reagent kit reads. A final hybrid assembly was then generated by polishing the ONT assembly with the high-quality Illumina reads using Pilon v1.24.²¹ Three iterative polishing steps were performed; in the first iteration, the Illumina was aligned to the Medaka-polished ONT assembly using Bowtie v2.5.4,²² and alignments were processed with Samtools v1.21. The resulting consensus contigs become assembly sources for the second iteration and continue until three iterations of polishing are reached. After these three polishing rounds constituted the final hybrid assembly for subsequent analyses.

Assembly statistics and completeness

Assembly metrics were obtained with the quality assessment tool for genome assembly (QUAST) v5.2.0,²³ reporting the number of contigs, N50, and GC content; gene predictions were disabled for this step (Data file 5).²⁴ Genome completeness was evaluated using the benchmarking universal single-copy orthologs (BUSCO) v5.8.2²⁵ in genome mode, employing the Streptomyces_odb12 lineage dataset (Data file 6).²⁶ The identified single-copy orthologs were visualised via the generate_plot.py script included with BUSCO. CheckM v1.2.3²⁷ was used with the lineage_wf workflow to estimate completeness and contamination, summarised in a table of reported statistics (qa) (Data file 7).²⁸ Assemblies with >95% completeness and <5% contamination were considered high-quality.

Bioinformatic analysis

Genome visualisation and annotation

Genome visualisation was performed using GenoVi v02.16,²⁹ which generated circular genome maps depicting sequence features, contig boundaries, and genomic architecture. Functional annotation of clusters of orthologous groups (COGs) categories was conducted using the integrated DeepNOG³⁰ within GenoVi, with gene sequences queried against the COG database to assign functional categories and identify orthologous proteins across genomes (Data file 8).³¹ The resulting visualisations provided both an overview of overall genomic organisation and detailed functional annotation profiles for downstream analysis.

TIR and plasmid identification

To examine the linear inverted repeats on the strain BTA 1-131 chromosome, a basic local alignment search tool (BLAST) was conducted using the parameters described previously by Jørgensen et al. (2024).³² While putative plasmid contigs were identified with PLASme v1.1³³ using the default database and the “balance” mode parameter. All other settings were kept at default values. Any candidate plasmid contigs were verified by examining the assembly graph in Bandage to confirm their circular topology or unique structural features.

Taxonomic inference and phylogenomic analysis

Taxonomic assignment was carried out on the type (strain) genome server (TYGS)³⁴ by submitting the polished genome assembly (Data file 9).³⁵ Additional inferences were made using GTBD-Tk v2.1.0,³⁶ incorporating Streptomyces kunmingensis strains and the type strain S. albus NRRL B-181 (ASM72588v1), as well as some selected Streptomyces strains associated with marine sponges reported by Xu et al. (2019)³⁷ for comparative analysis. Moreover, Deferribacter desulfuricans SSM1 was chosen as an outgroup. The resulting genome BLAST distance phylogeny (GBDP) tree was visualised in iToL v7.³⁸ The TYGS webserver also provided a DNA-DNA hybridization (dDDH) score.

To assess the genomic similarity between the polished genome assembly of the strain BTA 1-131 and one of the reference S. kunmingensis DSM 41681 (the nearest strain based on the genome-scale GBDP tree), average nucleotide identity (ANI) was calculated using FastANI.³⁹ The reference genome sequence of S. kunmingensis DSM 41681 (ASM3561610v1) was obtained from publicly available genomic repositories, and the strain BTA 1-131 genome assembly was prepared using the following command: FastANI -q assembly_genome.fasta -r reference_genome.fasta -o ani_output.txt, where -q specifies the query genome (the polished genome assembly of the strain BTA 1-131), -r denotes the reference genome (S. kunmingensis DSM 41681), and -o defines the output file storing ANI results. The ANI percentage from the analysis result was used to determine the genomic relatedness between the assembly and the reference strain. Results were interpreted based on established ANI thresholds (90–95%) for bacterial species delineation.³⁹

Complete assembled genome description

The whole genome sequence of Streptomyces sp. BTA 1-131 was obtained from a total of 567,347 sequence reads from the ONT long-read assembly, followed by polishing with the Illumina short-read assembly. This hybrid assembly resulted in a 10,234,869 bp genome size with an N50 median of 9,496,435 bp and a GC content of 71.62%, consisting of one large chromosomal contig and three small contigs. Based on the CheckM analysis, the genome was 99.89% complete with a contamination level of 0.21%. In addition, the overall average coverage of the complete assembled genome was noted to be approximately 170.7X.

Moreover, both linear and circular topologies were found in the genome of Streptomyces sp. BTA 1-131. The BLAST analysis revealed that a linear chromosomal genome of the strain BTA 1-131 contains a long-terminal inverted repeat (L-TIR) of up to 1.5 Mbp at the chromosome end (9,496,435–7,999,105; 1,497,730 bp), resulting in a loop-like feature ( Figure 1). Two additional small linear contigs, 386.1 Kbp and 301.8 Kbp in size, along with the smallest circular contig at 50.5 Kbp, were identified as plasmids. In total, there were 8,986 protein-coding genes (CDSs) with the average open reading frame (ORF) length of 995 bp. Within the genome Streptomyces sp. BTA 1-131, 21 rRNA and 92 tRNA operons were also predicted ( Table 1).

Figure 1. Streptomyces sp. BTA 1-131 assembly using Flye with high assembly depth.

(A) Schematic representation of the Streptomyces sp. BTA 1-131 genome detected, adding up to 1.5 Mbp in length to chromosome size. In addition, a circular and two linear plasmids were detected from the strain BTA 1-131 assembly. (B) L-TIR is indicated with a reverse arrow (pink colour) from the end of chromosome (9,496,435–7,999,105; 1,497,324 bp).

Furthermore, the general genome features of Streptomyces sp. BTA 1-131 are summarised in comparison to its closely related neighbour strain S. kunmingensis DSM 41681 and the type strain S. albus NRRL B-1811 in Table 1, while the schematic of circular maps of the genome, including the annotated COG categories, is presented in Figure 2.

Figure 2. A circular map of each contig from the Streptomyces sp. BTA 1-131 assembly.

Contig of L-TIR displayed together as one chromosome. Labelling from outside to the inside: contigs; COGs on the forward strand; CDS, rRNAs, and tRNAs on the forward strand; CDS, rRNAs, and tRNAs on the reverse strand; GC content; GC skew.

Phylogenomic analysis

As aforementioned, the strain BTA 1-131 was previously identified by sequencing the 16S rRNA gene and is closely related to S. kunmingensis NBRC 14463 with the sequence similarity of 98.4%. Indeed, the results of TYGS species identification suggested a potential novel Streptomyces species for the strain BTA 1-131, as indicated by the dDDH value below the threshold for classification as the same species as the nearest strain, S. kunmingensis DSM 41681 (<70%) ( Figure 3). Moreover, the ANI score between the BTA 1-131 genome assembly and S. kunmingensis DSM 41681 was 89.15%.

Figure 3. Phylogenomic tree construct based on whole genome sequence data of Streptomyces sp. BTA 1-131 and other Streptomyces strains.

Limitations

The main limitation in this study is that although the phylogenomic analysis showed a possibility of the strain BTA 1-131 as a new Streptomyces species, a follow-up taxonomical confirmation through polyphasic study is necessary to establish its species status. Additional taxonomical investigations, including chemotaxonomic, morphological, and biochemical characterisation, would be required to definitively confirm the novelty of this strain and validate its classification within the genus Streptomyces.