1 INTRODUCTION

As a result of global temperature and eutrophication, the Microcystis aeruginosa blooms are increasing in the freshwater ecosystems, becoming a global environmental problem (Wu et al., 2017). The direct effects of the mass growth of M. aeruginosa are unpleasant odor, loss of water transparency, and depletion of oxygen. This species produces microcystins (potent hepatotoxins to aquatic and terrestrial animals) that threaten public health regarding drinking water supply and recreational activities such as fishing (Otten and Paerl, 2015).

Many available methods tried to prevent the occurrence of blooms, including physical, chemical, and biological ways. One such a method is to reduce nutrient input into water systems to prevent nutrient enrichment and favor phytoplankton rather than cyanobacteria (Hrudey et al., 1999; Otten and Paerl, 2015), but this method is not particularly successful. Using chemical agents to inhibit or kill cyanobacteria may support ecological restoration and control the occurrence of blooms (Barrington and Ghadouani, 2008). Although the chemical methods may have serious negative impacts on the ecosystem, it has so far been considered an effective strategy to control M. aeruginosa blooms (Zhou et al., 2013; Wu et al., 2017).

At present, copper sulfate (CuSO₄) is a low-cost chemical reagent used to inhibit M. aeruginosa blooms, on which many papers have been published documenting toxicity of Cu²⁺ to Microcystis (Hadjoudja et al., 2009; Wu et al., 2017). Leng et al. (2015) reported that Cu²⁺ could cause rapid and excessive accumulation of reactive oxygen species (ROS) in M. aeruginosa and disrupt membrane integrity. Rascioa and Navari-Izzo (2011) suggested that the Cu²⁺ toxicity to plants might be caused by inactivating enzymes, blocking of functional groups of metabolically important molecules, and/or displacing or substituting for essential elements in numerous physiological processes. However, to the best of our knowledge, the published materials did not address the mechanisms of Cu²⁺ toxicity to M. aeruginosa from a transcriptome perspective. Presently, availability of complete M. aeruginosa genome and high-throughput RNA sequencing (RNA-seq) techniques offer tools for the study on the mechanisms of Cu-dependent stress response in M. aeruginosa (Kaneko et al., 2007; Wang et al., 2017).

A moderate dose of CuSO₄ can be used to inhibit the growth of Microcystis and control bloom formation, which is non-toxic to humans and aquatic animals (Zhang et al., 2007; Barrington and Ghadouani, 2008). Tsai (2015) reported that Cu²⁺ concentration of 0.16 mg/L is the minimum concentration required to inhibit the growth of M. aeruginosa. Zhang et al. (2007) found the 96-h EC₅₀ of M. aeruginosa was 1.02 mg/L Cu²⁺. The Cu²⁺ doses below 2 mg/L are considered safe according to the World Health Organization recommendation.

In this study, M. aeruginosa was treated in 0.5 mg/L Cu²⁺ (half of the 96-h EC₅₀ according to Zhang et al., 2007) for 72 h. We characterized the toxicological effects of Cu²⁺ at a physiological level. The goals of this study are to construct a M. aeruginosa transcriptome database, to screen the significantly expressed genes under Cu-toxic conditions, to understand Cu²⁺ toxicity mechanism in M. aeruginosa, and to provide a genetic resource of M. aeruginosa.

2 MATERIAL AND METHOD 2.1 M. aeruginosa culture and Cu²⁺ treatment

Microcystis aeruginosa was purchased from the Institute of Hydrobiology, Chinese Academy of Sciences. M. aeruginosa were cultured at 28℃ under 2 500-lx lighting in 12 h/12 h light/dark cycle in BG- 11 medium. The initial concentration of M. aeruginosa in both the control and the treatment was 10⁶ cells/mL, and 0.5 mg Cu²⁺ (as CuSO₄) per liter was added to the treatment group; the study lasted for 72 h. Each treatment was replicated three times. In this experiment, we aimed to monitor the Cu²⁺ stress responses rather than the algae death, so that the treatment of Cu²⁺ concentration was chosen carefully. In previous experiments, we found that 0.5 mg/L Cu²⁺ was the appropriate concentration, resulting in a sharp decline in several physiological traits, while the mortality of M. aeruginosa was below 20%.

2.2 Biochemical analysis

Photosynthetic efficiency was evaluated by using a handy photosynthetic efficiency analyzer (HandyPEA, Hansatech Instruments Ltd, Norfolk, UK). Parameters F_v (variation in fluorescence) and F_m (maximal fluorescence in the dark-adapted state) were measured at 0, 24, 48, and 72 h after commencement of the Cu²⁺ treatment. Subsequently, the value of F_v/F_m was calculated.

