1 INTRODUCTION

Microcystis is one of the most notorious phytoplankton species, acting as a dominant species on the lake eutrophication (Mohamed et al., 2003; Robson and Hamilton, 2004; Gan et al., 2012). Studies suggested that Microcystis could be the dominant species in the aquatic ecosystem, mainly due to its ecological advantages over other phytoplankton species (Paerl et al., 1985; Kruk et al., 2017; Esterhuizen-Londt et al, 2018). For example, Microcystis is characterized with sheaths and mucilage, which increases the floating velocity and resistance of cells (Cyr and Curtis, 1999; Kearns and Hunter, 2000; Wu and Kong, 2009; Yamamoto et al., 2011). In eutrophication water, Microcystis was able to float to water surface to capture light and CO₂, and outcompete the growth of other planktonic algae (Dokulil et al., 2000; Pfeifer, 2012). In addition, Microcystis could actively enrich phosphorus in nutrient-limited water by forming polyphosphate granules. They can be reproduced using such stored phosphorus under favorable conditions for growth (Jacobson and Halmann, 1982). Significantly, Microcystis overwinters in benthic sediments in the form of hypopus, outcompete in survival rate of other planktonic algae, and provide the inoculum for the next bloom (Verspagen et al., 2004, Misson et al., 2012). Several investigations have suggested that the germination of hypopus had provided the first propagules for development of new Microcystis populations (Preston et al., 1980, Head et al., 1998). Hence, recruitment from resting stages might be an important process of Microcystis population dynamics.

It is of particular importance to understand the recruitment of rest stages for the prediction and mitigation of Microcystis blooms. Recently, multiple abiotic and biotic factors have been recognized able to affect the recruitment of Microcystis. The abiotic factors, including altered temperature, light conditions, and concentrations of inorganic nutrients were contributed to the recruitment of Microcystis (Ståhl-Delbanco et al., 2003; Karlsson-Elfgren et al., 2004). Reynolds and Jaworski (1978) suggested that Microcystis would be active in sediments whentemperature of deep lakes is 7–8 ℃, and would grow rapidly at 15 ℃. Brunberg and Blomqvist (2003) indicated that enhanced light and temperature were important for triggering the growth of benthic Microcystis colonies, and the effects of biotic factorssuch as benthic fauna on the recruitment of Microcystis resting stages were also reported. The benthic invertebrates were found to stimulate the recruitment of freshwater cyanobacteria, by probably increasing nutrient availability in the pore water (Ståhl-Delbanco and Hansson, 2002). Sediment is a complex microenvironment that contains abundant benthic bacteria (Risgaard-Petersen et al., 2015; Zhou et al., 2017). Xing et al. (2011), Li et al. (2012), and Shao et al. (2014) reported significant association between bacterial communities and Microcystis blooms. Wu et al. (2019) found that the diversity and composition of the benthic microbial communities were changed after the addition of Microcystis biomass. In an earlier study of ours, three dominant strains of bacteria (named E.sp013, Ba.spD06, and Ba.spD24) were isolated from sediment of Chongtian Lake and identified able to accelerate the recovery of Microcystis aeruginosa (Zou et al., 2018). However, less is knownon the molecular response of Microcystis to those recruitment-promoting factors. In this study, we investigated the transcriptional response of M. aeruginosa triggered by benthic bacteria using theRNA-Seq analysis, and underlined the role of benthic bacteria in the recruitment process of M. aeruginosa in the early phase.

2 MATERIAL AND METHOD 2.1 Preparation of recruitment experiments

