Cultivation driven transcriptomic changes in the wild-type and mutant strains of Rhodospirillum rubrum

Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 Available online 20 June 2024 2001-0370/© 2024 The Authors. Published by Elsevier B.V. on behalf of Research Network of Computational and Structural Biotechnology. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Cultivation driven transcriptomic changes in the wild-type and mutant strains of Rhodospirillum rubrum Katerina Jureckova a , 1 , Marketa Nykrynova a , 1 , Eva Slaninova b , Hugo Fleuriot-Blitman c , V´ eronique Amstutz c , Kristyna Hermankova a , Matej Bezdicek d , e , Katerina Mrazova b , f , Kamila Hrubanova f , Manfred Zinn c , Stanislav Obruca b , Karel Sedlar a , * a Department of Biomedical Engineering, Faculty of Electrical Engineering and Communication, Brno University of Technology, Brno, Czech Republic b Department of Food Chemistry and Biotechnology, Faculty of Chemistry, Brno University of Technology, Brno, Czech Republic c Institute of Life Technologies, University of Applied Sciences and Arts Western Switzerland Valais-Wallis (HES-SO Valais-Wallis), Sion, Switzerland d Department of Internal Medicine – Haematology and Oncology, University Hospital Brno, Brno, Czech Republic e Department of Internal Medicine – Haematology and Oncology, Faculty of Medicine, Masaryk University, Brno, Czech Republic f Institute of Scientific Instruments of the Czech Academy of Sciences, v.v.i., Brno, Czech Republic ARTICLE INFO Keywords: Rhodospirillum rubrum RNA-Seq Transcriptome Genome Depolymerase knock-out Gene ontology Metabolism Fructose Acetate Polyhydroxyalkanoates ABSTRACT Purple photosynthetic bacteria (PPB) are versatile microorganisms capable of producing various value-added chemicals, e.g., biopolymers and biofuels. They employ diverse metabolic pathways, allowing them to adapt to various growth conditions and even extreme environments. Thus, they are ideal organisms for the Next Generation Industrial Biotechnology concept of reducing the risk of contamination by using naturally robust extremophiles. Unfortunately, the potential of PPB for use in biotechnology is hampered by missing knowledge on regulations of their metabolism. Although Rhodospirillum rubrum represents a model purple bacterium studied for polyhydroxyalkanoate and hydrogen production, light/chemical energy conversion, and nitrogen fixation, little is known regarding the regulation of its metabolism at the transcriptomic level. Using RNA sequencing, we compared gene expression during the cultivation utilizing fructose and acetate as substrates in case of the wildtype strain R. rubrum DSM 467 T and its knock-out mutant strain that is missing two polyhydroxyalkanoate synthases PhaC1 and PhaC2. During this first genome-wide expression study of R. rubrum, we were able to characterize cultivation-driven transcriptomic changes and to annotate non-coding elements as small RNAs. 1. Introduction Bacteria represent a remarkably diverse group of organisms, which is not surprising, as, according to estimations, there are 10 6 to 10 8 separate prokaryotic genospecies on Earth [1]. While the majority of these microorganisms living in a wide range of natural environments seem to be uncultivable with current techniques, some organisms are very versatile and can prosper under various cultivation conditions and grown on various substrates. The ideal example is Rhodospirillum rubrum, a purple, non-sulfur, Gram-negative facultative anaerobe from the class of Alphaproteobacteria, which was observed for the first time by Esmarch in 1887 [2]. It was shown to grow both aerobically and anaerobically. The absence of oxygen triggering the photosynthesis apparatus for synthesis of membrane proteins, bacteriochlorophylls, and the carotenoids turning the culture purple. Its metabolic versatility is further supported by the fact that it prospers under both heterotrophic and autotrophic conditions. Besides utilizing various organic substrates as saccharides or organic ions, R. rubrum can fix and metabolize inorganic compounds such as CO and CO 2 and is therefore considered a prospective strain for valorization of waste gases such as combustion products or syngas [3,4]. On the other hand, it is capable of producing other gaseous products utilizable as fuels, particularly hydrogen [5]. Last but not least, R. rubrum is a potent producer of biopolymers of the class of polyhydroxyalkanoates (PHAs) in the form of intracellular granules from many carbon sources. When supplemented by butyrate, for instance, it can accumulate up to 50 w/w % of dry weight as PHA [6]. Hence, it hosts a large variety of metabolic pathways that can be leveraged to produce sustainable carbon substrates. * Corresponding author. E-mail address: [email protected] (K. Sedlar). 1 Contributed equally Contents lists available at ScienceDirect Computational and Structural Biotechnology Journal journal homepage: www.elsevier.com/locate/csbj https://doi.org/10.1016/j.csbj.2024.06.023 Received 15 March 2024; Received in revised form 11 June 2024; Accepted 18 June 2024 Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 2682 R. rubrum is considered to be a model strain for studying the conversion of light energy to chemical energy [7], hydrogen biosynthesis [8], formation of photosynthetic membranes (PM) [9], and regulatory pathways of the nitrogen fixation system [10,11]. Therefore, it is not surprising that there are currently 10 genome assemblies of several R. rubrum in the GenBank database (accessed March 8th, 2024), including the genome presented in this study. The genes involved in various metabolic pathways are therefore known. Apart from synthetic pathways mentioned above, PHA synthesis occurs, and the three genes coding PHA synthases were identified in R. rubrum genome [12]. While one of these genes is as part of biosynthetic phaCAB operon, the other two are located separately in different locations in the genome. The main challenges related to R. rubrum and its capacity to produce value-added chemicals are a low specific growth rate, a low biomass and PHA volumetric productivity, and the need for an organic co-substrate to increase productivity of the autotrophic pathways in the case of PHAs synthesis. These hurdles could be overcome by adopting genome engineering strategies, such as the one already applied to R. rubrum, for example, consisting of overexpressing genes coding PHA synthases [13]. Despite relatively good genome characterization and even availability of various mutant strains, only a little is known on gene regulations in R. rubrum as studies exploring gene expression on a genomewide scale using RNA sequencing (RNA-Seq) are missing and only two studies based on microarrays are available [14,15]. Thus, in our study, we compared R. rubrum transcriptomes from cultivation on fructose and on acetate to describe basic changes in various pathways observed along the cultivation or between substrates. Moreover, to explore the stability of expression in engineered strains, we also analyzed transcriptomes of a mutant strain with two PHA synthases deleted and compared it to the wild-type strain for both same substrates. Furthermore, as the genome assembly of the type strain was created relatively long time ago, we re-sequenced the genome of R. rubrum DSM 467 T , which was used for experiments to exert influence of possible mutations that might be accumulated over time. We also sequenced a ΔphaC1ΔphaC2 mutant strain to confirm deletions of PHA synthases and to capture other genomic changes. Additionally, we improved genome annotation by small RNA (sRNA) gene inference using RNA-Seq data. 2. Material and methods 2.1. Growth conditions and experiments A freeze-dried bacterial culture of the type strain Rhodospirillum rubrum DSM 467 T (WT strain) was purchased from the Leibnitz Institute DSMZ-German Collection of Microorganism and Cell Cultures, Braunschweig, Germany. The double mutant R. rubrum ΔphaC1ΔphaC2 knock-out (KO) strain was obtained from the team of Professor Kevin E. O΄Connor and Professor Tanja Narancic (University College Dublin, Ireland). This PHA-negative mutant was designed and constructed as reported previously in [16]. The first cultivation step was incubation of both wild-type and mutant strains on LB broth (tryptone 10.0 g/L, yeast extract 5.0 g/L, NaCl 5.0 g/L) in Petri dishes at 30 ◦C in the dark for five days. The second part of cultivation was inoculated by loop of the bacterial culture from Petri dishes and the cultivation was performed in 500 mL Erlenmeyer flasks containing 100 mL of LB broth at 30 ◦C in the dark under shaking at 160 rpm for approximately 72 h till OD 660 =1.5. During the main part of the cultivation, the cultures were inoculated to OD 660 =0.1 by culture grown in liquid LB broth and cultivated in SYN medium. Its composition for 1 L of medium was: 250 mg of MgSO 4 ⋅7 H 2 O, 132 mg of CaCl 2 ⋅2 H 2 O, 10 g of NH 4 Cl, 21 g of MOPS buffer, 10 mL NiCl 2 (20 µM), 100 mL of a chelated iron-molybdenum solution (0.28 g H 3 BO 3 , 2.1 g of Na 2 EDTA, 0.7 g of FeSO 4 ⋅7 H 2 O and 0.1 g of Na 2 MoO 4 per Liter of distilled water). In this medium, two carbon sources were used: 2.4 mL of 1.5 M fructose (18 mM in total) with 4 mL of 50 g/L yeast extract (1 g/ L in total) and 11 mL of 1 M acetate (55 mM in total) with 4 mL of 50 g/L yeast extract (1 g/L in total). They will be referred to as FY and AY, respectively. Cultivation was performed in triplicates in 200 mL of medium in 1 L Erlenmeyer flask at 30 ◦C in the dark under 160 rpm until the last sample of culture in stationary phase was taken. All the cultivations were performed under aerobic conditions. During the cultivation of wild-type R. rubrum (WT) and R. rubrum knock-out (KO), samples were taken at three time-points, at different growth phases. In each sampling time-point, the biological triplicates were spectrophotometrically screened at λ =660 nm and used for further DNA and RNA analysis. The culture-specific growth rates and doubling time were calculated using optical density data. The OD660 measurements for calculating μ max were taken every hour. The μ max values were calculated using a standard method: OD values were transformed into ln(OD660), and the linear part of the ln(OD660) vs. time curve was used for regression analysis to determine μ max . 2.2. Transmission electron microscopy Cultures of wild-type R. rubrum cultivated on both fructose and acetate as carbon sources were fixed using the high-pressure freezing method. Samples were pipetted on the 200 µm side of 3 mm copper-gold carrier type A, which was closed using the flat side of the type B carrier. Both carriers were pretreated with 1 % solution of lecithin in chloroform. Vitrification of the samples was performed using a high-pressure freezer EM ICE (Leica Microsystems, Austria). Frozen samples were transferred under liquid nitrogen into a freeze substitution unit (EM AFS2, Leica Microsystems, Austria). The substitution solution contained 1.5 % OsO 4 in acetone and the protocol was set as previously described in [17]. Freeze substitution was followed by resin embedding (Epoxy Embedding Medium kit, Sigma Aldrich, Germany) and curing at 62 ◦C for 48 h. Cured samples were cut to ultrathin sections (~75 nm) using a diamond knife (ultra 45◦, DiATOME, Switzerland) and ultramicrotome (EM UC7, Leica Microsystems, Vienna, Austria). Since ultrathin sections were imaged using a low-voltage transmission electron microscope (LVEM 25 Delong Instruments, Czech Republic) at 25 kV voltage of electron beam, no post-staining procedure was necessary to achieve sufficient contrast [18]. 2.3. DNA and RNA extraction and sequencing High molecular weight genomic DNA of the WT strain for long-read sequencing was extracted using a MagAttract HMW DNA kit (Qiagen, NL) in accordance with the manufacturer’s instructions. The DNA purity was checked using NanoDrop (Thermo Fisher Scientific, USA), while the concentration was measured using Qubit 4.0 Fluorometer (Thermo Fisher Scientific, USA), and the length was measured using Agilent 4200 TapeStation (Agilent Technologies, USA). The Ligation Sequencing 1D Kit (Oxford Nanopore Technologies, UK) was used for library preparation, and the sequencing was performed using R9.4.1 Flow Cell on the Oxford Nanopore Technologies (ONT) MinION platform. Genomic DNA of WT and KO strains for the high-throughput shortread sequencing was extracted using GenElute Bacterial Genomic DNA Kit (Sigma-Aldrich, USA) in accordance with the manufacturer’s instructions. Sequencing libraries were prepared using the KAPA HyperPlus kit, and sequencing was carried out using MiSeq Reagent kit v2 (500 cycles) on the MiSeq platform (Illumina, USA). RNA isolation followed an optimized extraction protocol consisting of a combination of procedures focused on RNA isolation, where the crucial point was addition of 1 mL of TRIzol per 40 mg of wet biomass followed by incubation for 5 min to permit complete dissociation of the nucleoproteins complex. Then, 0.1 mL of 1-bromo-3-chloropropane per 1 mL of TRIzol™ Reagent used for cell lysis was added, and the tube was securely capped and incubated for 2–3 min at room temperature. Afterwards, the samples were centrifuged at 11,000 ×g at 4 ◦C for 15 min. Subsequently, the supernatant containing the RNA was transferred to a new tube where 70 v/v % EtOH was added at a ratio 1:1. Then, the K. Jureckova et al. Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 2683 samples were transferred to spin columns and the procedure continued according to the manual of the NucleoSpin RNA Plus isolation kit (Macherey-Nagel) with washing and drying steps of silica membrane and elution of RNA that were stored at −80 ◦C till the sequencing. Ribodepletion was performed with QIAseq FastSelect –5S/16S/23S Kit (Qiagen, NL) for WT samples cultivated on fructose (no. 1 – 9, see Supplementary Table S1) or with RiboCop rRNA Depletion Kit for Bacteria Mixed bacterial samples (Lexogen, AT) for remaining samples (no. 10 – 33). Strand-specific sequencing libraries were prepared with NEBNext Ultra II Directional RNA Library Prep Kit (New England Biolabs, USA) and sequenced with Illumina NextSeq550 to produce reversely stranded reads. For samples no. 10 – 33, Unique Molecular Identifiers (UMIs) were added using xGen Duplex Seq Adapters (IDT, USA). 2.4. Genome assembly Nanopore sequencing data of wild-type strain were basecalled using Guppy (v3.4.4), and the data quality was checked using PycoQC (v2.2.3, [19]). De novo assembly by Flye (v2.8.1, [20]) was performed, and the obtained contigs were polished using Minimap2 (v2.17, [21]) combined with Racon (v1.4.13, [22]) and then final polishing was performed by Medaka (v1.1.2). Both WT and KO strains were sequenced using Illumina MiSeq. Obtained reads were firstly checked for quality using FastQC (v0.11.5) combined with MultiQC (v1.7, [23]) and secondly, the low-quality ends of reads together with sequencing adapters were removed by Trimmomatic (v0.36, [24]) with following settings: ILLUMINACLIP:TruSeq3-PE. fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36. Next, the set of reads was checked for contamination by human DNA, and detected reads were removed using BBMap (v39.01) and human genome GCF_000001405.40 available in the NCBI RefSeq database. In the case of the WT strain, the reads were mapped to the nanopore contigs using BWA (v0.7.17, [25]), and for managing the obtained assembly Samtools (v1.10, [26]) was employed. Final assembly polishing was made by Pilon (v1.24, [27]). In the last step, the assembled genome and plasmid were rearranged so that the DnaA gene was the first gene in the genome and repB was the first gene in the plasmid. For the KO strain assembly, BWA and Samtools were used again, but in this case, the pre-processed reads were mapped to the previously assembled WT genome. The variant calling was conducted to find differences between both (WT and KO) analyzed strains. For this purpose, GATK4 (v4.3) was used. Underrepresented variants, variants with low coverage and false positive calls were filtered out, and the remaining ones underwent additional analysis to ascertain their presence within coding regions and, if applicable, their impact on the phenotype. 2.5. Genome annotation NCBI Prokaryotic Genome Annotation Pipeline (PGAP) [28] was used for the wild-type strain chromosome and plasmid annotation. Protein coding genes functional annotation was performed by classifying them into clusters of orthologous groups (COG) categories from the eggNOG database via the eggNOGmapper (v2.1.6, [29]). The DNAplotter [30], which is integrated into Artemis (v18.2.0, [31]), was used to create the chromosome and plasmid circular maps including GC content and CG skew plots, which were calculated in a window of size 10,000 bp with 200 bp step. The interspaced short palindromic repeats (CRISPR) arrays were searched by the CRISPDetect tool (v2.4, [32]), and Cas genes were manually searched in the annotated genome. Physpy (v4.2.6, [33]) and online tools Prophage Hunter and Phaster were used for prophage DNA identification. The restriction-modification systems were located using REBASE database [34]. The GenBank database was searched for genomes of R. rubrum using Entrez [35]. Then Roary (v3.13, [36]) was used to identify the core genome, which consists of genes located in all analyzed strains, the accessory genome formed by genes presented in at least two strains but not in all of them, and the number of unique genes. The minimum percentage identity to assess whether two genes are similar was 95 %, and all other parameters for Roary were left at their default settings. 2.6. Transcriptomic analysis The raw RNA-Seq reads for both strains were checked for their quality using FastQC and MultiQC. Next, the reads trimming was performed to discard low-quality bases and adapters using Trimmomatic with following parameters for samples no. 10 – 33: ILLUMINACLIP: TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 and for samples no. 1 – 9 HEADCROP parameter was added with value 5, because first five bases of the reads contained randomized 5 bp long adapters. Remaining contamination by rRNA was removed using SortMeRNA (v4.3.4, [37]) together with the default database (smr_v4.3_default_db.fasta) and the SILVA database [38] with known 16S and 23S rRNA sequences. Processed reads were mapped to the genome of the wild-type strain using STAR (v2.7.10a, [39]), and mapped reads were deduplicated using UMI-tools (v1.1.4, [40]). Finally, the reads were counted using featureCounts function [41] from Rsubread package (R/Bioconductor). Reads counting considered two options: Uniquely mapped reads and multimapping reads. For the multimapping reads the contribution of those reads to the final count was always divided by the number of genomic loci to which the read was mapped, therefore the number of reads remained unchanged. Created count tables for all samples were further normalized by calculating RPKM (Reads Per Kilobase per Million mapped reads) and by using built-in function in R/Bioconductor package DESeq2 [42], which was also used for differential expression analysis. Normalized count tables were used for dimension reduction using Barnes-Hut t-SNE implemented in R package Rtsne [43]. The results were visualized using ggplot2 [44] R package, which was also employed for creation of Volcano plots. Gene ontology (GO) enrichment analysis was performed by topGO [45] R/Bioconductor package together with GO map, that was created with custom Python script based on available GO annotation of wild-type strain in NCBI RefSeq database (NZ_CP077803.1 for chromosome and NZ_CP077804.1 for plasmid). Finally, the expression profiles of selected genes were visualized using heatmaps created using R packages gplots, RColorBrewer, magick and openxlsx. To gain insight into the regulatory processes of R. rubrum, mapped reads were further used for non-coding RNAs (ncRNA) prediction with baerhunter (v0.9, [46]) The sample depths were normalized using sizeFactors function from DESeq2 package and visualized using boxplots. Subsequently, to obtain a count table of a collapsed annotation file with newly predicted features, barhunter’s count_features script using featureCounts function was applied. Furthermore, attention was paid to putative sRNAs, and their length distribution was visualized by a histogram. Also, these predictions were categorized as trans, respectively cis-acting elements using IRanges R package [47]. Finally, differential expression analysis was performed on the count table using DESeq2. Counts of the differentially expressed (up or down) sRNAs between specific conditions were visualized by a bar chart using R package ggplot2 as well as for all mentioned graphs. 3. Results 3.1. Genome assembly The Oxford Nanopore Technologies MinION produced 135,804 reads; from them, 104,955 had Q >7 and were used for further WT strain assembly. The mean read’s length was approximately 3.5 kbp. The Illumina MiSeq provided 2.5 million 250 bp-long paired reads with an average Phred score of 35. From them, 296 reads were mapped to the human genome; thus, they were discarded from further analysis. The K. Jureckova et al. Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 2684 assembly process of wild-type strain resulted in the final assembly of one circular chromosome and one circular plasmid with a coverage of 370×. The sequences were deposited under accession numbers CP077803.1 for chromosome and CP077804.1 for plasmid at DDBJ/EMBL/GenBank. The obtained chromosome sequence length was 4,352,570 bp with a GC content of about 65.4 %, and the plasmid was 53,835 bp long with a GC content of around 59.8 %. In total, 3,968 open reading frames (ORFs) divided into 2,146 operons were identified for both sequences. Most of the genes were protein coding; however, 49 pseudogenes were also found. Furthermore, 521 genes overlapped with another neighboring gene, and 3 overlaps were found between a gene and a pseudogene. The overlap size was in 390 cases equal to three nucleotides and the longest overlap was 87 bp between genes KUL73_19715 and KUL73_19720 on plasmid. See Table 1 for complete statistics for chromosome and plasmid. Illumina MiSeq sequencing of KO strain provided about 2.2 million 250-bp long reads with an average Phred score of 33. Of these, 260 reads that mapped to the human genome were removed. Variant calling of the KO strain confirmed deletion of PHA polymerases. Moreover, it revealed three changes in its genome: two single nucleotide mutations and one insertion. The first mutation was found in position 823,861 (A → G) in the promoter of cimA gene encoding citramalate synthase, the second mutation was identified in position 1,305,962 (C → T) in the rpoH gene encoding RHA polymerase sigma factor RpoH and the last third change was localized in position 3,782,354 (C → CGCTTCAGGGGAAACACGTTATGAAG) in the promoter of gene encoding hypothetical protein. Ten genomes obtained from GenBank were used for R. rubrum core genome determination. In total, 1,732 genes were found in all analyzed strains and thus formed the core genome. Another 2,859 genes were identified in at least two genomes, and thus, they comprised accessory genome. Together these 4,591 genes form the R. rubrum pangenome. In addition, 920 unique genes only identified in one genome were found. The number of unique genes ranged from 4 to 467, with a median of 44 (see Supplementary Table S2). 3.2. Wild-type strain functional annotation Protein coding genes and pseudogenes were classified according to COG into 20 categories. Only 2 CDSs were not assigned to any COG; however, 1,093 CDSs were classified to class S with unknown function. The remaining 2,804 CDSs (out of all 3,850 protein coding genes and 49 pseudogenes) were classified to COG class. Individual chromosomal and plasmidic features are shown in Fig. 1 together with GC content and GC skew plots for the whole genome. Substantial drop in the average CG content was observed around position 3,810,000 bp on chromosome, where are located 5S rRNA (KUL73_17050), 23S rRNA (KUL73_17055) and 16S rRNA (KUL73_17070) genes. More detailed results can be seen in Supplementary Table S3. The genome of R. rubrum DSM 467 was annotated in terms of gene ontology. Together, 3,047 GO terms were assigned to 1,312 genomic loci on the chromosome, and 30 GO terms were connected to 14 genomic elements on the plasmid. The most common terms for the three GO categories were: GO:0006412 “translation” for biological process assigned 58 times, for molecular function GO:0005524 “ATP binding” with 86 genomic loci and for cellular component GO:0016020 “membrane” was assigned 107 times. CRISPR analysis showed 11 arrays with lengths from 337 to 2,467 bp and a number of spacers from 4 to 40 (see Supplementary Table S4). CRISPRDetect tool did not identify any Cas genes, but a manual search of the annotated genes found 34 Cas-like genes, but no Cas9 protein. In the case of prophage identification, Physpy did not identify any viral DNA, Prophage Hunter found one ambiguous prophage candidate, and Phaster localized four incomplete prophages, all of them presenting a low score (see Supplementary Table S5). The restriction-modification systems analysis revealed two systems containing restriction endonuclease and methylase where one belonged to type II and one to type III. In addition, another gene coding for methylase type II was also found (see Supplementary Table S6). 3.3. Cultivation & growth kinetics To achieve cultivation-driven transcriptome changes in the wild-type strain of R. rubrum and its knock-out strain, these microorganisms were cultivated on two substrates, namely acetate and fructose. Growth curves were determined throughout cultivation, with samples taken over time and their optical density measured at 660 nm. Except for the cultivation of the KO strain on acetate, samples were collected in each cultivation in the mid-exponential phase, at the end of the exponential phase, and in the stationary phase (for precise sampling times see Supplementary Table S1). For the KO strain on acetate, only two time points were selected for further characterization (see Fig. 2). From the obtained growth curves, the growth rate and doubling time for each cultivation were also calculated (see Table 2 and Supplementary Fig. S1). Additionally, pH was of course also analyzed during the cultivations. The pH of both cultures cultivated on fructose quickly reached values of 7.5 - 7.7 and remained almost constant until the end of cultivation. In the case of acetate, the pH reached more alkaline values of about 8.6 - 8.8 during the initial stage of cultivation and remained constant until the end of the experiment. 3.4. Ultrastructural analysis The presence of cytoplasmic PHA granules in the cultures cultivated on two different carbon sources as well as the overall ultrastructure of the cells of R. rubrum wild-type was determined using low voltage transmission electron microscopy. As seen in Fig. 3, both cultures contain in their spiral-shaped cells a substantial amount of small electron-lucent granules – chromatophores. However, only cultures grown on acetate as a carbon source were able to also produce PHA in their cells, which formed bigger electron-lucent granules. 3.5. RNA-Seq transcriptome RNA sequencing (Supplementary Table S1) produced 713 millions of single-end reads averaging on 21.6 million reads per sample. After the first part of pre-processing (adapter and quality trimming), the reads’ lengths were either 70 bp (samples no. 1 – 9) or 66 bp (samples no. 10 – 33) based on the used library construction method, and the average Phred score was 35 for all samples. The subsequent pre-processing step removed the remaining contamination of rRNA from the data. The proportion of reads corresponding to 16S and 23S differed among samples, e.g., all three samples for WT strain grown on fructose in the third time-point contained around 46 % of rRNA reads, and on the other hand 23 samples had less than 5 % of rRNA contamination (see Supplementary Fig. S2). Finally, reads were mapped to the reference genome of the wild-type strain. The majority of reads (81 - 99 %) were mapped uniquely, yet Table 1 Genomic features of Rhodospirillum rubrum DSM 467 T . Chromosome Plasmid Length [bp] 4,352,570 53,835 GC content [%] 65.4 59.8 Number of ORFs 3,919 49 Number of operons 2,116 30 Protein coding genes 3,807 43 Pseudogenes 43 6 rRNA genes (5S, 16S, 23S) 4, 4, 4 - tRNA 54 - ncRNA 3 - K. Jureckova et al. Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 2685 multimapping reads, persisting deduplication, were also considered due to the existence of overlapping genes in the genome (see Supplementary Fig. S3). The multimapping reads were in further analysis down heightened by the number of associated genomic features; therefore, the original number of reads in each sample remained the same. Overall, only two pseudogenes located on the plasmid (KUL73_19640 and KUL73_19670) remained completely silent with RPKM <1 in all tested conditions. Reproducibility of our experiments was verified by three biological replicates for each sampling time-point and all tested condition, and visualized through dimension reduction of normalized read counts by tDistributed Stochastic Neighbor Embedding (t-SNE) method. In most cases, the data points formed individual clusters representing specific cultivation condition and sampling time-point (see Fig. 4). Exceptions can be found in WT AY and KO AY samples, where two bigger clusters were formed based on the type of cultivated strain and used substrate, yet within them smaller clusters can be located and they are representing the time-points. The main differences between cultivation conditions were identified through differential expression analysis, which was performed between cultivations on fructose and on acetate separately for wild-type (WT AY vs WT FY) and for knock-out strain (KO AY vs KO FY), and between wildtype and knock-out strain, separately for cultivation on fructose (KO FY vs WT FY) and on acetate (KO AY vs WT AY). Therefore, four different combinations of cultivation conditions were tested and for each of them, comparisons between available time-points were conducted, e.g., for WT AY vs WT FY differences were found between WT AY T1 vs WT FY T1, WT AY T2 vs WT FY T2, and WT AY T3 vs WT FY T3, and similarly for the rest of comparisons (see Supplementary Figs. S4-S13 for corresponding Volcano plots). Results for each combination of cultivation Fig. 1. Chromosomal maps of Rhodospirillum rubrum DSM 467 T chromosome and plasmid. The first, second, and third outermost circles represent CDSs on the forward and backward strands, and pseudogenes, respectively. Classification of COGs is represented by colors. Next, RNA genes, distinguishing among tRNA, rRNA, and ncRNA, are represented in the fourth outermost circle. The inner area represents the GC content and GC skew (window size 10,000 bp, step size 200 bp). Fig. 2. Growth dynamics of R. rubrum WT and KO on acetate and fructose, obtained in SYN medium at 30 ◦C, 160 rpm in dark. Table 2 Growth rate and doubling time for R. rubrum WT and KO strains grown on acetate (AY) and fructose (FY). Sample µ max [h −1 ] Td [h] WT FY 0.057 ±0.001 12.082 ±0.280 WT AY 0.094 ±0.000 7.356 ±0.039 KO FY 0.037 ±0.002 18.992 ±1.299 KO AY 0.017 ±0.001 40.312 ±1.946 K. Jureckova et al. Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 2686 conditions were then separately searched for genes with any statistically significant changes in the gene expression levels (either upor downregulated) (adjusted p value <0.05, Benjamini-Hochberg correction) and they were sorted into two groups based on the number of significant regulations. The first group contained genes with at least one regulation, i.e., statistically significant differential expression, among tested pairs of time-points. The second group consisted of genes that had statistically significant change in every performed comparison of available timepoints. The first group was used as gene universum in GO enrichment analysis and the second group was used as a set of interesting genes, that Fig. 3. Morphology of R. rubrum wild-type grown on A) acetate or B) fructose as carbon source, imaged using low voltage transmission electron microscopy. PHA granules are marked with arrows, smaller electron-lucent granules represent chromophores. Fig. 4. Comparison of RNA-Seq samples through dimensionality reduction by t-SNE. All samples are represented as points color-coded according to strain (WT or KO), substrate (FY or AY), time-point (T1, T2 or T3) and text label indicates biological replicate (sfA, sfB or sfC). K. Jureckova et al. Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 2687 were used to identify GO terms that were significantly enriched (p value <0.05, Fisher’s exact test) between tested conditions. GO enrichment analysis revealed that for the comparison WT AY vs WT FY 18 GO terms were enriched in biological process (BP) category and 11 terms in molecular function (MF) category. BP terms were related to “disaccharide metabolic process”, “cell cycle”, “DNA recombination” and more, and MF terms corresponded to “acetyltransferase activity”, “transposase activity” or “isomerase activity” (see Supplementary File 1 and Table 3 for results of differential expression analysis of selected genes). Comparison between KO AY and KO FY showed 10 enriched BP GO terms connected to “regulation of biological process”, “maintenance of DNA repeat elements” and “organic acid catabolic process” terms and the only MF term “intramolecular transferase activity” (see Supplementary File 2). The difference between WT and KO strains grown on fructose where characterized with 12 BP terms and 5 MF terms, that where significantly enriched. The main BP terms described “archaeal or bacterial-type flagellum-dependent cell motility” or “vitamin metabolic process” processes and MF terms were related to “NAD binding”, “methyltransferase activity” or “oxidoreductase activity, acting on CH-OH group of donors” (see Supplementary File 3 and Table 4 for results of differential expression analysis of selected genes). Finally, GO enrichment analysis between KO AY and WT AY identified 21 BP and 6 MF enriched terms. The main BP terms belonged to “aspartate family amino acid metabolic process”, “regulation of biological process”, “regulation of cellular metabolic process” terms. MF terms were “pyrophosphatase activity” or “hydrolase activity, acting on acid anhydride” (see Supplementary File 4 and Table 5 for results of differential expression analysis of selected genes). 3.6. Small RNA prediction For the prediction of small RNAs from the strand-specific RNA-Seq data, a coverage-based detection with baerhunter tool, which requires setting three input values, was performed. While the low_cut_off threshold was left on default value of 5, the high_cut_off threshold was set to 25 as inferred from our dataset (see Supplementary Fig. S14). The last parameter, min _sRNA_length of predicted elements was set to 40 bp. The length distribution of putative sRNAs (see Supplementary Fig. S15) contains a wide variety of lengths ranging from ten to thousands bp. Table 6 then contains counts of predicted elements in the chromosome and plasmid sequences of R. rubrum. Small RNAs were further classified into two groups based on their overlap with other CDSs. Those predictions that did not overlap with the CDS on either strand, thus were in intergenic regions, were labeled as trans-encoded sRNAs. On the contrary, if there was even a slight overlap with any annotated element on the opposite strand as an inferred sRNA, the prediction was considered cis-encoded. Additionally, differential expression analysis revealed that majority (>99 %) of predicted sRNA were at least once differentially expressed among various combinations of available conditions, see Table 6. The Table 3 Differential expression analysis results of selected genes related to significantly enriched GO terms for comparison between cultivation on fructose and on acetate for wild-type strain (WT AY vs WT FY). WT AY T1 vs WT FY T1 WT AY T2 vs WT FY T2 WT AY T3 vs WT FY T3 Gene abbr. Putative physiological function Locus tag log2 fold change p-adj log2 fold change p-adj log2 fold change p-adj MF GO term: GO:0009055 electron transfer activity nuoB NADH-quinone oxidoreductase subunit NuoB KUL73_07380 -1.65 4.95E02 -1.52 5.11E02 -1.04 2.43E01 NADH-quinone oxidoreductase subunit C KUL73_08070 -0.33 3.31E02 0.54 2.38E04 -1.08 2.95E13 nuoF NADH-quinone oxidoreductase subunit NuoF KUL73_08085 0.66 6.42E04 2.15 5.34E31 1.51 3.40E15 NADH-quinone oxidoreductase subunit J KUL73_08105 -0.84 7.79E03 1.00 1.18E03 3.49 4.78E25 nuoL NADH-quinone oxidoreductase subunit L KUL73_08115 -1.36 2.39E09 0.61 8.42E03 1.98 2.17E17 NADH-quinone oxidoreductase subunit M KUL73_08120 -2.27 7.70E22 -0.35 1.58E01 1.23 6.61E07 nuoN NADH-quinone oxidoreductase subunit NuoN KUL73_08125 -1.12 1.87E06 0.87 1.75E04 2.02 6.22E17 grxD Grx4 family monothiol glutaredoxin KUL73_03670 -1.67 1.44E06 -3.93 1.46E31 -2.82 7.84E17 c-type cytochrome KUL73_05315 -1.51 1.16E12 -1.72 1.44E16 -2.66 1.14E37 petA ubiquinol-cytochrome c reductase iron-sulfur subunit KUL73_06225 0.14 7.21E01 -1.77 2.11E08 -2.18 9.80E12 cytochrome b/b6 domain-containing protein KUL73_06485 -0.55 8.36E02 -2.70 5.84E20 -2.79 4.51E21 ccoP cytochrome-c oxidase, cbb3-type subunit III KUL73_17205 -0.74 7.63E03 1.55 6.66E09 1.93 3.56E12 ccoO cytochrome-c oxidase, cbb3-type subunit II KUL73_17215 -1.91 1.81E09 1.25 8.00E05 1.03 1.97E03 ccoN cytochrome-c oxidase, cbb3-type subunit I KUL73_17220 -0.79 9.68E03 1.73 3.21E09 1.90 3.12E10 pseudoazurin KUL73_05930 -0.51 1.35E01 3.25 1.64E22 1.42 1.02E03 BP GO term: GO:0006950 response to stress ahpC peroxiredoxin KUL73_07360 -1.56 2.05E06 -1.80 2.39E08 -1.56 1.57E06 msrA peptide-methionine (S)-S-oxide reductase MsrA KUL73_12215 -1.64 1.64E04 -5.21 1.09E35 -4.50 5.38E27 nth endonuclease III KUL73_00795 -1.41 8.15E10 -3.08 5.56E43 -2.03 2.29E16 K. Jureckova et al. Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 2688 further statistics showing counts of differentially expressed sRNAs between these various combinations also distinguishes whether these elements are up/down regulated (see Supplementary Fig. S16). 4. Discussion 4.1. Genome and transcriptome R. rubrum DSM 467 T hybrid assembly combining long Oxford Table 4 Differential expression analysis results of selected genes related to significantly enriched GO terms for comparison between wild-type vs knock-out strain for cultivation on fructose (KO FY vs WT FY). KO FY T1 vs WT FY T1 KO FY T2 vs WT FY T2 KO FY T3 vs WT FY T3 Gene abbr. Putative physiological function Locus tag log2 fold change p-adj log2 fold change p-adj log2 fold change p-adj BP GO term: GO:0009110 vitamin biosynthetic process thiE thiamine phosphate synthase KUL73_05625 -1.19 1.32E06 -0.44 8.13E02 0.88 4.37E03 thiD bifunctional hydroxymethylpyrimidine kinase/ phosphomethylpyrimidine kinase KUL73_05735 0.09 4.91E01 0.75 5.81E10 -0.74 2.96E04 thiL thiamine-phosphate kinase KUL73_09425 1.17 1.44E07 3.08 5.67E43 1.23 1.96E04 thiC phosphomethylpyrimidine synthase ThiC KUL73_10385 1.06 1.49E04 0.34 2.45E01 -0.65 3.02E02 pdxY pyridoxal kinase KUL73_06260 1.05 4.41E04 1.91 5.21E11 1.74 3.95E07 pyridoxine 5 ′ -phosphate synthase KUL73_09595 1.05 4.41E04 1.91 5.21E11 1.74 3.95E07 pdxH pyridoxamine 5 ′ -phosphate oxidase KUL73_13555 0.05 8.79E01 0.89 5.05E04 -0.71 1.15E02 cobS cobaltochelatase subunit CobS KUL73_01095 -0.40 3.57E03 -0.79 3.57E09 -0.97 2.89E11 cobN cobaltochelatase subunit CobN KUL73_17390 0.46 8.99E02 2.40 5.11E21 1.07 1.19E04 cobS adenosylcobinamide-GDP ribazoletransferase KUL73_03465 1.58 1.63E08 2.38 5.65E18 1.42 1.66E04 cobU bifunctional adenosylcobinamide kinase/adenosylcobinamidephosphate guanylyltransferase KUL73_03475 0.42 5.80E02 2.44 7.94E29 0.52 1.77E01 cobW cobalamin biosynthesis protein CobW KUL73_03480 0.33 3.57E02 1.68 1.83E29 -0.93 1.43E08 cobF precorrin-6A synthase (deacetylating) KUL73_15320 -0.13 6.83E01 0.73 1.12E02 0.22 6.50E01 cobalt-precorrin-6A reductase KUL73_15410 0.36 1.54E01 0.90 2.06E04 -0.79 1.60E02 precorrin-2 C(20)-methyltransferase KUL73_15420 1.33 2.04E06 2.30 1.57E16 2.27 4.18E08 cobJ precorrin-3B C(17)-methyltransferase KUL73_15415 1.49 1.35E17 2.75 2.30E53 2.19 2.69E13 BP GO term: GO:0071973 bacterial-type flagellum-dependent cell motility BP GO term: GO:0097588 archaeal or bacterial-type flagellum-dependent cell motility fliJ flagellar export protein FliJ KUL73_02755 0.84 5.36E03 1.88 1.39E10 1.28 2.55E04 flagellin KUL73_13080 -1.35 1.00E05 -2.56 5.22E18 -2.46 2.30E16 flagellar basal body-associated FliL family protein KUL73_14640 1.06 2.11E02 1.13 1.12E02 2.72 1.03E09 flgG flagellar basal-body rod protein FlgG KUL73_14650 -5.77 2.36E33 -7.48 2.46E57 -1.41 5.68E03 flagellar basal body L-ring protein FlgH KUL73_14660 -1.21 5.01E03 -1.41 8.07E04 1.12 1.26E02 hypothetical protein KUL73_14700 -0.97 1.51E03 -3.52 2.04E33 -0.18 5.87E01 flagellin KUL73_14725 -1.46 7.41E39 -2.47 1.91E110 -2.37 3.48E101 MF GO term: GO:0016614 oxidoreductase activity, acting on CH-OH group of donors MF GO term: GO:0016616 oxidoreductase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor guaB IMP dehydrogenase KUL73_01290 -1.21 1.53E20 -0.64 6.71E07 -0.99 5.12E11 NADP-dependent isocitrate dehydrogenase KUL73_01875 -2.10 3.10E50 -1.24 1.05E18 -1.75 1.80E28 mdh malate dehydrogenase KUL73_06315 -1.13 2.58E06 -0.72 2.82E03 -0.04 8.98E01 3-hydroxybutyryl-CoA dehydrogenase KUL73_15880 -0.68 4.10E04 -1.41 4.66E14 -0.92 1.85E06 MF GO term: GO:0051287 NAD binding NADP-dependent isocitrate dehydrogenase KUL73_01875 -2.10 3.10E50 -1.24 1.05E18 -1.75 1.80E28 nuoF NADH-quinone oxidoreductase subunit NuoF KUL73_08085 0.41 3.86E02 2.17 2.66E31 0.45 2.96E02 K. Jureckova et al. Computational and Structural Biotechnology Journal 23 (2024) 2681–2694 2689 Nanopore and short Illumina reads reveal one circular chromosome and one plasmid contig as was expected based on the previously published type strain genome [10]. The chromosome exhibits GC content of 65.4 %, and the plasmid has 59.8 %, which is more than the average for Gram-negative bacteria [48]; nevertheless, it again corresponds to the previously published type strain as well as the chromosome length 4.35 Mbp and plasmid length 53.84 kbp. In comparison with type strain, the overall number of genes is similar; there is only a slight difference between the number of predicted protein-coding genes (DSM 467 T :3807 vs S1 T :3850), number of RNAs (DSM 467 T :71 vs S1 T :83) and predicted pseudogenes (DSM 467 T :49 vs S1 T :9). The differences may be due to the use of different sequencing platforms for strain sequencing and various tools for their assemblies and annotations. RNA-Seq samples for WT cultivation on fructose (samples no. 1 – 9) were prepared slightly differently from the remaining ones (samples no. 10 – 33). Particularly, different rRNA depletion kit was used for the latter group, which resulted in very low contamination of reads belonging to 16S and 23S rRNA genes (see Supplementary Fig. S2) and UMIs were added to allow their deduplication which caused the difference in read lengths after pre-processing. While the reads without UMIs collected from the WT cultivation on fructose were 70 bp long, the remaining ones were only 66 bp long. Yet, these dissimilarities did not cause any significant differences in the overall high number of reads per sample or in the quality of the reads. Furthermore, reads were mapped to the reference wild-type strain and the results also proved a high quality of our data as we were able to uniquely map the majority of the reads. In the worst case, WT_FY_sfB_T3 sample had 81 % of uniquely mapped reads. This particular sample also had the highest contamination of rRNA in reads with almost 47 %. On the other hand, many samples had over 98 % of the reads with unique hit within the genome, and the multimapping reads were identified only in 1–2 % of the reads. Multimapping reads could be a consequence of 523 existing overlaps between neighboring genes within the genome of R. rubrum (see Supplementary Fig. S3). Nevertheless, our RNA-Seq data proved to be of a very high quality, thus enabling further genome-wide study of changes in the expression levels between different cultivation conditions. As main differences were observed for particular substrates, we firstly used RNASeq data to detect differential expression on a genome-wide scale in particular time-points for acetate and fructose cultivation for WT (Supplementary Figs. S4 – S6) and KO (Supplementary Figs. S7 and S8) strains. Additionally, we compared differences between KO and WT strains in particular time-points on fructose (Supplementary Figs. S9 – S11) and acetate (Supplementary Figs. S12 and S13). All volcano plots showed expected relation among fold changes and statistically significant differences in expression of protein coding genes where number of down-regulated and up-regulated genes is roughly the same. 4.2. Comparative analysis of acetate and fructose growth The GO enrichment analysis comparing acetate and fructose cultures Table 5 Differential expression analysis results of selected genes related to significantly enriched GO terms for comparison between wild-type vs knock-out strain for cultivation on acetate (KO AY vs WT AY). KO AY T1 vs WT AY T1 KO AY T2 vs WT AY T2 Gene abbr. Putative physiological function Locus tag log2 fold change p-adj log2 fold change p-adj MF GO term: GO:0016887 ATP hydrolysis activity bchI magnesium chelatase ATPase subunit I KUL73_02545 -1.00 4.95E04 -0.68 1.98E02 AFG1 family ATPase KUL73_06305 -0.51 2.44E03 -0.28 9.87E02 arsA arsenical pump-driving ATPase KUL73_07505 1.89 4.00E10 1.50 9.17E07 tsaE tRNA (adenosine(37)-N6)-threonylcarbamoyltransferase complex ATPase subunit type 1 TsaE KUL73_17750 1.32 1.32E05 1.58 2.37E07 AAA family ATPase KUL73_19335 -0.63 3.15E05 -0.68 7.10E06 mfd transcription-repair coupling factor KUL73_08940 -0.96 6.99E10 -1.07 7.44E12 uvrA excinuclease ABC subunit UvrA KUL73_09080 1.09 1.59E04 1.53 1.04E07 mutL DNA mismatch repair endonuclease MutL KUL73_15195 -1.16 1.75E09 -1.02 1.40E07 htpG molecular chaperone HtpG KUL73_00370 1.06 1.54E02 1.00 2.13E02 groL chaperonin GroEL KUL73_00840 1.71 4.34E04 1.88 1.06E04 clpB ATP-dependent chaperone ClpB KUL73_03910 1.89 3.75E11 2.55 2.95E19 dnaK molecular chaperone DnaK KUL73_18350 0.77 1.04E02 1.42 1.44E06 lon endopeptidase La KUL73_08040 1.04 3.86E12 1.14 3.20E14 hslU ATP-dependent protease ATPase subunit HslU KUL73_18580 0.76 1.54E05 1.06 9.12E10 Table 6 Number of predicted small RNAs in Rhodospirillum rubrum DSM 467 T . Feature Predicted # in chromosome Differentially expressed Predicted # in plasmid Differentially expressed Total sRNAs 2,329 2,321 10 10 trans 55 55 0 0 cis 2,274 2,266 10 10 K. Jureckova et al.