The ROS concentration was measured at 0, 24, 48, and 72 h according to Guo et al. (2016a) after staining the cells by DHR123 (dihydrorhodamine 123) for 1 h. The final concentration of DHR123 in the cell cultures was 5 μmol/L. After incubation, fluorescence (for excitation/emission at a wavelength of 485/530 nm) was measured using a fluorescence spectrophotometer (Model LS-55, PerkinElmer, CA, USA) to quantify the ROS concentration.

At the end of the experiment, chlorophyll a (Chl a) concentration was determined at room temperature using a PHYTO-PAM phytoplankton analyzer (Water-PAM, Walz, Germany). This was performed after keeping the cells in the dark for 7.5 min prior to measurement based on the published method (Gera et al., 2012)

Copper concentration in M. aeruginosa was measured as per Foster et al. (2008). Briefly, approximately 50.0 mg freeze-dried algal sample was digested with 2 mL HNO₃ at 120℃. Total Cu²⁺ was determined by inductively coupled plasma-optical emission spectroscopy (ICP-OES, Optima 7000, PerkinElmer, USA).

Antioxidant enzyme activities were measured using the commercial kits (purchased from Nanjing Jiancheng Bioengineering Research Institute of China, Nanjing, China) (Gui et al., 2010). The total superoxide dismutase (SOD) assay was based on the SOD-mediated inhibition of nitrite formation from hydroxylammonium in the presence of an O₂^- generator (Hong et al., 2008). The catalase (CAT) activity was defined as the amount of enzyme required to transform 1-μmol H₂O₂/min at 28℃ and pH 7.0. The peroxidase (POD) activity was defined as a 0.01-absorbance unit change per min at 470 nm. The malondialdehyde (MDA) content was measured based on the chemical reaction between MDA and 2-thiobarbituric acid (TBA), which is expressed as nmol MDA/mg protein (Ohkawa et al., 1979). Total soluble protein was measured using the Bradford method (Bradford, 1976).

2.3 Transcriptome analysis 2.3.1 RNA extraction, cDNA library construction and Illumina deep sequencing

Total RNA samples of Cu²⁺-treated and control group were extracted using commercial kits (Trizol reagent, Invitrogen, USA), and subsequently used for mRNA purification and library construction following the manufacturer's instructions (the Ultra RNA library prep kit for Illumina, NEB, USA). Then, the samples were sequenced on an Illumina Hiseq 2500, and each sample yielded > 2 Gb of data. Sequencing was completed by the Nanjing Genepioneer Biotechnologies Company (Nanjing, Jiangsu, China). Three repetitions were used for RNA-Sequencing.

2.3.2 Transcriptome data analysis

