Synechocystis sp. PCC 6803 was obtained from American Type Culture Collection (ATCC), and used as a wild type to construct single-gene knockout mutants of TR genes. A total of 32 knockout mutants of putative TR-coding genes were constructed, confirmed and described previously (Zhu et al., 2015). Briefly, for the gene target selected, three sets of primers were designed to amplify a linear DNA fragment containing the chloramphenicol resistance cassette (amplified from a plasmid pACYC184) with two flanking arms of DNA upstream and downstream of the targeted gene. The linear fused PCR amplicon was used directly for transformation into Synechocystis by natural transformation. The chloramphenicol resistant transformants were obtained and passed several times on fresh BG11 plates supplemented with 10 μg mL^-1 chloramphenicol to achieve full chromosome segregation (confirmed by PCR). The mutants and the wild type were grown in the BG11 medium (pH 7.5) in 100-mL flasks each with 25 mL medium, the light intensity was approximately 50 μmol photons m^-2s^-1 and the illuminating incubator was 130 rpm, the temperature was controlled at 30°C (HNY-211B Illuminating Shaker, Honour, China). All mutants were first cultivated in BG11 culture with 10 μg/ml chloromycetin for 48 h and then inoculated into BG11 culture without chloromycetin. The growth was determined by cell density measured at OD₆₃₀ on a UV-1750 spectrophotometer (Shimadzu, Japan) every 12 h. For each mutant, three biological replicates were established independently, and each sample was measured in triplicates (Zhu et al., 2015). To confirm the growth patterns, growth experiment of every knockout mutant was repeated at least three times independently, then the growth rates of all mutants were calculated (Supplementary Table S1). Only the growth rates between different RR mutants and wild type with p-value <0.005 by t-test were considered a significant growth difference.

Liquid chromatography-mass spectrometry based targeted metabolomics was performed according to the protocol described previously (Wang et al., 2014). All chemicals used for LC-MS metabolomics analysis were obtained from Sigma–Aldrich (Taufkirchen, Germany). For metabolomic analysis, the wild type and the mutant cells were collected at 48 and 72 h, respectively, and each sample was prepared with three biological replicates. Due to the large amount of cultivation needed to finish the comparative experiments of 32 TRs, the samples had to be cultivated, prepared, and analyzed in five batches. A separate cultivation and analysis of the wild type as control was conducted for every batch to minimize possible batch difference. Briefly, the cells were collected by centrifugation at 7500 × g for 8 min at 4°C (Eppendorf 5430R, Hamburg, Germany), quenched, and extracted rapidly with 900 μL of 80:20 MeOH/H₂O (-80°C) and then frozen in liquid nitrogen. The samples were then frozen-thawed three times to release metabolites from the cells. The supernatant was collected after centrifugation at 15,000 × g for 5 min at -4°C and then stored at -80°C. The remaining cell pellets were re-suspended in 500 μL of 80:20 MeOH/H₂O (-80°C), and then the above extraction process was repeated. The supernatant from the second extraction was pooled with that from the first extraction and stored at -80°C until the LC-MS analysis was conducted. LC-MS analysis was conducted on an Agilent 1260 series binary HPLC system (Agilent Technologies, Waldbronn, Germany) using a Synergi Hydro-RP (C18) 150 mm × 2.0 mm ID, 4-μm 80-Å particle column (Phenomenex, Torrance, CA, USA), coupled to an Agilent 6410 triple quadrupole mass analyzer equipped with an electrospray ionization (ESI) source. Data was acquired using Agilent Mass Hunter workstation LC/QQQ acquisition software (version B.04.01), and chromatographic peaks were subsequently integrated via Agilent Qualitative Analysis software (version B.04.00). A total of 24 metabolites were selected for LC-MS-based targeted metabolite analysis in this study. All data of metabolomic profiling was first normalized by the internal control and the cell numbers of the samples.

The 24 targeted metabolites include acetyl coenzyme A (AcCOA), adenosine 5′-diphosphate (ADP), adenosine-5′-diphosphoglucose (ADP-GCS), α-ketoglutaric acid (AKG), adenosine 5′-monophosphate (AMP), adenosine 5′-triphosphate (ATP), coenzyme A hydrate (CoA), dihydroxyacetone phosphate (DHAP), D-fructose 1,6-bisphosphate (FBP), D-fructose 6-phosphate (F6P), sodium fumarate dibasic (FUM), DL-glyceraldehyde 3-phosphate (GAP), D-glucose 6-phosphate (G6P), L-glutamic acid (GLU), α-nicotinamide adenine dinucleotide (NAD), reduced α-nicotinamide adenine dinucleotide (NADH), nicotinamide adenine dinucleotide phosphate (NADP), reduced nicotinamide adenine dinucleotide phosphate (NADPH), uridine 5′-diphosphoglucose (UDP-GCS), oxaloacetic acid (OXA), phosphor (enol)pyruvic acid (PEP), D-(-)-3-phosphoglyceric acid (3PG), D-ribose 5-phosphate (R5P), and D-ribulose1,5-bisphosphate (UDP-GCS), uridine 5′-diphosphoglucose (RiBP).