The Chongtian Lake (Changde City, Hunan Province, China) is often suffered from blooms caused by M. aeruginosa in the past few years (Zou et al., 2017). As described in Zou et al. (2018), three dominant benthic bacteria (E.sp013, Ba.spD06, and Ba.spD24) that separated from the sediment of Chongtian Lake were found able to induce the recruitment of M. aeruginosa. In this study, the simulation experiments on the recruitment of M. aeruginosa cells were performed in the presence orabsence of the dominant benthic bacteria. As described in our previous study (Zou et al., 2018), the sediment, lake water, and M. aeruginosa cell were sampled from the lake. The sediment was sampled using the KC Kajak sediment corer (6.5-cm diameter, KC, Denmark), and screened by a filter (sieve size 125 μm). The lake water was modified by adding NaNO₃ (0.16 g/L) and K₂HPO₄ (0.07 g/L) to provide enough nutrition for the growth of bacteria and M. aeruginosa. The M. aeruginosa cell was expanding incubated in BG11 liquid media (provided by the Chinese Academy of Sciences, 25 ℃, 50 μmol photons/(m²·s)), then placing at low temperature and dark for 3.5 months, and finally obtained the hypopus of M. aeruginosa. Three dominant benthic bacteria were cultured by the beef peptone liquid medium and mixed with the volume ratio of 1꞉1꞉1 (final volume: 10 mL), then added with M. aeruginosa cell solution in the same volume (10 mL). 10 mL of lake water instead of mixed bacteria was added as the control. Both the whole mixed solution and sediment were transferred into sterilized cylindrical glassware (SCG: diameter 15 cm, height 50 cm) and diluted to 7 L with sterilized lake water. All SCGs were placed in the artificial chamber (MMM Climacell, Germany) at the irradiance of 15 μmol photons/(m²·s), the temperature of 15 ℃, and illumination time 12 h꞉12 h (light꞉dark). The M. aeruginosa recruit for 3 days (control: CK3;treated: TR3), 5 days (control: CK5, treated: TR5), and 7 days (control: CK7; treated: TR7) were collected by centrifugation (2 000 r/min, 5 min) and used for RNA-Seq analysis.

2.2 Determination of chlorophyll a

The M. aeruginosa were collected from the sediment and overlying water with the different recruited days respectively, according to the method of Verspagen et al. (2004). The chlorophyll-a concentration of M. aeruginosa was determined according to the method of Berrendero et al. (2016). All data were presented as mean±SD, and conducted for one-way ANOVA by SPSS 17.0 software.

2.3 RNA-Seq analysis 2.3.1 Sequencing and annotation

RNA extraction and library construction were conducted using the TRzol^R Reagent (Thermo Fisher Scientific, USA) and NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, USA), respectively. All samples were sequenced on the Illumina sequencing platform (HiSeqTM 2500), and generated paired-end reads with a length of 150 bp/125 bp. Raw reads were assembled to transcript sequences by Trinity (Grabherr et al., 2011). The public databases were applied for functional annotation of unigenes, including Nr, Swiss-Prot, KEGG, and GO.

2.3.2 Gene expression analysis

All of the clean reads were mapped back onto the assembled transcriptome reference sequences through RSEM, after which the read count value of each gene was calculated (Li and Dewey, 2011). Then the FPKM algorithm was used to normalize the read count value of each gene (Trapnell et al., 2010). Finally, differentially expressed genes were screened by software DESeq based on the negative binomial distribution. A threshold for adjust P-value (q)<0.05 and absolute value of Log2 fold change value >1 was used to determine significant differences in gene expression. Moreover, the KEGG enrichment analysis was performed to describe the function of differentially expressed genes (DEGs) with threshold q<0.05.

2.3.3 Protein to protein interaction analysis (PPI)