The procedure of transcriptome data analysis was as follows. Firstly, low quality reads and reads containing adapter or ploy-N were removed from the raw data, generating the high-quality clean data. Secondly, the clean reads were mapped to the M. aeruginosa NIES-843 genome. Thirdly, the numbers of reads mapping to each gene were counted, and FPKM of each transcript was calculated according to the published method (Trapnell et al., 2010). Fourthly, differential expression genes (DEGs) were analyzed by the DESeq R package (Anders and Huber, 2010). The significant DEGs among different samples were identified with P < 0.05 and log2 foldchange > 1. Finally, we used the GOseq R package, COGPipeline and KOBAS software to implement, respectively, the Gene Ontology (GO), Cluster of Orthologous Groups of proteins (COG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DEGs. In our study, the degree of KEGG enrichment was measured using the Richfactor, Q-value, and the number of genes enriched in a KEGG pathway. The Richfactor refers to the ratio of the number of differentially expressed to the total number of annotated genes. The larger the Richfactor value, the greater the degree of enrichment. The Q-value is the P-value after multi-hypotheses test correction in a range of 0–0.05; the closer the value is to zero, the more significant the enrichment (Wang et al., 2018).

2.4 Real-time PCR

Purified RNA was reverse-transcribed into cDNA using a First Strand cDNA Synthesis Kit according to the manufacturer's instructions (TaKaRa, Japan). To adjust the quantity of input cDNA, 16S rRNA was used as the internal control (Chen et al., 2015). The primers (as listed in Table 1) for RT-PCR were designed according to the known sequences of M. aeruginosa in National Center for Biotechnology Information (NCBI).

RT-PCR was performed using the SYBR Premix EX Taq (TaKaRa, Japan) on an ABI 7500 Real-Time PCR System. An aliquot of 2.0-μL template cDNA was added to the final volume of 20-μL reaction mixture. The RT-PCR protocol consisted of initial denaturation at 95℃ for 30 s, followed by 35 cycles of 95℃ for 3 s each, and annealing at 60℃ for 34 s. The final dissociation stage included one cycle of 95℃ for 15 s, 60℃ for 1 min and 95℃ for 15 s. In this study, we calculated the relative quantity of the target gene transcripts with respect to the 16S rRNA (Chen et al., 2015) using the 2^-ΔΔCt method.

3 RESULT 3.1 Physiological parameters of M. aeruginosa

As shown in Table 2, Cu²⁺ content in M. aeruginosa increased after Cu²⁺ exposure for 72 h, while Chl a content in M. aeruginosa sharply decreased. Activities of antioxidant enzymes (SOD, CAT and POD) and MDA concentration were all increased after the 72-h Cu²⁺ exposure.

3.2 Photosynthetic efficiency and ROS accumulation of M. aeruginosa

The value of F_v/F_m was decreased in the Cu²⁺ treatment with increasing exposure time (Fig. 1a). The cellular ROS production was detected by measuring the specific fluorescence intensity that increased with increasing exposure time in the Cu²⁺ treatment (Fig. 1b).

3.3 Overview of the sequencing

In total, there were 1 994 017 500 raw reads generated from the Cu²⁺-treated and 2 364 087 300 from the control library (Table 3). After discarding the low quality reads, 6 646 725 and 7 880 91 clean reads were obtained for the Cu²⁺-treated and control libraries, respectively. Of these clean reads, the Q30 percentage was above 84% and GC content was approximately 46% for the two libraries.

3.4 Functional annotation of the DEGs

For annotation, all DEGs were aligned to databases such as NR, Swiss-Prot, GO, KEGG, and COG by blastx with an E-value cut-off of 10^-5. The results indicate that 594 DEGs in the transcriptome libraries were successfully annotated in at least one of the above databases. Specifically, all DEGs were annotated in the NR database, whereas only 292 had the significant matches in KEGG (Table 4).

3.5 Functional classification of transcripts by GO

GO is a useful program for annotating and analyzing the functional categorization of annotated genes (Conesa et al., 2005). The total of 595 genes (343 up-regulated and 252 down-regulated) with the GO annotations were classified into the following three functional categories: cellular component, molecular function and biological process. In the cellular component category, the majority of the annotated DEGs were assigned to cell, cell part, and membrane. In the molecular function category, the widest distribution of the DEGs was associated with catalytic activity and binding. The majority of the DEGs were involved in the metabolic, cellular and single-organism processes in the category biological process (Fig. 2).

3.6 Functional classification of transcripts by COG

COG is a database for the classification of orthologous gene products. The COG analysis of the DEGs provided information about their possible functions. The total of 418 DEGs with 25 functional definitions were obtained (Fig. 3). Among them, 'replication, recombination and repair', 'energy production and conversion', and 'general function prediction only' were the largest three groups of transcripts.

3.7 KEGG pathway analysis

As shown in Fig. 4, the KEGG enrichment analyses of the DEGs showed that the majority were assigned to three specific KEGG categories, including "environmental information processing", "metabolism" and "genetic information processing". In addition, the KEGG pathway enrichment analysis showed that, among the top 10 enriched pathway for the Cu²⁺-treated vs control, the "photosynthesis" and "oxidative phosphorylation" pathway was significantly enriched (with Q-value < 0.05 and the largest number of genes enriched).

3.8 The genes related to photosynthesis and oxidative phosphorylation based on KEGG enrichment

In the present study, 29 genes related to photosynthesis (Table 5) and 17 genes related to oxidative phosphorylation (Table 6) were identified. Among 29 down-regulated photosynthesis genes, most genes were involved in photosystem I (PS I) (4 genes) and photosystem Ⅱ (PSII) (12 genes). Among 17 genes related to oxidative phosphorylation, the ndh series genes (including that coding for NADH dehydrogenase) were the most prevalent.

3.9 Validation of RNA-Seq by real-time PCR

Real-time PCR was used to validate the mRNA expression of eight genes (from those included in Tables 5 and 6). Regardless of the Cu²⁺ treatment, the mRNAs of rbcL, psbA, psbD, psaC, ATPb, and ndhC were all significantly down-regulated compared to the control. While, the mRNAs of Cyt-c and SDH were significantly up-regulated compared to the control. Hence, the expression patterns were similar to those in the transcriptome (Fig. 5).

3.10 RT-PCR analysis of eight genes during periods of Cu²⁺ stress

In M. aeruginosa, over the entire 72 h of Cu²⁺ stress, the 8 validated genes displayed various expression patterns, and the transcription of rbcL, psbA, psbD, psaC, ATPb, and ndhC were all decreased at the tested time points. However, the transcription of Cyt-c and SDH were enhanced during the periods of Cu²⁺ stress (Fig. 6).

4 DISCUSSION

There are few related reports on the influence of chemical reagents on M. aeruginosa, and most of these papers are focused on the analysis of physiological parameters such as algae mortality, chlorophyll content and antioxidant enzymes (Hong et al., 2008; Lu et al., 2014; Chen et al., 2015). In this study, we investigated the effect of Cu²⁺ on the above physiological parameters of M. aeruginosa, and on its transcriptome too.

Copper is a micronutrient essential for the growth of organisms by participating in many physiological processes, such as the photosynthetic and respiratory electron-transport chains, the C and N metabolism, hormone perception, cell wall metabolism, and oxidative stress protection. However, Cu²⁺ poisoning can occur if organisms take excessive amounts of Cu²⁺ (Garcia-Molina et al., 2011). In the present study, Cu²⁺ content in M. aeruginosa increased by 1 356% after Cu²⁺ exposure for 72 h, with Chl a content in M. aeruginosa sharply decreasing and antioxidant enzymes (SOD, CAT and POD) all increasing. The results indicate that the Cu²⁺ concentration in BG-11 medium was above the M. aeruginosa growth needs, and led to an oxidative stress in M. aeruginosa. Similar results were reported by Lu et al. (2015) and Wu et al. (2017).

Photosynthesis is one of the most important activities in cyanobacterial cells and is also regarded as a major target of allelopathic substances. In this study, the F_v/F_m ratio indicates the efficiency of photosynthetic apparatus (Shao et al., 2011) by representing the maximum quantum efficiency of PSⅡ photochemistry (Baker, 2008). We find the F_v/F_m value was decreased in the Cu²⁺ treatment with increasing exposure time, indicating Cu²⁺ was responsible for the inhibition of photosynthesis. Similar results were reported by Boucher and Carpentier (1999) and Chen et al. (2015) using different metals (e.g. Hg²⁺, Pb²⁺, and Na²⁺).

In order to understand the potential regulatory mechanisms of Cu²⁺ toxicity stress in M. aeruginosa, RNA-seq technology was used as an important tool for global analysis of gene discovery and molecular mechanisms under environmental stress (Wang et al., 2017). With RNA-seq technology, we sequenced the transcriptomes of Cu²⁺-treated and untreated M. aeruginosa, respectively, and obtained 6 646 725 and 7 880 291 clean reads matched to one or more public databases. According to the results of GO classification in our study, most DEGs were involved in the important category of biological process, such as metabolic process, cellular process and singleorganism process. These findings are similar to the result of Wang et al. (2016). However, in another study, 'metallochaperone activity' and 'translation regulator activity' were the most representative (Gai et al., 2012), two categories containing the fewest genes in the study presented here. This difference may be on due to different species used. In addition, the second explanation for functional class differences may be better gene sequence information, existing a series of low-abundance transcripts, being provided by Illumina sequencing platform (this study) than 454 GS-FLX platform (Gai et al., 2012).

In addition to GO analysis, assignments of COG were performed to reveal the diversity of the assembled transcripts while reflecting the global landscape of the transcriptome (Wu et al., 2014; Wang et al., 2016). In this paper, a total of 418 DEGs were obtained, of which 25 have functional definitions. 'Replication, recombination and repair', 'energy production and conversion', and 'general function prediction only' were the largest three groups of transcripts. Hanawalt (2007) reported that the replication and maintenance of the genome are absolute requirements for life. A cell must duplicate its entire complement of DNA with astonishing precision in the face of a barrage of deleterious endogenous and environmental genotoxic agents, as well as the intrinsic chemical instability of the DNA molecule itself. In our study, many DEGs might still need to be functionally verified in our findings, as nucleotide sequencing and protein metabolism of M. aeruginosa could be active.

The KEGG pathway database involves a systematic analysis of inner-cell metabolic pathways and functions of gene products (Wang et al., 2016). In this study, KEGG pathway analysis annotated about half of DEGs, whereas NR database annotated all DEGs identified. Similar results were reported by Wang et al. (2016), in which more DEGs were identified in NR database than other database. However, KEGG pathway database is the widely used database to find enriched pathway. KEGG pathway database is a valuable collection of manually constructed maps representing the current knowledge on the molecular interaction and reaction networks for metabolism, processing of genetic and environmental information, cellular processes, etc. Moreover, KEGG (as a part of a system biology approach) can be considered as a virtual biological system including various types of information essential for the recreation of a living creature (Kanehisa and Goto, 2000). The results of this study indicate that most affected metabolic pathways were in photosynthesis and oxidative phosphorylation. Among the down-regulated photosynthesis genes, most genes were involved in photosystem I (PS I) and photosystem Ⅱ (PSII). Zou et al. (2015) found PS I and PSII in photosynthetic organisms are sensitive to environmental variation. Our transcriptome and RT-PCR results proved that core proteins expressions in the PS I and PSII reaction centers was depressed in response to Cu²⁺ stress. These findings are consistent with Lu et al. (2014), who examined the effect of the epigallocatechin-3- gallate treatment on the gene expression of M. aeruginosa.

Through the research, the 29 genes related to photosynthesis can be considered a target for the subsequent Cu²⁺ stress analysis. Among 16 genes involved in PS I and PSII is the psbA gene that encodes D1 protein as the key protein in the PSII reaction center. The D1 protein is involved in the PSⅡ repair cycle by binding to all the redox active components, while participating in the PSII repair cycle with degrading and resynthesizing in a multistep process (Leu et al., 2002; Leng et al., 2015). As Mulo et al. (2012) reported, in this study the time-dependent down regulation of psbA transcription can impede electron transport and suppress the electron transport activity in PSⅡ. The rubisco protein is encoded by the rbcL gene and acts as a key protein in the lightindependent reactions, involved in the carbon assimilation process (Mulo et al., 2012; Chen et al., 2015). The experimental results reveal that electron transport, energy conversion and photosynthesis of M. aeruginosa are necessary for CO₂ fixation, but they are inhibited under Cu²⁺ conditions. In summary, the hindrance of carbon assimilation under Cu²⁺ stress suppressed the photosynthesis of M. aeruginosa. Interestingly, 12 differentially expressed transcripts linked to PSII compared to only 4 in PS I, demonstrating that PSII was more sensitive than PSI to Cu²⁺ toxicity stress, similar results were reported by Boucher and Carpentier (1999) and Leng et al. (2015).

Oxidative phosphorylation is a metabolic pathway in which cells use enzymes to oxidize nutrients, and the energy released during oxidation of the substance in the body provides the coupling reaction process for the synthesis of ATP from ADP and inorganic phosphate (Powell and Somero, 1986). This metabolic pathway was taken place in the mitochondria and releases energy more efficiently than alternative fermentation processes such as anaerobic glycolysis, so that oxidative phosphorylation process is present in most aerobic organisms (Guo et al., 2016b). Although oxidative phosphorylation pathway is a primary part of the metabolism, it produces ROS, which in turn leads to the proliferation of radicals, destroying cells and causing aging and disease (Kuno et al., 1985; Powell and Somero, 1986). In the present study, ROS concentration and antioxidant enzyme activity (e.g. SOD, CAT) increased in M. aeruginosa after 72 h Cu²⁺ exposure (as above), with excess ROS causing cell damages. Leng et al. (2015) reported NADPH oxidase has been suggested to play a crucial role in the generation and accumulation of ROS in cyanobacteria under stress conditions.

It was worth nothing that the 17 differentially expressed genes related to oxidative phosphorylation could be used as the targets genes for oxidative analysis. Among these genes, ndh series genes represented the largest proportion. The ndh could encode NADH dehydrogenase, a main component of the electron transport chain that transfers electrons from NADH to coenzyme Q (Lu et al., 2014). In this study, the expression of ndhC gene was timedependent decreased, which might give rise to serious damage to the cells. Furthermore, the ATP synthase genes were also time-dependent down (such as ATPb). Lu et al. (2014) pointed out that the down-regulation of photosynthesis and ATP synthase were markers of the photoinhibition.

5 CONCLUSION

This study discovered numerous genes that were differentially expressed between the control and Cu²⁺- toxicity-stressed M. aeruginosa. Transcriptome analysis results indicate that exogenous Cu²⁺ regulated numerous genes related to photosynthesis and oxidative phosphorylation. The photosynthesis in M. aeruginosa was hindered by changing carbon assimilation under Cu²⁺ stress. Our research will help understand and explain the complex molecular mechanisms of M. aeruginosa responses to Cu²⁺ toxicity stress.

6 DATA AVAILABILITY STATEMENT

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

7 ACKNOWLEDGEMENT

We thank Genepioneer Biotechnologies (Nanjing, China) for the assistance in analysis of the RNA-seq data.