Protein sequences of all the 32 TR genes were downloaded from NCBI¹. To define potential TR families, we used BLAST software for homology identification, only those TR with 80% aligned coverage with E-value < 1e-20 were consider as same families.

A protein–protein interaction (PPI) dataset of Synechocystis was downloaded from the STRING database (http://www.string-db.org/) (Jensen et al., 2009). STRING aggregates data and predictions stemming from a wide spectrum of cell types and environmental conditions, and aims to represent the union of all possible protein–protein links. In the STRING database, several types of evidence for the association, including genomic context, high-throughput experiments, conserved co-expression and previous biological knowledge were used to calculate a single combined score for each gene in the genome. In this study, only those experimentally validated were applied to construct the PPI network to cover potential protein-protein connections, and the notes of all the proteins in this study were renamed using gene IDs (Szklarczyk et al., 2011).

Although bioinformatics analysis of TRs in Synechocystis based on the sequence similarity has been conducted previously (Huffman and Brennan, 2002; Los et al., 2010), their functional classification using experimental approaches is still insufficient. A library of single deletion mutants for 32 TR-encoding genes in Synechocystis was constructed and confirmed previously in our laboratory (Zhu et al., 2015), and the majority of them have not been functionally characterized. To seek more functional information for these Synechocystis TRs, we first measured differential growth of all TR mutants in normal BG11 medium in flask cultivation, in parallel with the wild type Synechocystis. While most of the mutants grew equally well as the wild type (Supplementary Figure S2), the comparative analysis showed that four TR mutants, Δslr0395, Δslr1871, Δslr0741, and Δsll1872, grew poorly in the BG11 medium when compared with the wild type (Figure 1 and Supplementary Table S1), suggesting that the function of these four TRs might be related to key metabolism necessary for normal growth in the BG11 autotrophic growth medium. Among them, deletion of slr0395 (ntcB) caused the most significant growth arrest, with approximately 40% of the growth compared to the wild type after cultivation of 72 h, while deletion of slr1871 (pirR) resulted in the least growth defect, with only 15% of the growth compared to the wild type after 72 h cultivation.

Slr0395 has previously been annotated as nitrite-responsive transcriptional enhancer NtcB in Synechocystis, on the basis of the inability of the Δslr0395 mutant to rapidly accumulate the transcripts of the nitrate assimilation genes upon induction and to respond to nitrite. In the ntcB mutant, activities of the nitrate assimilation enzymes were 40 to 50% of the wild type level, and the cells grew on nitrate at a rate approximately threefold lower than that of the wild type (Aichi et al., 2001).

Slr1871 was previously annotated as PirR of a LysR family, whose encoding gene is located immediately upstream of pirAB encoding an ortholog of pirin in the divergent direction, and DNA microarray analysis indicated that PirR repressed expression of closely located ORFs, slr1870 and mutS (sll1772), in addition to pirAB and pirR itself (Hihara et al., 2004). Slr0741 was previously found to encode a negative regulator of the Pi regulon and its insertional inactivation in Synechocystis led to increase of the intracellular polyP level (Morohoshi et al., 2002); Slr0741 was also found involved in transduction of the phosphate-limitation signal in Synechocystis (Juntarajumnong et al., 2007), and it was also up-regulated upon ethanol stress as revealed by RNA-seq analysis (Wang et al., 2012). Currently no functional information is available for Sll1872.

Liquid chromatography-mass spectrometry-based metabolomics analysis has been recently used to investigate cyanobacterial metabolism due to its advantages toward chemically unstable metabolites, such as the hydrolytically unstable nucleotides (i.e., ATP, GTP, cAMP, and PEP) and the redox active nucleotides (i.e., NADPH, NADP) whose determination could be important in deciphering metabolic responses to genetic or physiological changes. Using a protocol optimized in our previous studies (Wang et al., 2014), a LC-MS-based comparative metabolomics analysis was conducted on all TR mutants and the wild type, with 24 key metabolites involved in central carbon metabolism, cellular energy charge and redox monitored in all samples at 48 and 72 h. The cell samples of 32 TR mutants and the wild type used for LC-MS-based metabolomic analysis were cultivated in BG11 media under autotrophic growth condition and collected at both 48 and 72 h, which were corresponding to earlier and latter exponential phases of cell growth. Each sample was prepared with three biological replicates. As for some, metabolite levels might change during the 8 min centrifugation procedure, a faster centrifugation method may be considered. To reduce the sampling time and maintain the metabolites as much as possible, a higher rotation rate and a lower temperature should be considered. However, with higher rotation rate, more severe damages may occur to cells. A protection agent may be considered in the future. When the rotation rate was higher than 7500 × g, the cells of Synechocystis would be easily broken, leading to less metabolites being preserved.

Several patterns were observed in the PCA plots of the metabolomic data (Figure 2): (i) As a large number of cultivation was needed to finish the comparative experiments of 32 TRs, the samples had to be cultivated, prepared and analyzed in five batches. To minimize possible batch difference, a separate cultivation and analysis of the wild type as control was conducted for every batch. In the PCA plots, the big dots of five different colors representing the controls of five batches were found clustered together after data normalization, demonstrating the systematic errors resulting from the experimental design and different cultivation batches were not significant (Figure 2); (ii) metabolic profiles of TR mutants were in general well separated at both time points, demonstrating that the LC-MS-based methodology we utilized in this study is sensitive enough to investigate possible differences between controls and all the TR mutants; (iii) except for the four mutants that were grew poorer than the wild type, almost no growth difference was observed between the remaining 28 TR mutants and the wild type when cultivated in the BG11 medium; however, PCA analysis of metabolic profiles showed that these TR mutants were well separated from the wild type in the plots, suggesting that the deletion of these TR-encoding genes has caused significant changes to the cells at the metabolite level. At 48 h, seven TR mutants with the most significant metabolic changes from the wild type control were Δslr1871, Δsll1392, Δsll0690, Δslr1666, Δsll1872, Δsll0782, and Δsll0792 (Figure 2A, with two components proportions of 15.94 and 11.23%); while at 72 h, six TR mutants, Δslr1937, Δsll1594, Δslr1529, Δsll0792, Δslr1871, and Δslr0449 (Figure 2B, with two components proportions of 14.70 and 13.10%), displayed the most significant metabolic changes, suggesting a significant functional diversity of TRs in Synechocystis, as revealed by the metabolic profiling of selected metabolites related to central carbohydrate metabolism; (iv) only two TR mutants, Δsll0792 andΔslr1871, were found significantly regulated at both 48 and 72 h at the metabolite level, suggesting that time- or growth phase-dependent regulation may be involved for most of the responsive TRs; (v) a close examination of the mutant Δslr1871 showed that almost all the metabolites involved in central carbohydrate metabolism were down-regulated when compared with the wild type, consistent with the previous results that the slr1871 (pirR) gene had a reduced transcript level during a light-limited linear growth when compared to the exponential growth (Foster et al., 2007). Although it has been reported that the sll0792 gene encodes ZiaR, a Zn²⁺-responsive repressor of ziaA encoding a polypeptide with sequence features of heavy metal transporting P-type ATPases in Synechocystis (Thelwell et al., 1998), its regulatory function on cellar metabolism has not yet been established.

To confirm the analysis with PCA, another approach, the Euclidean distances calculated based on the different metabolite profiles between each mutant and its control (the wild type), were also determined. In Supplementary Figure S1, the TR mutants were ordered based on the degree of metabolic changes when compared with the wild type in descending order. The analysis was conducted separately with metabolomic profiling data of two time points (i.e., 48 and 72 h), and showed that the top changed mutants were Δsll0690, Δslr1871, Δsll1670, Δslr1666, Δslr1245, Δsll1957, Δsll1594, and Δsll1872 at 48 h; Δslr0449, Δslr1871, Δssl0564, Δsll1937, Δslr1666, Δslr0115, Δslr0724, and Δsll0690 at 72 h, respectively. The results were consistent with those of the PCA analysis.

Comparative growth analysis showed that four TR mutants were grown poorly in BG11 medium when compared with the wild type. A detailed analysis of the metabolite abundance of the 24 metabolites in these mutants was then conducted (Figure 3). Slr0395 is involved in regulation of nitrate assimilation gene (ntcB) (Burnap et al., 2015), so it is expected that the deletion of the slr0395 (ntcB) gene would decrease nitrogen metabolism. Accordingly, the metabolomics analysis showed Glu was up-regulated by 24.2, and 11.7% at 48 h and 72 h, respectively. In addition, abundance of CoA in the Δslr0395 was found increased.

Slr0395 and Slr1871 are LysR-type transcriptional regulator protein (LTTRs) that have been found important for regulation of the carbon concentration mechanism (CCM) in cyanobacteria (Daley et al., 2012). In Calvin–Benson cycle, F6P and G3P were converted to GAP, R5P, and then Ru5P and RiBP for CO₂ fixation (Wang et al., 2011). The metabolomic analysis showed that metabolite RiBP, F6P, GAP, and R5P associated with RiBP were down-regulated in the Δslr0395 and Δslr1871 mutants (Figure 3), which could be responsible for the slow growth in the mutants.

Regulatory function of gene sll1872 (lytR) has not been determined previously. According to our metabolomic analysis, several metabolites related to energy metabolism, including NADP, NADH, NAD and ATP were all down-regulated at both sampling time points, suggesting that the absence of the sll1872 (lytR) may negatively affect the energy metabolism and then reduced the growth of Δsll1872 mutant (Figure 3).

Interestingly, although the growth of Δslr0741 was clearly slower than the wild type, the metabolomic analysis showed no obvious difference at the metabolite level between the Δslr0741 mutants and the wild type, implying that its regulatory function may be not directly related to the central carbohydrate metabolism.

The finding that a range of metabolic changes occurred in single-deletion TR mutants led to several immediate questions. First, whether the TR mutants with similar metabolic profiles have a close relationship on the evolutionary tree for their encoding genes? Second, whether the TR-encoding genes with a close evolutionary relationship have similar metabolic changes to the gene deletion? Answers to these questions could provide clues to possible function, evolution and origin of the TR genes. To seek answers to the questions, a gene families’ analysis was conducted using full protein sequences of 32 TR genes. The homology analysis showed that 32 TR genes were classified into several different families, indicating different origin during the evolution of TR genes (Figure 4 and Supplementary Table S3). Meanwhile, several pairs of TRs were found in the same gene family, suggesting possible events of gene duplication in recent evolutionary courses. Interestingly, while comparing the trees resulted from the hierarchical clustering analysis of metabolomic data and the gene families, we found that several pairs of TRs were clustered together in same gene family and the relationship tree generated using metabolomic profiles, suggesting possible functional conservation of the genes during genome evolution (Figure 4). At 48 h, one pair of TR genes, slr1489 (pchR) and sll1408 (pcrR), with similar function and evolution was identified; however, when metabolic profiles of 72 h were used, five pairs of TRs were identified. The pairs were slr0895 (prqR) and sll1286. slr1489 (pchR) and sll1408 (pcrR); sll1712 and sll1670 (hrcA); slr0115 (rpaA) and slr0449 (dnr); slr1871 (pirR) and slr1245. Only one pair of TRs, slr1489 (pchR) and sll1408 (pcrR), was identified when using metabolomic profiling data of both time points, suggesting a very conserved role that they might be playing. The difference between 48 and 72 h was probably due to the phase-dependent regulation of TRs in Synechocystis, which has been commonly reported in various microbes (Zhu et al., 2015).

Meanwhile, the results also showed a total of 17 pairs of TRs that were clustered together in the metabolomic trees but did not locate in the same gene family (Figure 4). Although experimental measurements of more metabolites are still necessary, the preliminary analysis pointed to the possibility of functional convergence of TR genes during the genome evolution, as the results showed that deletion of TR genes of different evolutionary origins caused similar metabolic responses in the mutant cells. For example, the slr0395 (ntcB) and slr0724 (sohA) genes were clustered together in the tree resulted from metabolic profiles of 48 h, deletion of these two TR genes caused CoA and GAP were down-regulation while ADP and ADP-glucose were up-regulated, although they belong to different gene families.

Experimental and computational data from genome-scale PPI analysis has contributed significantly to the understanding of the gene function (Marcotte et al., 1999; Ikeuchi and Tabata, 2001; Sato et al., 2007). In this study, an attempt was also made to apply PPI network of Synechocystis to determine possible gene targets of the five pairs of TRs identified above with similar clustering patterns in metabolomic profiles and matching gene family membership. As TRs of each pair have similar evolutionary and metabolomics patterns, it was expected that they might function through the same targets. Based on this hypothesis, we then implemented a strategy to first determine regulatory targets of each TR using PPI network analysis, and then identify the common target genes for every pair of TRs. The analysis allowed inference of possible regulatory gene targets for four out of five pairs of TRs, although no common target gene was identified for the slr1871 (pirR) and slr1245 pair (Supplementary Table S2).

In this study, 32 knockout mutants for putative TR-encoding genes of Synechocystis were constructed and comparatively analyzed via LC-MS-based metabolomics. Four mutants, sll1872 (lytR). slr0741 (phoU). slr0395 (ntcB), and slr1871 (pirR), showed differential growth patterns in BG11 medium when compared with the wild type. In the remaining TR mutants that did not show growth difference compared with the wild type, metabolomic profiling showed that they were clearly different at the metabolite level, suggesting significant functional diversity of TRs in Synechocystis. Finally, protein-protein interaction network analysis predicted possible regulatory targets of TRs.