The associations between the proteins and the proteins encoded by DEGs were determined using the STRING database (http://string-db.org/), with setting the string score at 700 (Seebacher and Gavin, 2011). Cytoscape software was applied for visualization and analysis of the network. The Cytohubba plugin of Cytoscape was applied to identify the hub genes in the network, based on the MCC algorithm. MOCDE plugin of Cytoscape was applied to screen modules from the PPI network with degree cutoff=2, node score cutoff=0.2, k-core=2, and max. Depth=100 (Talebi et al., 2018) Also, the functional modules were chosen with the cutoff values of intra connection nodes≥6 and node score≥4.0.

2.4 RT-qPCR analysis

A total of 5 differentially expressed transcription factors (TFs) were selected and quantified by real-time RT-PCR (RT-qPCR). The First Strand cDNA Synthesis Kit was used for inverse transcription according to the manufacturer's instructions (TINGEN, China). Real-time PCR was performed using the SYBR Premix EX Taq (TINGEN, China) by the ABI 7500 Real-time PCR System (Thermo Fisher Scientific, USA). In this study, the expression levels of the target unigene were normalized by 16S ribosomal RNA and quantified using the 2^-ΔΔCt method (Livak and Schmittgen, 2001).

3 RESULT 3.1 Overview of transcriptome from the M. aeruginosa

As described in our previous study (Zou et al., 2018), the concentration of chlorophyll a (Chl a) in the sediment was decreased with time, while those in the overlying water was increasing, which suggested that M. aeruginosa floated from sediment into the overlying water and proliferating rapidly there. The recovery efficiency of M. aeruginosa in the benthic bacteria-treated group was higher than that in control.

RNA-Seq analysis was conducted for M. aeruginosa at the early phase of recruitment(3–7 days). We obtained 5 803 803 unigenes with an average length of 404 bp. The N50 value of unigenes was 428 bp. The unigenes in length of 200 to 500-bp sequences represented the highest proportion, followed by 500 to 1 000-bp sequences (Fig. 1).

3.2 Gene expression and function analysis

Differential expression analysis was conducted between the samples from the different recovery days. For control ones, there were 1 and 76 843 unigenes up-regulated in the comparisons CK5 vs. CK3, and CK7 vs. CK5 respectively. For benthic bacteria treated ones, there were 71 570 and 104 351 unigenes up-regulated in the comparisons TR5 vs. TR3, and TR7 vs. TR5 respectively (q<0.05, |log2.Fold_ change|>1), showing the more active status in expression of genes. Results from the RT-qPCR showed that the expression of selected unigenes was consistent with transcriptome analysis (Supplementary Table S1).

Interestingly, a set of unigenes (54 982 unigenes) was up-regulated in the comparison CK7 vs. CK5 and was also up-regulated in comparison TR5 vs. TR3, which suggested that these unigenes would express earlier in M aeruginosa under the influence of benthic bacteria (named benthic bacteria responsive-unigenes, marked with grayness, Fig. 2a). KEGG enrichment analysis was conducted for these benthic bacteria responsive unigenes, found the "Ribosome" was the most enriched pathway and followed by "Cell cycle-Caulobacter" and "Oxidative phosphorylation". This suggested that functions of these unigenes were involved in the ribosome, mitosis, and energy metabolism of M. aeruginosa cells (Fig. 2b).

3.3 Protein to protein interaction analysis (PPI) for benthic bacteria responsive unigenes

The hub genes and core functions of benthic bacteria-responsive unigenes were characterized based on the PPI analysis. The PPI network was consisted of 1 092 interactions and 191 proteins, providing the overview of functional links between proteins encoded by those responsive unigenes (Fig. 3a). From the total PPI network, transcription factor sigma A (sigA, MAE_54470, marked with red), translation initiation factors IF-2 (MAE_14330) and IF-3 (MAE_30180), elongation factors lepA (MAE_39840), fusA (MAE_42760), and tufA (MAE_42770, marked with orange) were all up-regulated, which were essential regulatory factors for gene transcription and translation (Fig. 3a). Meanwhile, hub genes (top10) were identified including 30S ribosomal proteins and 50S ribosomal proteins and marked with red series among the network (Fig. 3b).

A module is a group of closely related proteins that act in concert to perform specific biological functions through the PPI network (Lin et al., 2015). There were 3 functional modules extracted from the PPI network (Fig. 3c). Among them, module 1 that contained translation initiation factors and elongation factors was the largest functional module with a significant score (31 proteins, 409 interactions). KEGG enrichment analysis was conducted for proteins from module 1, showed that the 'Ribosome'' was also the most enriched pathway (Table 1). This demonstrated that the core function of benthic bacteria-responsive unigenes was involved in ribosome in M.aeruginosa.

3.4 Expression analysis of Chl a synthesis related unigenes

In total of 5 enzymes that catalyzed the biosynthesis of Chl a were matched through the Swiss-Prot database, including magnesium chelatase (106 unigenes), magnesium-protoporphyrin O-methyltransferase (30 unigenes), anaerobic magnesium-protoporphyrin IX monomethyl ester cyclase (37 unigenes), divinyl chlorophyllide a 8-vinyl-reductase, chloroplastic (40 unigenes), and light-dependent protochlorophyllide reductase (17 unigenes, Fig. 4a). Results from the hierarchical cluster analysis showed that Chl a-related unigenes were more likely to high express in sample CK7 and TR5. These two samples were clustered together with similar expression patterns of unigenes. It may be contributed by synthetic activities of proteins from the ribosomes (Fig. 4b).

4 DISCUSSION

The cyanobacterium M. aeruginosa is a dominant species in many lakes, which tends to aggregate at the surface to form bloom on the water (Wilson et al., 2005). They were hibernated on lake sediments in winter, due to poor conditions like limited light and low temperature (Cirés et al., 2013). Bacteria constitute a crucial component of aquatic ecosystems, and the close interactions between the Microcystis bloom have received considerable attention (Grossart et al., 2006; Parveen et al., 2013; Landa et al., 2016). Buchan et al. (2014) indicated that Flavobacteria and Proteobacteria were dominant bacterial taxa that related to Microcystis blooms; Xu et al. (2018) suggested that Enterobacter aerogenes could enhance M. aeruginosa recovery, and increase the expressionof nicotinamide adenine dinucleotide (NADH) synthase and polysaccharide export-related genes. In our previous study, three benthic bacteria that promoted the recovery of M. aeruginosa were isolated from the sediment of Chongtian Lake (Zou et al., 2018). In this study, transcriptional response of M. aeruginosa to the recruitment promoting-benthic bacteria at the early stage of recovery was further investigated, underling the important role of benthic bacteria on the Microcystis bloom.

Ribosome was the ancient and highly conserved machinery that translates the genetic code into functional proteins. Studies suggested that the ribosomes of cyanobacteria would be in an inactive state during the stationary phase and in response to stress (ribosome hibernation), which was one prominent molecular strategy for post-stress survival (Wilson and Nierhaus, 2007; Trösch and Willmund, 2019). Hence, the activated ribosome was an essential part during the recovery of Microcystis. Here, a set of unigenes were identified to be expressed early during the recovery of M. aeruginosa, when in response to benthic bacteria. Results from the KEGG enrichment analysis showed that the function of these unigenes was related to ribosome. During the translation of M. aeruginosa cell, the mRNA first binds to the smallribosomal subunit (30S subunit), forming a complex together with the three initiation factors IF-1, IF-2, and IF-3, then the initiation complex interacts with the large ribosomal subunit (50S subunit) to assemble the full 70S ribosome. The peptide chain is continuously extended by the action of elongation factors lepA, fusA, and tufA. (Green and Noller, 1997). In this study, the hub genes of the network mediated by benthic bacteria-responsive unigenes all were ribosomal proteins of the 30S and 50S. Meanwhile, regulatory factors including translation initiation factors (IF-2 and IF-3) and elongation factors (lepA, fusA, and tufA) were up-regulated in the largest functional module that related to the ribosome. These demonstrated that the benthic bacteria had a positive influence on activating the ribosome during the early recovery stage of M. aeruginosa.

An active ribosome could provide abundant enzymes that catalyze various cellular processes. In this study, expression profiles of unigenes encoding the enzymes that catalyze the biosynthesis of Chl a was characterized by the heat map. These unigenes were more likely to be highly expressed in sample CK7 and TR5, which may be indirectly influenced by the activated ribosome in M. aeruginosa. This study underlined the role of benthic bacteria in the recovery of M. aeruginosa, and provided new thinking for understanding the regulation of Microcystis bloom.

5 CONCLUSION

In the presented study, we characterized the transcriptional response of M. aeruginosa to the recruitment promoting-benthic bacteria. There were 54 982 unigenes identified to be expressed earlier in response to benthic bacteria based on the expression level analysis. Results from the PPI analysis showed that hub genes of benthic bacteria-responsive unigenes mediated network were ribosomal proteins of 30S and 50S, and the most significant functional module of the network was also related to the ribosome. Unigenes encoding the translation initiation factors (IF-2, IF-3) and elongation factors (lepA, fusA, and tufA) were both up-regulated triggered by benthic bacteria. These indicated that the benthic bacteria had a positive influence on activating the ribosome during the early recovery stage of M. aeruginosa.

6 DATA AVAILABILITY STATEMENT

The datasets generated during the current study are available from the corresponding author on reasonable request.

Electronic supplementary material

Supplementary material (Supplementary Table S1) is available in the online version of this article at https://doi.org/10.1007/s00343-021-0387-4.