CC-BY 4.0
Introduction
Biomineralization is a process through which some living organisms form mineral phases. It is found in all the domains of life (Boskey, 2003). In some cases, biomineralization is genetically controlled, i.e. it is tightly controlled by specific genes, as documented in diverse eukaryotic organisms (Brownlee et al., 2015; Gilbert et al., 2022; Knoll, 2003; Marin & Luquet, 2004; Monteiro et al., 2016). For example, key genes encoding, e.g., carbonic anhydrases (Chen et al., 2018; Le Roy et al., 2014), Ca²⁺, H⁺ and HCO3- transporters (Mackinder et al., 2011), and/or molecules inhibiting mineralization (Marin et al., 2000, 1996), play pivotal roles in the control of calcium carbonate (CaCO3) precipitation by eukaryotic organisms. However, there are much fewer documented cases of genetically controlled biomineralization by prokaryotes (Cosmidis & Benzerara, 2022; Görgen et al., 2020). For a long time, the only well documented case was the intracellular formation of magnetites (Fe3O4) and/or greigites (Fe3S4) by magnetotactic bacteria (Monteil et al., 2018). This process is regulated by a suite of more than 40 genes involved in the control of the size, shape and crystallochemistry of the mineral byproducts (Liu et al., 2023; Taoka et al., 2023; Uebe & Schüler, 2016). More recently, the formation of intracellular amorphous calcium carbonates (iACC) by several species of cyanobacteria has been evidenced as an additional case of genetically controlled biomineralization in bacteria (Benzerara et al., 2014; Couradeau et al., 2012). Before this discovery, cyanobacteria were usually thought to induce the precipitation of CaCO3 as an indirect by-product of photosynthesis, with no apparent involvement of specific genes (Altermann et al., 2006). Cyanobacteria biomineralizing iACC are cosmopolitan (Ragon et al., 2014) and can sequester alkaline earth elements such as Ca, Sr, Ba and Ra, including their radioactive isotopes (Blondeau et al., 2018a; Cam et al., 2016; Mehta et al., 2019). As a result, these bacteria may be interesting candidates for bioremediating some radioactive pollutions (Mehta et al., 2022a). Morever, the biomineralization of iACC has been evidenced in phylogenetically diverse strains within the phylum of Cyanobacteria as well as in other bacterial phyla (Liu et al., 2021; Monteil et al., 2021). It is also found in several strains (but not all) of the well-studied Microcystis genus (Gaëtan et al., 2023). This cyanobacterium is particularly interesting since it can form blooms, with some genotypes being toxic for human and wildlife (Harke et al., 2016), causing numerous environmental issues in eutrophic freshwater ecosystems (Dick et al., 2021). Whether the formation of iACC by some Microcystis strains may impact their ecological success under certain conditions remains to be deciphered. Some putative biological functions have been suggested for iACC such as intracellular storage of inorganic carbon, intracellular pH buffering or their use as a ballast modifying cell buoyancy (Cosmidis & Benzerara, 2022). About the later hypothesis, we note that buoyancy is not just controlled by cell density but also by additional parameters such as EPS adhering to the cells, enhancing cell aggregation. For example, Gu et al. (2020) showed that Ca induces EPS production by M. aeruginosa, which can increase buoyancy through the increased size of cell aggregates. The relative contribution of these opposing parameters on the buoyancy should be assessed in the future. Moreover, we can speculate that the fulfilment of these functions may involve some kind of homeostasis which could require both iACC precipitation and dissolution processes, although iACC dissolution has not been evidenced yet in Cyanobacteriota. Last, the formation of iACC, a concentrated intracellular sink of Ca, by the large biomasses of Microcystis during blooms may impact the local geochemical cycle of Ca in these environments, although future studies will have to further document this possibility (Gaëtan et al., 2023). Indeed, the fraction of a natural Microcystis bloom with Ca-concentrating capabilities remains to be determined. Moreover, the dependence of Ca-sequestration by cells on extracellular Ca concentrations needs to be specified. While Ca concentrations may be lower in many lakes than in the BG11 medium (less than 50% of the lakes with a pH higher than 7.4 in the database investigated by Weyhenmeyer et al. 2019), we note, however, that concentrations in lakes vs laboratory cultures cannot be simply compared. Indeed, cultures in the laboratory of iACC-forming cyanobacteria have been conducted so far under batch conditions with a Ca supply limited by the volume of the culture, a condition that may be very different from those encountered in an open system such as a lake.
Despite the potential environmental importance of this process, the molecular mechanisms driving the biomineralization of iACC in cyanobacteria remain poorly known. It has been shown that it costs some energy to the cells (Cam et al., 2018; De Wever et al., 2019). Moreover, a genetic control has been proposed (Benzerara et al., 2022). Indeed, based on a comparative genomics approach, one family of genes, named ccyA, was found in phylogenetically diverse iACC-forming cyanobacteria and was absent from cyanobacteria not biomineralizing iACC. The ccyA gene was found in about one third of all publicly released genomes of Microcystis i.e., 93 out of 282 genomes assemblies available at the time of publication (Gaëtan et al., 2023). The protein encoded by ccyA, called calcyanin, lacks detectable full-length homologs with a known function in public databases. Therefore, its function remains unknown. Several clues suggest that calcyanin may play a role in calcium homeostasis (Benzerara et al., 2023, 2022). However, since the expression of the ccyA gene in mutants of cyanobacteria not forming iACC did not result in the formation of iACC, this biomineralization process likely involves additional genes (Benzerara et al., 2022). For example, some of these genes may play a role in the fact that iACC inclusions form within an intracellular compartment, the envelope of which remains of unknown chemical composition (Blondeau et al., 2018b). Such a compartmentalization may serve to achieve chemical conditions allowing the precipitation and the stabilization of ACC. Additional examples of molecular processes that may favor iACC precipitation are carbon-concentrating mechanisms (CCMs) (Li et al., 2016; Görgen et al., 2020), which have been abundantly described in cyanobacteria and involve diverse genes and multiple mechanisms (Badger & Price, 2003; Kupriyanova et al., 2023). Last, calcium transport genes may also be involved in the accumulation of high amounts of calcium within iACC (De Wever et al., 2019). Consistently, the presence/absence of the ccyA gene in genomes was correlated with the presence/absence of genes involved in the transport and homeostasis of calcium and inorganic carbon (Benzerara et al., 2022). Moreover, in some genomes, ccyA was even colocalized with such genes.
Overall, a first step to gaining a better understanding about the molecular mechanisms involved in iACC biomineralization is to look at the expression of the ccyA gene together with other genes, including those neighboring ccyA in the genome and those mentioned above. The diel (day/night) cycle is an important parameter driving large variations of gene expression in cyanobacteria (Stöckel et al., 2008; Welkie et al., 2019; Zinser et al., 2009). For example, the diel cycle influences the expression of numerous genes involved in key cellular processes in Microcystis aeruginosa such as photosynthesis, carbon fixation, nitrogen metabolism, and circadian rhythms as shown, for example, by Huang et al. (2014) in M. aeruginosa PCC 7820. Similarly, the genes involved in a circadian cycle were studied in M. aeruginosa PCC 7806 (Straub et al., 2011) under a 24-hour light/dark cycle. They found that while the kaiA gene showed no significant variations during the cycle, the transcription patterns of the kaiB and kaiC genes, as well as the sasA gene encoding the two-component sensor histidine kinase, a KaiC-interacting protein, exhibited significant changes. These findings suggested that light is not the sole factor triggering the transcription of genes involved in photosynthesis and respiration; instead, their transcription may also be regulated by an endogenous circadian clock. Therefore, we have investigated the expression of the whole genome of Microcystis aeruginosa PCC 7806, which contains the ccyA gene and forms iACC under laboratory culture conditions, over a 24 h day/night cycle. In particular, we assessed whether the expression of genes either neighboring ccyA in the genome and/or involved in CCM and calcium transport might be correlated with that of ccyA.
Methods
Cultures and sampling
The study was conducted using the cyanobacterial strain Microcystis aeruginosa PCC 7806 available at the Pasteur culture of cyanobacteria (PCC) collection at the Institut Pasteur. The culture, sampling, nucleic acid extraction, and RNA sequencing were carried out at the the Institut Pasteur in 2017. The photosynthetic photon flux density provided during the day (50 μmol photons.m-2.s-1; cool white OSRAM L 18/640) was measured using a LICOR LI-185B quantum/radiometer/photometer equipped with a LICOR LI-193SB spherical sensor. A preculture of M. aeruginosa strain PCC 7806 was grown at 22 °C in 40 mL of the BG11 growth medium (Rippka et al., 1979) amended with 10 mM NaHCO3 under a 13/11 h day/night cycle. A 13-hour light period corresponds to what is observed in September at midlatitudes in the Northern hemisphere, a month during which Microcystis proliferates in numerous ecosystems.Then, the culture was transferred into three separate flasks (named A, B and C). The triplicate cultures of 40 mL were progressively transferred to a final set of three cultures with a volume of 1250 mL each, in fresh BG11 medium supplemented with 20 mM NaHCO3 and synchronized in an INFORS incubator with 1% CO2 and agitation at 50-100 rpm under the same day/night cycle and temperature conditions. Synchronization involved three consecutive transfers, aiming for an initial optical density (OD) at 750 nm of 0.05 and growing for a month. Sampling was performed during the mid-exponential phase (OD of 0.6-0.7) for each biological replicate, over a period of 24 h at eight time points: (i) 7 am, 60 min before the transition from night to day (t1_N); (ii) 9 am, 60 min after the transition from night to day (t2_D); (iii) 12 am (t3-D); (iv) 4 pm (t4_D); (v) 8 pm, 1 h before the transition from day to night (t5_D); (vi) 10 pm, 60 min after the transition from day to night (t6_N); (vii) 2 am, i.e. the middle of the night (t7_N), and (viii) in the end of the night period (t8_N), 24 h after the first sampling point. The sampling procedure involved the removal of 150 mL of culture, which were centrifuged, rinsed with sterile water, and centrifuged again. The resulting pellets were flash-frozen in liquid nitrogen and stored at -80 °C until further processing. To minimize batch effects, the samples were randomized before subsequent processing. We note that relatively high phosphate concentrations such as those found in BG11 (~180 μmol.L-1) can inhibit the precipitation of extracellular carbonates and/or induce that of Ca-phosphates as shown by, e.g., Rivadeneyra et al. (2006; 2010). However, by contrast, it has been shown that this does not prevent at all the precipitation of intracellular carbonates, which has been experimentally studied in BG11 before (Cam et al., 2018; De Wever et al., 2019).
Nucleic acid extraction
RNA extraction was carried out using a Trizol reagent (Life Technologies) and the NucleoSpin miRNA kit (Macherey-Nagel). Firstly, 1 mL of Trizol preheated at 65 °C was added to the frozen pellets and incubated at 65 °C for 12 min with vortexing three times. After centrifugation (10 min at 4 °C and 12,000 rcf), the upper phase was transferred to a new tube and mixed with 500 µL of Trizol preheated at 65 °C. RNA extraction was achieved by adding 400 µL of chloroform, and incubating samples for 3 min at room temperature before centrifugation (10 min at 4°C and 12,000 rcf). The upper phase was again transferred to a new tube for extraction with 400 µL of Trizol and 400 µL of chloroform with similar incubation and centrifugation conditions. Then, the upper phase obtained from this step was further processed using the NucleoSpin miRNA kit according to the manufacturer’s instructions. The final purified RNA was eluted in 20 µL of RNase/DNase-free sterile water. The absence of genomic DNA contamination was confirmed by PCR targeting the 16S rRNA gene using specific primers. The quality and quantity of RNA were assessed using a Nanodrop spectrophotometer (ThermoFisher Scientific) and a Bioanalyzer (Agilent) with RNA Nano chips.
Library Preparation and RNA Sequencing
Library preparation was conducted using the TruSeq Stranded mRNA sample preparation kit (Illumina, San Diego, California) following the manufacturer’s instructions. Ribosomal RNA (rRNA) depletion was performed using the Ribo-Zero rRNA removal kit (bacteria, #MRZB12424, Illumina) with 5 µg of total RNA. The rRNA-depleted RNA was then fragmented using divalent ions at 94°C for 8 min. The resulting fragmented RNA samples were reverse-transcribed using random primers, followed by complementary-strand synthesis to generate double-stranded cDNA fragments. There was no need for an end repair step in this process. A 3’-end adenine was added, and specific Illumina adapters were ligated to the cDNA fragments. The ligated products were amplified by PCR. The quality of the resulting oriented libraries was assessed using Bioanalyzer DNA1000 Chips (Agilent, # 5067-1504) and quantified using spectrofluorimetry (Quant-iT™ High-Sensitivity DNA Assay Kit, #Q33120, Invitrogen). Sequencing was performed on the Illumina HiSeq2500 platform of the Institut Pasteur, generating single-end 65 bp reads with strand specificity. Each lane of sequencing accommodated a mixture of 18 multiplexed samples.
Read processing
Raw RNA-seq reads (available online, see section Data, scripts, code, and supplementary information availability) were cleaned from adapter sequences and low-quality sequences using an in-house program (https://github.com/baj12/clean_ngs). Only sequences of at least 25 nt in length were considered for further analysis. Transcript abundance was quantified from the cleaned read dataset using the Salmon software (Patro et al., 2017). Pseudo-mapping of reads was done considering the coding sequences (CDS, including pseudogenes) of the complete RefSeq genome assembly of M. aeruginosa PCC 7806 as a reference (GCF_002095975.1). For each replicate, the estimated transcript abundances were normalized by Salmon according to the transcript size, genome size (number of CDS) and sample size (number of reads).
Differential expression analysis
Total abundance data were further processed using the DiCoExpress R-script-based tool (Lambert et al., 2020) in order to detect genes expressed differentially between different time points. First, normalization was achieved using Trimmed Mean of the M-value (Robinson and Oshlack, 2010). During this procedure, 220 low-count transcripts were discarded, resulting in a total of 4,440 genes with normalized counts. Principal component analysis (PCA) was performed on this normalized dataset using the R packages FactoMineR and Factoshiny (Lê et al., 2008). Then, the differential expression between time points was tested using a generalized linear model. For comparison between each time point, statistics of differential gene expression are provided as Fold Change (FC) and p-values.
Gene functional annotation
In addition to the annotations provided by the NCBI prokaryotic genome annotation pipeline for the genome assembly of M. aeruginosa PCC 7806 (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_002095975.1), protein sequences translated from CDS were processed using DeepNOG protein orthologous groups assignment (Feldbauer et al., 2021) as implemented in the Genovi tool (Cumsille et al., 2023). An orthologous group with a COG accession number is assigned to each gene, as well as a confidence score ranging within [0, 1]. The stringent threshold was set at 0.8. Orthologs of genes involved in carbon-concentrating mechanisms (CCM) and calcium transport were specifically detected as bidirectional best BLASTp hits using protein sequence queries described in (Tang et al., 2022) (15 genes from Synechocystis sp. PCC 6803, and the ecaA gene from Anabaena PCC 7120) and (De Wever et al., 2019) (13 genes from Chroococcidiopsis thermalis PCC 7203), respectively. The following thresholds were used: E-value cut-off of 1E-6, ≥30% identity and 70% coverage.
Sequence and structure analysis of unknown proteins
General information on the three domains of unknown function (DUF) encoded by genes located in the genomic neighborhood of the ccyA gene was extracted from their respective InterPro entries (Paysan-Lafosse et al., 2023; DUF5132: IPR033456; DUF454: IPR00740; DUF1269: IPR009200). AlphaFold2 (AF2) 3D structure models of the three DUF-containing proteins and two hypothetical proteins neighboring ccyA were build using ColabFold (Mirdita et al., 2022). Foldseek (van Kempen et al., 2024) was used to search for structural similarities against the Protein Data Bank (PDB) and AlphaFold Database (AFDB). 3D structures were manipulated using Chimera (Pettersen et al., 2004). Hydrophobic cluster analysis (HCA) (Callebaut et al., 1997) was used to analyze protein secondary structure features.
Results
Overall expression pattern
A principal component analysis (PCA) was performed on the transcript abundances of the 4,440 genes, for each replicate at each sampling time (Figure 1). All the data points were projected onto a two-dimensional plane defined by the first two principal components, which accounted for 64.69% of the total variance (Figure 1). Among the triplicates (i.e. three independent cultures for each time-point) only minor variations were observed, at least along axis 1, and their values were clearly separated from the samplings at other time points. Furthermore, the resulting plot exhibited a clear separation along axis 1 between night (N) and day (D) points.
Expression of the ccyA gene
In order to analyze the expression profile of the ccyA gene during the 24 h day/night cycle, we plotted the normalized expression of ccyA in all replicates at the different time points (Figure 2). A clear day/night pattern can be observed in the expression of ccyA, with maximum transcript abundances during the second half of the night (t7_N and t8_N, Table 1) and a minimum transcript abundance at 12 am during daytime (t3_D, Table 1). Pairwise comparisons of ccyA transcript abundance at different time points show that the expression of ccyA at daytime points (t2_D to t5_D) did not significantly differ between each other, but were significantly different (p-value < 0.01) from some of the nightime points, t7_N and t8_N, respectively (Table S1). The 10 pm nighttime point (t6_N) did not significantly differ from any others except from the 12 am daytime point (t3_D; p-value < 0.05) and can therefore be considered as a transition point between day and night. When compared with the mean expression of the 4,440 genes of the whole transcriptome, the mean expression of the ccyA gene was significantly higher (p-value < 0.05) during the second half of the night, i.e. at the two-time points 2 am (t7_N) and 7 am (t8_N) (Table 1).
Figure 1 - Principal component analysis of the transcript abundances of 4440 M. aeruginosa genes during a day/night cycle. Abundances are expressed as log base 2 of the normalized counts (see Methods section for normalization procedure), for each replicate (a-c) at each sampling time (t1-t8). The 24 samples are projected onto the first two principal components (Dim. 1 and Dim. 2), explaining 48.9% and 15.8% of the variance of the dataset, respectively. Samples/replicates collected at different times over the day/night cycle appear with a different color code. The “_N” and “_D” codes in the sample names refer to night and day, respectively.
Figure 2 - Time course of the abundance of ccyA transcripts during a day/night cycle. Abundances are expressed as log base 2 of the normalized counts. The abundance mean for the three replicates is shown by a black filled circle. Replicates a, b and c are represented by different symbols. For each sampling time (t1-t8), error bar represents the standard error on the mean of the three replicates. Grey shaded areas outline the night periods.
Table 1 - Comparison between the transcript abundances of ccyA and the whole transcriptome, for each time point. Abundances are expressed as log base 2 of the normalized counts. When means significantly differ, p-values of Wilcoxon signed rank test with continuity correction are depicted as follows: (**): 0.01 – 0.05, (*): 0.05 – 0.1. Night points are presented in bold.
t1_N(*) 7 am | t2_D 9 am | t3_D 12 am | t4_D 4 pm | t5_D(*) 8 pm | t6_N(*) 10 pm | t7_N(**) 2 am | t8_N(**) 7 am | |
|---|---|---|---|---|---|---|---|---|
Mean (ccyA gene, n=3 replicates) | 10.31 | 9.07 | 8.88 | 9.12 | 9.31 | 9.73 | 10.39 | 10.43 |
Standard deviation (ccyA gene, n=3 replicates) | 0.15 | 0.22 | 0.23 | 0.19 | 0.14 | 0.34 | 0.33 | 0.12 |
Mean (4440 genes, n=3 replicates) | 6.87 | 7.35 | 7.39 | 7.44 | 7.44 | 7.36 | 7.23 | 7.18 |
Standard deviation (4440 genes, n=3 replicates) | 3.79 | 2.45 | 2.57 | 2.67 | 2.61 | 2.87 | 3.21 | 3.33 |
Search of differential gene expression during the day/night cycle
Subsequently, a comprehensive analysis was conducted throughout the genome to identify all genes exhibiting diel variations in their expression, using the minima and maxima of ccyA gene expression (t3_D and t8_N) as reference time points (Figure 3). Overall, 1,799 of the 4,440 total genes showed a significant difference in expression (at least 2-fold) between 12 am (during the day, t3_D) and 7 am (during the night, t8_N) (Figure 3 and Table S2). Among these 1,799 genes, 906 were significantly more expressed during the day (t3_D), and conversely, 893 were over-expressed during the night (t8_N).
Figure 3 provides the annotations for the genes with the most substantial expression differences between these two time points (-log10(p-value)>60): genes overexpressed during the day included a NAD(P)H-quinone oxidoreductase and a protein of the UDP-glucose/GDP-mannose dehydrogenase family. The genes overexpressed at night included: two genes encoding hypothetical proteins, a pseudo-gene encoding a protein containing a helix-turn-helix domain, a gene encoding a protein from the SDR oxidoreductase family, and a gene encoding a dihydroorotate dehydrogenase-like protein.
Genomic neighborhood of ccyA
Figure 4 displays the map of the 20 genes (18 genes and 2 pseudo-genes) located upstream (10 genes) and downstream (10 genes) to ccyA. This number of 20 genes was set arbitrarily and considered as a close neighborhood. When available, the functional annotations were extracted from the NCBI RefSeq database (Table 2). In this database, the ccyA gene is incorrectly annotated as a gene coding for the cell envelope biogenesis protein OmpA (Benzerara et al., 2022). Interestingly, three transposition elements were found near the ccyA gene: an ISL3 family and an IS1 family transposases, both pseudo-genes located at -6 and -5 genes compared to ccyA, and an IS1634 family transposase located at +2 genes compared with ccyA. The gene located directly upstream of ccyA was annotated as encoding a domain-of-unknown-function (DUF5132)-containing protein, and the one directly downstream a calcium/proton exchanger (cax gene), previously characterized in eukaryotes (Hirschi et al., 1996; Waight et al., 2013). At position -3 (upstream of ccyA), another DUF-containing protein (DUF454) was detected.
Figure 3 - Differential expression of the whole transcriptome between sampling time t3_D (12 am, during the day period) and t8_N (7 am, end of the night period). Log base 2 of the fold change (FC) values (x-axis) are defined for 3116 CDS. CDS with a negative log2(FC) value correspond to transcripts significantly more abundant at 12 am (t3_D) compared with 7 am (t8_N) (gold area; 906 CDS). CDS with a positive log2(FC) value correspond to transcripts significantly more abundant at 7 am (t8_N) compared with 12 am (t3_D) (grey area; 893 CDS). The ccyA gene is indicated by a red dot. The horizontal red dashed line indicates the significance p-value threshold after Bonferroni correction (-log10(P-value)=5). The vertical red dashed lines indicate FC=2 and FC=0.5, e.g., 2-fold abundance ratio between the two sampling times. RefSeq accessions and annotations are indicated for CDS with P<2.0e-66). See Suppl. Table S2 for details.
Figure 4 - Schematic representation of the ccyA gene neighborhood in M. aeruginosa PCC 7806 complete genome assembly (NCBI genome browser display). The ccyA gene is outlined by a green box. RefSeq protein accession, gene name (if defined), gene orientation (red for minus; green for plus) and locus tag are indicated for the 10 CDS located upstream (-10 to -1) and the 10 CDS located downstream (+1 to +10) with respect to ccyA. Protein functional annotation is provided in Table 2. This genomic region spans within the nucleotide position interval [4515650, 4535259] in the NCBI genome assembly NZ_CP020771.1. CDS boundaries and features correspond to RefSeq annotations (GCF_002095975.1). Gene symbols stand for shc: squalene—hopene cyclase; tsf: translation elongation factor Ts; rspsB: 30S ribosomal protein S2 ; cax: calcium/proton exchanger; acs: acetate—CoA ligase; psaC: photosystem I iron-sulfur center protein PsaC.
Table 2 - Neighborhood of the ccyA gene in M. aeruginosa PCC 7806 complete genome assembly. General features and annotations are listed for the ten CDS located upstream (-10 to -1) and the ten CDS located downstream (+1 to +10) of the position of ccyA (in bold, protein improperly annotated in RefSeq as OmpA). RefSeq annotations correspond to the assembly accession GCF_002095975.1 (NZ_CP020771.1).
Protein accession (NCBI RefSeq) | CDS position | Gene orientation | Start (nt) | End (nt) | Gene name | Protein functional annotation (NCBI RefSeq) | COG (DeepNOG annotation) | |||
|---|---|---|---|---|---|---|---|---|---|---|
NCBI RefSeq | relative to ccyA | Accession | Confidence | Category | ||||||
WP_002741645.1 | 4262 | -10 | + | 4515650 | 4516447 | rpsB | 30S ribosomal protein S2 | COG0052 | 1.0 | J |
WP_002741643.1 | 4263 | -9 | + | 4516564 | 4517319 | tsf | translation elongation factor Ts | COG0264 | 1.0 | J |
WP_002741642.1 | 4264 | -8 | + | 4517577 | 4519484 | shc | squalene--hopene cyclase | COG1657 | 1.0 | I |
WP_051048450.1 | 4265 | -7 | - | 4519545 | 4519886 | transcriptional regulator | COG5499 | 0.99 | V | |
pseudo (partial=3') | 4266 | -6 | + | 4520389 | >4520916 | ISL3 family transposase | ||||
pseudo (partial=3') | 4267 | -5 | - | <4520947 | 4521127 | IS1 family transposase | ||||
WP_002753396.1 | 4268 | -4 | - | 4521343 | 4521633 | DUF5132 domain-containing protein | COG4980 | 0.99 | R | |
WP_036399769.1 | 4269 | -3 | + | 4521771 | 4522358 | DUF454 family protein | COG2205 | 0.28 | T | |
WP_228036328.1 | 4270 | -2 | + | 4522489 | 4522689 | hypothetical protein | COG2161 | 0.20 | V | |
WP_002741624.1 | 4271 | -1 | + | 4522868 | 4523182 | DUF5132 domain-containing protein | COG1826 | 0.96 | U | |
WP_084990073.1 | 4272 | 0 | + | 4523229 | 4524278 | ccyA | cell envelope biogenesis protein OmpA | COG3134 | 0.95 | S |
WP_084990074.1 | 4273 | +1 | + | 4524268 | 4525377 | cax | calcium/proton exchanger | COG0387 | 1.0 | P |
WP_084990075.1 | 4274 | +2 | - | 4525733 | 4527469 | IS1634 family transposase | COG5421 | 0.99 | X | |
WP_002734319.1 | 4275 | +3 | - | 4527740 | 4527985 | psaC | photosystem I iron-sulfur center protein PsaC | COG5703 | 1.0 | C |
WP_002741619.1 | 4276 | +4 | + | 4528237 | 4530207 | acs | acetate--CoA ligase | COG0365 | 1.0 | I |
WP_084990076.1 | 4277 | +5 | + | 4530562 | 4531158 | hypothetical protein | COG1196 | 0.40 | D | |
WP_002741615.1 | 4278 | +6 | + | 4531198 | 4531722 | GNAT family N-acetyltransferase | COG1670 | 0.99 | J,O | |
WP_002731545.1 | 4279 | +7 | - | 4531850 | 4532356 | DUF1269 domain-containing protein | COG4803 | 1.0 | S | |
WP_084990077.1 | 4280 | +8 | - | 4532613 | 4533059 | DUF1269 domain-containing protein | COG4803 | 0.99 | S | |
WP_002756663.1 | 4281 | +9 | - | 4533361 | 4533978 | CBS domain-containing protein | COG5792 | 1.0 | C,T | |
WP_002741611.1 | 4282 | +10 | - | 4534069 | 4535259 | glycosyltransferase family 4 protein | COG0438 | 1.0 | M | |
AlphaFold 2 (AF2) 3D structure models of two hypothetical proteins, as well as five proteins containing DUFs present in the ccyA genomic neighborhood were analyzed (Figure 5) and compared with experimental 3D structures using Foldseek, in order to highlight possible structural and functional relationships. (i) Analysis of the DUF454 family protein (WP_036399769.1; position -3 relatively to ccyA) indicates that its C-terminal end matches the C-terminal end of the DUF454 IPR profile (inner membrane protein YbaN from E.coli). Meanwhile, its N-terminal part matches the Pfam HMA_2 profile (PF19991), corresponding to a CoBaHMA domain. CoBaHMA (after Conserved Basic residues HMA) is the name of a recently identified family of domains belonging to the HMA superfamily (Gaschignard et al., 2024). Within this superfamily, the CoBaHMA domain family shows distinct features such as the presence of an additional strand β0 at the N-terminus of the domain and the presence of a charged patch on one side of the β-sheet. The AF2 3D structure model also consistently highlights this two-domain architecture, with the C-terminal, DUF454-like domain folding as a helical hairpin, while the CoBaHMA domain possesses the basic conserved amino acids characteristic of this family (Figure 5-A). (ii) The DUF5132-containing protein (WP_002741624.1; position -1 relatively to ccyA) is covered by the DUF5132 profile (aa 48-89) over a little less than half of its length and contains a predicted transmembrane segment in its N-terminal part (aa 12-32). The AF2 3D structure model, characterized overall by good pLDDT values, indicates a putative assembly, in a hairpin-like conformation, of the helical transmembrane segment and another alpha-helix preceding a long, soluble helix matching the DUF5132 profile (Figure 5-B). (iii) A second member of this family, with similar structural features is also present in the ccyA neighborhood (WP_002753396.1; position -4 relatively to ccyA). No significant relationship with any known 3D structures could be detected by FoldSeek. (iv) A hypothetical protein (WP_228036328.1; position -2 relatively to ccyA) corresponds to a truncated part of a hypothetical protein referenced under RefSeq WP_271988195.1, which is also composed of a CoBaHMA domain sharing the basic signature of this family (Figure 5-C). (v) Another hypothetical protein (WP_084990076.1; position +5 relatively to ccyA) is predicted as a long coiled-coil (data not shown). (vi) and (vii) Finally, two DUF1269-containing proteins (WP_002731545.1 and WP_084990077.1; genes at +7 and +8 relatively to ccyA), corresponding to two isoforms differing by their C-terminal end, match the DUF1269 IPR profile over their whole lengths and include a long glycine-zipper motif embedded in a soluble domain, which based on Foldseek searches, corresponds to a ferredoxin fold, as found for instance in the copper tolerance CutA1 protein (pdb 1V6H, Bagautdinov, 2014) and in the nickel-responsive transcription factor NikR (pdb 2BJ3, Chivers and Tahirov 2005) (Figure 5-D). This glycine-zipper motif is similar, albeit more hydrophobic, to the GlyZip motif in the C-ter domain of calcyanins, composed of two parts separated by a proline preceded by a small amino acid (here a serine, while it is a glycine in calcyanins) (Figure 5-D). Overall, this analysis highlighted some features in the neighborhood of the ccya gene, which were not captured by automatic annotations, including motifs that are found in calcyanins of all cyanobacteria (GlyZip and hairpin-like helices) as well as calcyanins of some species (CoBaHMA domains).
Correlated expression of other genes with the ccyA gene
The expression patterns of genes exhibiting a significant difference in their circadian expression between daytime (t3_D, 12 am) and nighttime (t8_N, 7 am) were compared with that of ccyA to figure out which ones may have a maximum expression at the same time (Figure 6; Table S3). We note that such gene expression variations over a day/night cycle are expected in an autotrophic organisms and at least some correlations are likely not causal. However, their identification offers one possibility to better understand which other pathways may be activated at the same time as ccyA is expressed. As a result, two subsets were identified: on the one hand, 773 genes showed a highly positive correlation with the expression pattern of ccyA (with a Pearson correlation coefficient rho > 0.75, this threshold being arbitrarily fixed; grey subset in Figure 6; see also Figure S2 and Table S3); on the other hand, 706 displayed a negative correlation with the expression pattern of ccyA (rho < -0.75, yellow subset in Figure 6; see also Figure S2 and Table S3). Figure 6 displays the significance level (p-value) of the correlations between the expression profiles of these 1479 genes and that of ccyA, and provides annotation for the genes with the most significant correlation.
Figure 5 - AlphaFold2 3D structure models of some proteins encoded by the genomic neighbors of ccya. The 3D structure models are represented as ribbons colored according to the plDDT values (color code given at bottom). In panels A and C, the basic amino acids constituting the signature of the CoBaHMA family are shown as balls and sticks. (A) The CoBaHMA_DUF454-Cter architecture (encoding gene in position -3 relatively to ccyA) is compared with a protein entirely covered by the DUF454 IPR profile (inner membrane protein YbaN (YABN_ECOLI, UniProt P0AAR5); insert). (B) DUF5132 protein (encoding gene in position -1 relatively to ccyA). TM stands for TransMembrane. (C) Hypothetical protein (encoding gene in position -2 relatively to ccyA) with a typical CoBaHMA domain. (D) DUF1269 domain-containing protein (encoding gene in position +7 relatively to ccyA) containing a GlyZip motif, predicted to fold as a helical hairpin (see HCA plot and sequence alignment in Figure S1).
The genes for which the expression patterns were the most significantly positively correlated with that of the ccyA gene were: (i) an intermembrane mitochondrial space (IMS) domain-containing protein, (ii) a GTP-binding protein, and (iii) a DUF-containing protein (DUF4278). The genes the most significantly anticorrelated with ccyA expression profile were annotated as (i) an ABC transporter ATP-binding protein/permease, (ii) an adenylyl-sulfate kinase, and (iii) a hypothetical protein.
Among the 20 neighboring genes of ccyA, eight showed an expression profile either positively or negatively correlated with that of ccyA (see Figure S3 for a time plot of expression of all 20 genes). Interestingly, the expression profile of the two genes adjacent to ccyA, located on the same DNA strand, i.e. encoding for a calcium/proton exchanger (CAX) and DUF5132 (Table 2), displayed a strong positive correlation with the expression profile of ccyA (Figure 6 and 7). The expression profile of the gene coding for a DUF454 family protein (position -3 relatively to ccyA; Table 2), exhibited one of the strongest positive correlations with the expression profile of ccyA (Figure 6), with a nearly identical expression profile (Figure 7). Finally, the expression profiles of two other genes, located at +10 and +5 genes with respect to ccyA showed a weaker but significant positive correlation with the expression profile of ccyA. The gene at +10 was annotated as a glycosyltransferase family 4 protein. The gene at +5 codes for a hypothetical protein predicted as a long coiled-coil (Table 2).
Moreover, proteins involved in calcium transport and in CCM were searched in the genome of PCC 7806 (Table S4). In total, 24 genes involved in CCM and/or their homologs and 24 genes involved in calcium passive or active transport and their homologs were identified. The expression patterns of seven out of the 24 genes involved in CCM were correlated negatively with the expression pattern of ccyA, whereas the expression patterns of five of CCM genes were positively correlated with the expression pattern of ccyA (Figure 6, in magenta; see Figure S4 for a time plot of the expression of all CCM genes). Among the 24 genes involved in calcium transport, the expression patterns of five of them were negatively correlated with that of ccyA, while the expression patterns of ten of them were positively correlated with that of ccyA (Figure 6, in cyan except the cax gene in green; see Figure S5 for a plot of time expression of all “Ca genes”).
Figure 6 - CDS with a day/night expression pattern strongly (anti-)correlated with that of ccyA. The Pearson correlation coefficient (rho) is highly negative for 706 CDS (left). The Pearson correlation coefficient is highly positive for 773 CDS (right). The distribution of correlation coefficients over the 4440 CDS is shown in Figure S2. RefSeq accessions and annotations are indicated for the 3 CDS with the strongest negative and positive correlations (see Table S3 for details). CDS belonging to the carbon concentrating mechanism or the calcium transport system are represented as magenta and cyan dots, respectively. CDS located in the [-10, +10] neighborhood of ccyA (Table 2 and Figure 4) are represented by green and red dots (for positive and negative DNA strands, respectively). The cax gene located at a genomic position of +1 from the ccyA gene also belongs to the calcium transport system and appears in green (see Table S4 for details).
Discussion
The identification of a marker gene of iACC biomineralisation in cyanobacteria was a first step towards understanding the molecular mechanisms controlling the formation of these biomineral phases (Benzerara et al., 2022). However, up to now, only the presence/absence of the gene and the structural features of the encoded protein have been investigated (Benzerara et al., 2022; Gaëtan et al., 2023; Gaschignard et al., 2024). Here, we explored whether this marker gene is transcribed in an organism capable of forming iACC, how this may be influenced by day/night cycles and how the expression of neighboring genes compares.
Several tools can be used for assessing the transcription of a large number of genes within a genome, including DNA microarray and high-throughput RNA sequencing (RNAseq) (Wang et al., 2009). Microarrays are limited by the use of probes: the number of studied genes remains limited by the available knowledge of the targeted genome. By contrast, RNAseq allows targeting all the mRNAs in a cell and has the advantage of measuring the level of transcriptional activity in the organisms (Mantione et al., 2014). Moreover, RNAseq can achieve higher resolution (i.e. higher accuracy in the estimation of gene expression) than microarrays. Because this technique does not suffer from microarray-based limitations such as background noise and saturation, it can also achieve a lower detection limit, and thus detect very low gene expression (Zhao et al., 2014). Straub et al. (2011) already studied the expression of the whole genome of M. aeruginosa PCC 7806 over a 24-h day/night cycle, using DNA microarray. Using their publicly available dataset, we determined that the expression level of the ccyA gene (MIC_2500) was not significantly different from the background noise. Here, the use of RNA-seq on the same organism enabled the detection and analysis of the ccyA gene expression dynamics. Indeed, the expression level of ccyA in our dataset was above the mean expression of the 4,440 genes of Microcystis at all time and this was statistically significant for several time points during the night. We can thus conclude that the ccyA gene is clearly expressed in an iACC-forming strain.
Figure 7 – Time course of the abundances of three genomic neighbors of the ccyA gene during a day/night cycle. The CDS located at positions -3 (blue), -1 (purple) and +1 (orange) correspond to the ccyA neighbors with transcript abundance profiles the most correlated to ccyA transcript abundance profile: Pearson correlation coefficient of 0.97 (p-value=1.7e-15), 0.95 (p-value=8.0e-13), 0.94 (p-value=1.2e-11), respectively. See Methods section for more details about how counts were normalized.
In line with the study by Straub et al. (2011) on the dynamics of gene expression in M. aeruginosa PCC 7806, we found, as expected, that the day/night cycle is an important factor driving variations in global gene expression in PCC 7806. Cyanobacteria have a highly defined circadian clock, and this day/night gene expression pattern has been demonstrated in several strains (Pattanayak & Rust, 2014; Stöckel et al., 2008; Zinser et al., 2009). Straub et al. (2011) showed that the metabolism of M. aeruginosa undergoes significant changes between the light and dark periods. Specifically, during the light period, processes including carbon uptake, photosynthesis, and the reductive pentose phosphate pathway result in glycogen synthesis. Conversely, during the dark period, glycogen degradation, the oxidative pentose phosphate pathway, the tricarboxylic acid (TCA) branched pathway, and ammonium uptake promote amino acid biosynthesis. Furthermore, the biosynthesis of secondary metabolites, such as microcystins, aeruginosin, and cyanopeptolin, predominantly occurs during the light period.
A clear day-night expression pattern was also observed for the ccyA gene: the expression was significantly higher than the average genome expression from the middle (2 am, t7_N) to the end of the night (7 am, t8_N). In all other time points of the 24-h cycle, the expression level of this gene did not statistically differ from the average expression of all genes. The gene is therefore overtranscribed during the second half of the night. While mRNA abundance provides valuable information about gene expression, it does not always directly correlate with the protein levels in the cell (Ingolia et al., 2009; Riba et al., 2019; Schwanhäusser et al., 2011). Several factors, including translation, initiation and elongation rate (Riba et al., 2019), post-translational modifications, protein stability, and degradation rates (Christiano et al., 2014; Gancedo et al., 1982) may impact the relationship between mRNA abundance and protein content. As such, it is difficult to infer the abundance of calcyanin at a given time from the transcript expression level only. However, one can speculate that the protein may be more produced during the second part of the night, and remain at an optimal abundance during the day when it may most function. In the future, direct measurements of the protein abundance should be performed using proteomics in order to quantify the calcyanin content in cells. Last, Straub et al. (2011) mentioned that the metabolism of M. aeruginosa strongly changes between the light and the dark period. For example, during the latter, amino acid biosynthesis is promoted by glycogen degradation, the oxidative pentose phosphate pathway, the TCA branched pathway and ammonium uptake. This illustrates that seeing ccyA expression at its strongest at the end of the dark period is not particularly surprising energetically-wise.
Interestingly, we showed that a Ca2+/H+ antiporter (cax gene) and two DUF-containing proteins were not only directly co-localized upstream and downstream of ccyA, on the same DNA brand but they were also co-expressed with ccyA. This suggests that those three genes could be part of a functioning unit of DNA that could be transcribed together into a single mRNA strand, i.e. an operon, and may play a concerted role in iACC formation.
The process of iACC formation necessitates the accumulation of calcium ions at concentrations higher than the cytosolic concentrations typically measured in cells (Cam et al., 2015). Therefore, it may be suggested that the protein encoded by the cax gene is involved in the build-up of these accumulations. However, the functioning of CAX proteins raises questions. Indeed, the single Ca2+/H+ exchanger found in Synechocystis PCC 6803 has been shown to be located in the plasma membrane and to catalyze Ca2+ efflux out of the cells (Waditee et al., 2004). By contrast, it has been suggested that CAX proteins in plants may have different locations and efflux Ca2+ from the cytoplasm either to the extracellular space or to intracellular vacuoles (Shigaki et al., 2006). Investigating further the similarities/dissimilarities of the CAX protein of M. aeruginosa PCC 7806 with already known CAX proteins and locating these transporters within M. aeruginosa PCC 7806 cells would therefore be crucial to better understand their potential involvement in iACC formation.
The gene encoding a DUF454 protein that co-localizes (position -3 relatively to ccyA) and co-expresses with ccya is particularly interesting. Indeed, we show here that this DUF454 protein bears a CoBaHMA (after Conserved Basic residues HMA) domain (Gaschignard et al., 2024). Strikingly, the calcyanins of many cyanobacteria are similarly composed of a CoBaHMA domain in their N-terminal region and a C-terminal three-fold repeat of a long glycine zipper motif conserved in all iACC-forming cyanobacteria (Benzerara et al 2022). This is not the case for the calcyanins of Microcystis, which does not contain a CoBaHMA domain but instead displays an unknown domain in their N-terminal region. Since the discovery of the calcyanin family, it has been surprising to observe such a broad diversity in their N-terminal domains despite an expected common function. Here, although we still ignore the function of calcyanin, the discovery of a CoBaHMA-encoding gene neighboring the ccyA gene suggests that Microcystis calcyanin may work in concert with other proteins, including a DUF454 protein and resulting in a protein complex with a possibly similar architecture in Microcystis and other iACC-forming cyanobacteria. Moreover, the presence in the close genomic neighborhood of ccyA of genes encoding for another protein of the CoBaHMA family (-2 position relative to ccyA), as well as proteins possessing a GlyZip-like motif (+7 and +8 positions relative to ccyA) relatively similar to the GlyZip motif in the C-ter domain of calcyanins additionally suggests the possibility of complementary roles of these proteins for a function to identify.
Moreover, building upon the approach followed by De Wever et al. (2019), we conducted a comprehensive search for annotated genes involved in CCM and calcium transport, as well as their homologous counterparts in the genome of M. aeruginosa PCC 7806. Firstly, a greater number of genes associated with CCM or their homologs appeared to be expressed preferentially during the day instead of night (nine vs four, respectively), compared with genes associated with calcium transport (three vs nine, respectively) (Table S4). Similarly, Straub et al. (2011) showed that the transcript abundance of CCM genes showed a peak in expression just after the night/day transition and then a rapid decrease after the day/night transition. The elevated concentration of CO2 in carboxysomes resulting from the CCM is a way to improve carbon fixation (Price and Howitt, 2011). The expression of CCM genes during the day is therefore expected. Yet, among the 24 genes identified as involved or potentially involved in CCM, 5 had an expression pattern positively correlated with that of ccyA, i.e. maximum at night (Table S4). Notably, two of these genes (cmpC and cmpD) encode two subunits of the inorganic carbon transporter BCT1 involved in the cyanobacterial CCM (Price & Howitt, 2011). BCT1 is a high-affinity but low-flux HCO3- transporter. It is composed of four subunits, and among them, CmpC and CmpD are extrinsic cytoplasmic proteins that possess binding sites for ATP (Price & Howitt, 2011). It has been shown that the expression of these genes in M. aeruginosa PCC 7806 depends on the ambient pCO2 (Sandrini et al., 2015). Notably, under high pCO2 conditions (1450 ppm in their experiments), the expression of high-affinity bicarbonate uptake systems (notably cmpA, cmpB, and cmpC genes) is down-regulated, and cells shift towards CO2 and low-affinity bicarbonate uptake compared with cells grown under low pCO2 conditions (200 ppm) that the authors considered as Ci-limited. For the production of the transcriptomes that we studied here, pCO2 was not regulated nor monitored, and the change of expression in CCM-involved genes such as cmpD and cmpC could be due to changes in pCO2 conditions during day/night transitions, even if there are expected to be lower than in the Sandrini et al. (2015)’s study. Moreover, we note that cmpA and cmpB were more expressed during the day. Therefore, the potential variation of the BCT1 activity over a diel cycle remains to be ascertained. In future studies, exploring whether the co-expression of the five genes involved in CCM with ccyA may play a role in iACC biomineralization would be interesting, and should include a special focus on the pCO2 variations during the experiments. Last, Walter et al. (2016) demonstrated that in Anabaena sp. PCC 7120, the concentration of extracellular Ca2+ can also impact the expression of genes, including those involved in bicarbonate uptake. While we do not expect a dramatic change of the extracellular Ca2+ concentration over a 24-hour period based on previous studies (Cam et al., 2018; De Wever et al., 2019), this parameter should also be considered in future studies. In contrast to the bicarbonate transporters and their homologs, the number of genes possibly linked with calcium transport (24 genes) and having an expression pattern positively correlated with that of ccyA (10 genes) was higher than genes with an expression pattern negatively correlated with ccyA (5 genes) (Table S4). The transport of Ca2+ ions into the cytosol is commonly considered as a passive process facilitated by channels that exhibit low ionic specificity (Domínguez et al., 2015). Here, only one out of the ten genes whose expression was positively correlated with that of ccyA was identified as a passive transporter (hB1-1). The other nine genes include the cax gene, neighboring ccyA, two apnhaP genes (Na+/H+ antiporter acting as a Ca2+/H+ antiporter at alkaline pH as shown by Waditee et al. (2001) and three of their homologs, a UPF0016 gene (putative Ca2+/H+ transporter) and a Ca2+-ATPase and one of their homologs. They could thus also be candidate genes contributing to the intracellular accumulation of calcium with the same question arising about the polarity of the transport as for cax (towards the exterior of the cells or inside compartments containing iACC).
In a comparative genomics approach performed at the scale of the whole Cyanobacteriota phylum, ccyA was the only gene identified as shared by iACC+ strains and absent in iACC- strains (Benzerara et al., 2022). However, some genes from both phenotypes (iACC- and iACC+), and therefore not highlighted by the above-mentioned comparative genomics approach, could also be involved in iACC formation. Indeed, these genes may be recruited in iACC+ strains for iACC formation, whereas they would have a different function in iACC- strains. Moreover, there could also be genes specific to Microcystis involved in iACC formation, with other non-homologous genes but of similar functions involved in other cyanobacterial genera. The present study is a step forward towards identifying possible partners in iACC formation. However, only future studies specifically targeting these genes will be able to provide definitive answers about their implications in the iACC formation process. Furthermore, investigations into the comprehensive response of Microcystis to fluctuations in calcium and/or pCO2 levels could provide valuable insights into the molecular mechanisms underlying the biomineralization of iACC. These studies should encompass monitoring of the physicochemical parameters of the environment and the quantification of iACC formation. Notably, certain semi-quantitative techniques, such as Fourier transform infrared spectroscopy (FTIR) (Mehta et al., 2022b), scanning transmission x-ray microscopy (STXM) (Benzerara et al., 2023) and X-ray absorption near edge structure (XANES) spectroscopy (Mehta et al., 2023), now offer the capability to detect cyanobacterial iACC and assess their relative abundance. In the future, studying control strains not hosting the ccyA gene might be an additional interesting perspective. However, for this purpose, one cannot consider any ccyA- strain as a robust control since its genome may differ from the one of PCC 7806 in several genes. A more reliable control would consist in a yet-unavailable mutant of PCC 7806 with a desactivated/deleted ccyA gene. Last, one may speculate that cells may sometimes, also direct iACC dissolution to regulate their abundance. While this process has been shown in Achromatium (e.g., Yang et al., 2019), it still needs to be unambiguously evidenced in Cyanobacteriota before assessing which genes might be involved.
Finally, it can be noted that the correlation between the ccyA gene and genes involved in calcium transport and the CCM is not restricted to Microcystis. For example, the combined presence of homologs of a Ca2+/H+ antiporter gene and the Na+-dependent bicarbonate bicA gene was shown to be significantly correlated with that of ccyA in the Cyanobacteriota phylum (Benzerara et al., 2022). Moreover, the ccyA gene was co-localized with a Ca2+/H+ antiporter gene in several other cyanobacteria (C. fritschii PCC 9212 and PCC 6912, and Fischerella sp. NIES-4106), mirroring the situation observed in M. aeruginosa PCC 7806. Overall, this offers some target genes for future genetic studies to better depict the molecular mechanisms of iACC biomineralization.
Acknowledgements
We thank Varet H, Sismeiro O, Dillies MA, Legendre R and Coppe JY who run the RNA platform of the Institut Pasteur and provided us with the processed trancriptomics data. We also thank two reviewers for comments that significantly improved the manuscript. Preprint version 3 of this article has been peer-reviewed and recommended by Peer Community In Microbiology (https://doi.org/10.24072/pci.microbiol.100153; De Wit, 2025).
Data, scripts, code, and supplementary information availability
Raw transcriptomics data are available online on NCBI, Gene expression Omnibus, GEO accession number: GSE255450 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255450). Data (output data from Salmon, DICOEXPRESS and DeepNOG analyses as well as gene expression correlation), supplementary information, supplementary figures and tables are available online on available on Zenodo (https://doi.org/10.5281/zenodo.14739957, Duprat & Benzerara 2025). Links to the tools used in this study are noted in the Material and methods section.
Conflict of interest disclosure
The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.
Funding
This work was supported by the Agence Nationale de la Recherche (ANR Harley, ANR-19-CE44-0017-01; ANR PHOSTORE, ANR-19-CE01-0005). The PCC collection is funded by the Institut Pasteur. The PhD of C.P. was supported by the Ile-de-France ARDoC Grant. Juliette Gaëtan PhD grant was funded by the Learning Planet Institute and University Paris Cité grant Frontiers of Innovation in Research and Education.