In silico screening for candidate chassis strains of free fatty acid-producing cyanobacteria

Establishing properties that are favorable for cyanobacterial FFA cell factory

The common procedures used to enhance the biotechnological production of FFA include the introduction of heterologous pathways, as well as the modification of the candidate cell factory metabolism via deletion of genes or enhancing gene expression. However, genetic engineering was not based on the consideration of the collective effects of different criteria that characterize a good cyanobacterial cell factory for FFA production, even though experimental outcomes have shown that not all cyanobacteria are suitable producers [13, 14, 32]. Criteria that would potentially characterize the natural candidate cyanobacterial FFA cell factory include the presence of endogenous FA biosynthesis pathway enzymes [11, 34], as well as associated enzymes that have been modified and tested (through the insertion, overexpression, knockout or knockdown of protein-encoding genes) to increase FFA production in organisms such as algae, cyanobacteria, yeast, E. coli and diatoms [11–17, 32, 35–44]. Through a literature search, we identified 64 proteins that are relevant for FFA production. We further classified these 64 proteins into 49 OGs (Table 1, Additional file 1: Table S1), defined here as sets of proteins that are homologous with sufficient domains in common adequate to assume that they affect FA production similarly. To illustrate how these 49 OGs (into which 64 proteins are classified) affect FFA production, in Fig. 1 we show the link of the 49 OGs with the associated metabolic pathways and links to processes associated with energy, carbohydrate and lipid metabolism. Although these 64 proteins cannot be considered complete, they represent the majority of engineering considerations. Based on the results we obtained, it appears these proteins capture many of the relevant characteristics of the organism.

KEGG Orthology	Definition	Effects	Method	Organism	Ref.
rOGs
K00873  pyk	pyruvate kinase	Carbohydrate metabolism		cyan.	[11, 34]
K01007  pps	pyruvate, water dikinase
K00161  pdhA	pyruvate dehydrogenase E1 component alpha subunit
K00162  pdhB	pyruvate dehydrogenase E1 component beta subunit
K00627  pdhC	pyruvate dehydrogenase E2 component (dihydrolipoamide acetyltransferase)
K00382  phdD (ipdA)	dihydrolipoamide dehydrogenase
K00648  fabH	3-oxoacyl-[acyl-carrier-protein] synthase III	Lipid metabolism
K00645  fabD	[acyl-carrier-protein] S-malonyltransferase
K09458  fabF	3-oxoacyl-[acyl-carrier-protein] synthase II
K02372  fabZ	3-hydroxyacyl-[acyl-carrier-protein] dehydratase
K00208  fabI	enoyl-[acyl-carrier protein] reductase I
K01046  E3.1.1.3	triacylglycerol lipase	Increase chance of strain to secrete FA	secretion & extraction		[15]
pOGs
K01962  accA	acetyl-CoA carboxylase carboxyl transferase subunit alpha	Enhance FFA production (Increase supply of desired substrate)	secretion	cyan.	[12, 14, 34, 44]
K01963  accD	acetyl-CoA carboxylase carboxyl transferase subunit beta
K01961  accC	acetyl-CoA carboxylase, biotin carboxylase subunit
K02160  accB	acetyl-CoA carboxylase biotin carboxyl carrier protein
K00432  gpx	glutathione peroxidase	Reduce the toxic effect of FFA production and improve cell growth, physiology and FFA production	secretion	cyan.	[13]
K04564  SOD2	superoxide dismutase, Fe-Mn family
K06198  coiA	competence protein CoiA
K03782  katG	catalase-peroxidase
K03621  plsX	glycerol-3-phosphate acyltransferase PlsX	Lead to higher lipid levels		plant	[11, 35]
K08591  plsY	glycerol-3-phosphate acyltransferase PlsY
K00655  plsC	1-acyl-sn-glycerol-3-phosphate acyltransferase
virNOG10454 PDAT1	IQ-domain	Enhancing FA synthesis and diverting FA from membrane lipid to Triacylglycerol	accu.		[36]
virNOG19439  OLEO1	oleosin 1
K14457  MGAT2	2-acylglycerol O-acyltransferase 2	Enhance acyl-CoA-dependent triacylglycerol TAG			[39]
virNOG24576 LCIA	Anion transporter	Help regulate CO2 intake and increase biomass		algae	[32, 37]
virNOG22763  LCIB	Low-CO2 inducible protein
K00006  GPD1	glycerol-3-phosphate dehydrogenase (NAD+)	Increase glycerol and neutral lipid content (16- and 18-carbon monounsaturated FA significantly increased)		diatom	[38]
K01601  rbcL	ribulose-bisphosphate carboxylase large chain	Improve FFA production		cyan.	[14, 88]
K01602  rbcS	ribulose-bisphosphate carboxylase small chain
K01648  ACLY	ATP citrate (pro-S)-lyase	Enhance biofuel precursor production		yeast	[40]
K10804  tesA	acyl-CoA thioesterase I	Remove feedback inhibition and increase production of FFA	secretion	cyan.	[12, 14, 17, 34, 43, 44]
K10781  FATB	fatty acyl-ACP thioesterase B (Plant thioesterase)	Modify the chain length of FFAs for better fuel quality			[11, 12]
K10782  FATA	fatty acyl-ACP thioesterase A	Release FFA			[88]
K14075  PLRP2	pancreatic lipase-related protein 2	Degrade the membrane lipids into FFA with collapse of cell	extraction		[15]
nOGs
K01595  ppc	phosphoenolpyruvate carboxylase	Increase the lipid content		cyan.	[11]
K01897  aas(fadD)	long-chain acyl-CoA synthetase	Channel needed substrates for synthesis of FFA into divergent or reverse pathways and preventing degradation of desired product	secretion	cyan.	[11–14, 17, 34]
K00059  fabG	3-oxoacyl-[acyl-carrier protein] reductase	Divert energy into production of substantial by-products that would compete with production of FFA			[12]
K00626  E2.3.1.9	acetyl-CoA C-acetyltransferase
K11003  hlyD	hemolysin D	Enhance secretion of FFA by weakening cell walls
cyaNOG01264  (PBP2)	penicillin-binding protein	Enhance secretion of FFA by weakening peptidoglycan layer
K13788  pta	phosphate acetyltransferase	“Channel needed substrates for synthesis of FFA into divergent or reverse pathways and preventing degradation of desired product”
K13282  cphB	cyanophycinase	“Divert energy into production of substantial by-products that would compete with production of FFA”
K03802  cphA	cyanophycin synthetase
cyaNOG01069  porin protein	Carbohydrate-selective porin OprB	Enhanced extracellular FFA concentration			[13]
K13535  CLD1	cardiolipin-specific phospholipase	Increase lipid yields without affecting growth or biomass	accu.	diatom	[41]
K00030  IDH3	isocitrate dehydrogenase (NAD+)	Increase intracellular citrate level which enhance biofuel precursor production		yeast	[40]
K03603  fadR	GntR family transcriptional regulator, negative regulator for fad regulon and positive regulator of fabA	Fatty acid biosynthesis is feadback-inhibited at the transcriptional level by fadR		bacterium	[11, 42]

Abbreviations: rOGs required OGs, pOGs, OGs that positively impact FFA production, nOGs, OGs that negatively impact FFA production, FFA Free Fatty Acid, accu. Accumulation, cyan. Cyanobactia

Classification: nOG (based on reported knockout or knockdown) and pOGs (based on reported inserted or overexpressed) during genetic engineering experiments on that organism in order to secretion, extraction, or accumulation fatty acid

In total, we identified 13 OGs (based on reported knockout or knockdown experiments) whose presence in the organisms negatively impacts FFA production. These proteins we collectively named nOG (‘negative OG’; Additional file 1: Table S2). Acyl-ACP synthetase/long-chain-fatty-acid CoA ligase (AAS/FadD) is an example of one of the cyanobacterial proteins from this group. Kaczmarzyk and Fulda [45] demonstrated AAS is capable of incorporating exogenous FFA from the culture medium into membrane lipids, an opposite process that reduces FFA production. AAS is also responsible for recovering endogenous FFA released from membrane lipids. aas knockout mutants for Synechocystis sp. PCC 6803 and S. elongatus PCC 7942 (strain SE01) exhibited increased secretion of FFA into the culture medium compared to the wild-type strains [45]. The data suggests that the detected FFA is detached from membrane lipids, and also suggests that AAS plays a role in recycling the released FA, explaining why the presence of the aas gene negatively impacts the efficiency of the candidate cell factory.

Based on reported gene insertion and overexpression experiments, we also identified 24 OGs that contain proteins whose presence in the organisms positively impacts FFA production capability (named pOG; Additional file 1: Table S2). Thioesterase (TesA) is an example from this group. It was previously demonstrated that TesA cleaves the acyl-carrier-protein from the FA moiety, and in this manner increases FA biosynthesis in E. coli by reducing feedback inhibition [46]. Thus, Ruffing and Jones [17] cloned the E. coli-derived truncated thioesterase (‘tesA) and inserted it into the S. elongatus PCC 7942 genome along with the aas knockout, thereby generating a mutant strain SE02. SE02 produced a higher percentage of saturated FFA and a lower percentage of unsaturated FFA compared to the wild type [17]. Thus, the presence of ‘tesA positively impacted the efficiency of the biofuel production. The remaining 12 OGs identified are required for FA production, but are not included in pOG, and we named them rOG (‘required OGs’). The difference between these two groups is that rOGs are essential for FFA production, while pOGs can be considered as ‘enhancers’.

Based on these 49 OGs and their subgrouping to nOG, pOG and rOG, we derived criteria for assessment of suitability of an organism for FFA production (see Materials and Method section, subheading FFASC). In order to estimate an organism’s potential for FFA production, we used all of these derived criteria to generate an overall score that reflects FFA potential. For this purpose we developed FFASC. Our optimization process through which we estimated the optimized weights of the criteria used, is based on two species, Synechocystis sp. PCC 6803 and Arthrospira (Spirullina) platensis NISE-39. Thus, our estimated weights are skewed and not optimal. However, they still provide better qualitative ranking of species for FFA production potential than in the case when all weights are assumed to be equal (see Additional file 1: Table S10). These weights could be improved when more confirmed FFA-producing strains become available for this type of study.

Screening cyanobacterial proteomes by FFASC

Even though a limited number of cyanobacterial strains have been engineered as FFA/biofuel producers, several trends can be identified. Wild type Synechococcus sp. PCC 7002, Synechocystis PCC 6803 and Synechococcus elongatus PCC 7942 are reported to produce approximately 2.5 [14], 1.8 [12] and 0.3 [14] mg/L of FFA, respectively. However, these criteria are generally not sufficient to identify the putative chassis strains. Ruffing [14] has demonstrated that Synechococcus sp. PCC 7002 is a superior host strain compared to S. elongatus PCC 7942 regarding biomass growth rate, environment tolerance, FFA tolerance and production. The ‘tesA-expressing aas-deficient mutants’ of Synechococcus sp. PCC 7002, Synechocystis PCC 6803 and Synechococcus elongatus PCC 7942, showed an increase in FFA concentration of 40 [14], 83.6 [12] and 29.3 [14] mg/L, respectively, indicating that the increase in FFA concentration depends on the favorable traits in each organisms overall genetic make-up. An additional genetic manipulation, that is, the overexpression of Rubisco, in Synechococcus sp. PCC 7002 further increased the FFA concentration to 103 mg/L. To-date the strain with the most genetic manipulations is Synechocystis PCC 6803, which yields the highest FFA concentration of 197 mg/L. However, its genetic modifications include weakening of the cell wall layers that may affect survival capabilities under adverse conditions [12]. It was also demonstrated that while engineered S. elongatus PCC 7942 strains successfully produce and secrete FFA, these cells are compromised with a decrease in Chl-a content and photosynthetic yield, as well as changes in pigment localization that may be partially attributed to the unsaturated FFA being oxidized into toxic products [17]. Such cell physiology associated ramifications are not known for engineered Synechocystis sp. PCC 6803. However, engineered Synechocystis PCC 6803 were reported to mainly produce saturated FFA. These potential differences in the host metabolism suggest that Synechocystis sp. PCC 6803 may be a better chassis strain for FFA production than S. elongatus PCC 7942. Nonetheless, both Synechocystis PCC 6803 and S. elongatus PCC 7942 are fresh water strains. On the other hand, marine strain Synechococcus sp. PCC 7002 has been shown to endure salt concentrations up to 1.7M [47], making it an attractive target for large-scale production using marine water based media. Synechococcus sp. PCC 7002 may also be the superior chassis strain, compared to both Synechocystis sp. PCC 6803 and S. elongatus PCC 7942, owing to its short doubling time and remarkable light and temperature tolerance [14]. Additionally, Lyngbya sp. PCC 8106 was shown to produce less FFA/biodiesel than S. elongatus PCC 7942 [48], while A. platensis NIES.39 showed resistance to genetic manipulation [19, 49]. Thus, the positive reference chassis strains include Synechococcus sp. PCC 7002 and Synechocystis sp. PCC 6803 as they are easily genetically modified and show superior FFA production followed by S. elongatus PCC 7942. Thus, Lyngbya sp. PCC 8106 and A. platensis NIES.39 are considered in this study as negative reference hosts. Due to the limited number of candidate cyanobacterial FFA producers. Moreover, taking into account the reported outcomes for five cyanobacterial species included in our reference dataset, Synechococcus sp. PCC 7002 is expected to perform better than both Synechocystis sp. PCC 6803 and S. elongatus PCC 7942, followed by Lyngbya sp. PCC 8106 and A. platensis NIES.39.

The subsequent list of ranked cyanobacterial strains demonstrates that the positive reference strains rank above the negative reference strains. However, they are not the top ranked strains. The positive reference strains Synechococcus sp. PCC 7002, Synechocystis sp. PCC 6803 and S. elongatus PCC 7942, ranked at position 21, 23 and 37, respectively, while negative reference strains A. platensis NIES.39 and Lyngbya sp. PCC 8106 ranked at positions 101 and 123, respectively (Table 2). Thus, 36 cyanobacterial strains were ranked above the lowest ranked positive control reference strain at position 37, of which 20 strains (denoted as top ranked strains) ranked above all positive reference strains. All 20 top ranked strains are unicellular. We further observed that the reference strains were ranked as per experimental outcomes reported in the literature. Additionally, weights assigned to criteria after optimization show that 21 of the 49 criteria have the greatest impact on the score and thus the ranking of the strains for FFA production potential (Table 3). However, the criteria impact the score of every strain differently as this impact depends on the composition of the strain’s proteome. We point out that since we are interested in the organism’s natural potential to produce FA, we did not normalize the results for the genome size. We further provide heatmap visualization of the cyanobacteria screened for their potential as FFA producers against the 49 OGs (Fig. 2). The heatmap shows that the majority of the top ranked strains (above Synechococcus sp. PCC 7002) are placed in one major clade along with cyanobacterial positive reference strains, while the diatoms, used as an outgroup needed for hierarchical clustering, are placed in a clade of their own. Also, the negative reference strains do not mix with the clade that contain the top ranked strains, that is, the heatmap shows a clear separation between these clades. Moreover, the major clade that contains the top ranked strains generally has a higher number of pOGs (represented by the reddish shaded area) and lower numbers of nOGs (represented by the greenish shaded area), which contrasts with the clade in which negative reference strains are placed. Taken together, the clade with top ranked strains displays more favorable traits for FFA production based on the 49 OGs assessed.

Importance of features
Features	Weight
overexpression_K00432_Synpcc7942_1214	0.999999981
overexpression_K04564_Synpcc7942_0801	0.999999101
overexpression_K03782_Synpcc7942_1656	0.999998942
overexpression_K02160_accB	0.999998794
present_K00873_pykf	0.999998794
knockout_K11003_hemolysin	0.999998724
underexpression_K13535_Thaps3_264297	0.997856841
knockout_cyaNOG01069_porin	0.946931624
knockout_K01897_fadD	0.921718924
present_K09458_fabF	0.456133041
overexpression_K00006_GPDH	0.396694273
present_K00208_fabI	0.387646822
present_K00161_pdhA	0.314952182
present_K02372_fabZ	0.288150995
overexpression_virNOG24576_LCIA	0.228675187
present_K00648_fabH	0.17462096
present_K00627_odhB	0.168613677
present_K00645_fabD	0.160058392
insert_K01602_rbcS	0.150753918
present_K01046_lipase	0.14966023
overexpression_K06198_Synpcc7942_0437	0.119438174
present_K01007_pps	0.020105541
insert_K14075_gpl	0.013465425
overexpression_K00655_plsC	0.008613511
overexpression_virNOG10454_PDAT1	0.00833575
overexpression_K01963_accD	0.008246186
overexpression_K01961_accC	0.008089475
knockout_K00059_fabG	0.007999865
overexpression_K08591_plsY	0.007682015
overexpression_K01962_accA	0.007630664
insert_K10804_tesA	0.005907112
knockout_K00626_thi	0.004833629
knockout_cyaNOG01264_PBP2	0.004590024
knockout_K03802_slr2002	0.004303632
knockout_K03603_fadR	0.004102976
insert_K01601_rbcL	0.003963175
knockout_K00030_idh	0.003153309
knockout_K13788_pta	0.001763091
overexpression_K03621_plsX	0.001763091
overexpression_virNOG22763_LCIB	0.001763091
present_K00162_pdhB	0.001763091
present_K00382_phdD	0.001763091
underexpression_K01595_ppc	0.001763091
overexpression_virNOG19439_oleosins	0.001299491
knockout_K13282_slr2001	0.00115274
insert_K01648_acl	0.001045169
insert_K14457_DGTT2	0.001001378
insert_K10781_fatB	0.001000657
insert_K10782_fat1	0.001000152

A more in depth assessment of the weights assigned to the 49 OGs (see Table 3) revealed that the medium ranked group (with optimized weights in the range 0.12-0.46) contains mostly the core enzymes of the general fatty acid biosynthesis pathway. These core enzymes are necessary for any producer strain, and their presence cannot be expected to distinguish weak from strong producers. By contrast, the top ranked group (optimized weights in the range 0.92-0.99) contains two principal categories of enzymes: those ensuring precursor availability for biosynthesis of lipids and those involved in handling the oxidative stress associated to FFA synthesis. Belonging to the first category are acetyl-CoA carboxylase [12, 14], pyruvate kinase [11], and acyl-ACP synthetase/long-chain acyl-CoA synthetase [11]. These key enzymes have been validated as metabolic engineering targets for increasing the flux of lipid production [12], and it is not surprising that they have been ranked in the top group. Recently, it was shown that the production of FFAs in cyanobacteria entails the creation of high levels of reactive oxygen species (ROS) which causes oxidative stress, and ultimately loss of membrane integrity [13]. Several enzymes identified in the top group provide relief from oxidative stress and/or are related to membrane permeability: glutathione peroxidase, superoxide dismutase, catalase and porin. Under light, photosynthesis is known to induce the production of ROS which cause lipid peroxidation [50], and the activity of the above-mentioned enzymes can thus also ensure quality control of the produced lipids. A multifunctional lipase was also identified in the top group, coherent with the finding by [51] that stimulating lipid catabolism is required to balance lipid accumulation with efficient growth. The composition of the top group therefore reflects the requirement for the producing cell to handle the flux control points (precursors, lipid accumulation versus biomass accumulation) and to possess enzymes enhancing stress tolerance related to lipid accumulation (ROS/membrane stress tolerance). The weight values obtained during the optimization procedure thus reflect the importance of these two types of key markers for affecting the strain’s potential as cell factories that can be expected to reach a high titer of lipids.

Comparison between FFASC and Model SEED

Since, Model SEED [52] automatically produces annotations and draft genome-scale metabolic models, we used it here to compare its results with the proposed FFASC approach using the EC numbers corresponding to the 49 OGs that affect FFA production. We found that 41 of the 49 OGs in FFASC can be used for a comparison with Model SEED, as it only focuses on enzymes required for metabolic model reconstruction. Thus, the eight OGs omitted from this analysis include one enzyme that does not have a defined EC number such as EC 3.1.1.-, while other OGs are proteins that do not function as enzymes. For the 41 OGs (Fig. 3), we found Model SEED and FFASC have 28 identical OG hits (68%) for all 25 cyanobacterial strains screened (these are the 20 top-ranked cyanobacterial strains and the five control reference strains). FFASC showed the presence of nine OG hits (22%) that were not present in Model SEED for some species. Similarly, Model SEED showed the presence of four OGs (10%) that were not found to be present using FFASC.

To analyze this data, we tabulated the engineered genes in model organisms Synechocystis sp. PCC 6803, Synechococcus sp. PCC 7002 and S. elongatus PCC 7942, to show the set of genes known to be present in these organisms (see Additional file 1: Table S5). Liu et al. [12] made six successive generations of genetic modifications for Synechocystis sp. PCC 6803, these modifications include the knockout of slr2001 and slr2002, which encode the cyanophycin synthetases [53]. This shows that slr2001 and slr2002 are known to be present in Synechocystis sp. PCC 6803, and is reported as present by FFASC, but absent in Model SEED. We further verified that RAST [54] correctly annotated both slr2001 and slr2002 in the Synechocystis sp. PCC 6803 genome. However, it was omitted from Model SEED, due to the lack of gene-protein-reaction (GPR) association required for incorporation into SEED models. For the four enzymes missing from FFASC, another modification made by Liu et al. include the knockout of the slr1710 (PBP2) gene responsible for peptidoglycan layer assembly [55]. This shows once again that slr1710 is known to be present in Synechocystis sp. PCC 6803, and is correctly found by both Model SEED and FFASC. However, we found that Model SEED identified slr1710 in 22 additional cyanobacterial strains, whereas FFASC only identified slr1710 in 11 additional cyanobacteria screened. We found FFASC filtered out the other slr1710 hits as a consequence of the stringent protein-domain condition applied to increase the accuracy underlying FFASC predictions, that is, only homologous protein sequences that have all domains of the associated protein from the group of 64 proteins were recorded as OG hits. Moreover, all the core enzymes of the general fatty acid biosynthesis pathway were identified using FFASC, whereas Model SEED did not identify FabZ due to the lack of GPR association required for incorporation into SEED models. Here, the differences between Model SEED and FFASC are a consequence of: 1/Model SEED is a generic method in which all pathways are treated equally, whereas FFASC is specialized and focuses on FFA production and is built based on proteins known to either positively or negatively affect FFA production; 2/Model SEED provides the presence or absence of the enzymes, whereas FFASC takes the copy number into account when assessing potential for FFA production; and 3/FFASC include all proteins (not just enzymes) that directly or indirectly affect FFA production. Taken together, FFASC is more refined in assessing the “natural” cyanobacterial strains potential for FFA production, whereas Model SEED was developed for a more generic purpose.

Additional in silico support for estimated FFA production potential of cyanobacteria

To provide additional support that the predictions obtained by FFASC are reasonable, we used K-means clustering [56] based on the same 49 criteria. To cluster the 128 target species into k clusters, where distance of species within a single cluster is minimized and distance between clusters or cluster centers is maximized, a value for k has to be set in away that reflects the natural groupings. That is, if k is too small, the clustering algorithms will reduce the total number of groups to the specified value of k, which forces some natural clusters to combine, thereby producing artificial fusions [57]. Likewise, if the value of k is too large, natural clusters will start dividing in an artificial way, to match the specified k value.

To determine the appropriate number of clusters, we take into account that diatoms are eukaryotes and thus act as a type of outlier. When they fall into the same cluster this would indicate the point at which the artificial grouping is omitted [57]. Thus, the clustering will be considered good when diatoms fall into a separate cluster. The number of clusters where diatoms start to group together is k = 6 and k = 7, the point at which diatoms start to separate is when the number of clusters is k = 8. Additionally, using an average silhouette width as the measure of ‘natural’ clustering [57], we found that when considering k = 6, 7 or 8, the highest average silhouette width of 0.41 (Fig. 4) was associated with k = 6. To further verify the appropriate number of clusters, we also calculated the Calinski-Harabasz (CH) index for k = 6 (67.43), k = 7 (56.91) and k = 8 (61.89) (starting from the point when diatoms cluster together without cyanobacteria, to the point where the diatoms start to separate into different clusters). CH index results verify that k = 6 is the appropriate cluster number. A visual illustration of the case k = 6 (Fig. 5) shows that cluster 3 is the most distant from the other clusters. This cluster includes the 3 diatoms alone as the outliers, while the negative reference host Lyngbya sp. PCC 8106 and A. platensis NIES.39 were placed in cluster 5. Top ranked strains, above Synechococcus sp. PCC 7002, were all placed in cluster 6. Moreover, all positive reference chassis strains; Synechococcus sp. PCC 7002, Synechocystis sp. PCC 6803 and S. elongatus PCC 7942 were grouped together in cluster 4. Additionally, all strains that ranked below Synechococcus sp. PCC 7002 but above S. elongatus PCC 7942, were either placed in cluster 6 or 4. The placement of cluster 4 was closest to cluster 6; these clusters slightly overlap one another, but are separate from the other clusters. This indicates that even though K-means clustering does not rank strains, it is still able to discern the potential FFA producers identified with FFASC by clustering them primarily in cluster 6 based on the OG criteria.

Phylogenetic relationships of cyanobacteria

We explored phylogenetic groupings of 124 cyanobacterial strains used in this study. We found that several of our top ranked candidate cyanobacterial strains are grouped together based on their 16S rRNA. Some exceptions include two Thermosynechococcus sp., two Synechococcus sp. JA* and Candidatus Atelocyanobacterium thalassa (isolate ALOHA) (Fig. 6).

This result is supported by literature, since the top ranked cyanobacterial strains primarily include Prochlorococcus (order Prochlorales) and marine Synechococcus (order Chroococcales), which are reported to have diverged from common ancestry [58]. Following the divergence, the Prochlorococcus genome is further thought to have ‘streamlined’ [59], thus, the genome size of Synechococcus and other cyanobacteria is larger than Prochlorococcus genome sizes [60]. Another key feature that differentiates Prochlorococcus from Synechococcus is their divergent light-harvesting strategies [61]: Synechococcus uses the phycobilisome as their light-harvesting antenna that are not found in Prochlorococcus. These phycobilisome antenna systems are used by Synechococcus to adjust to changes in temperature, likely contributing to its greater geographical occupancy range [62, 63]. Instead, the Prochlorococcus main light-harvesting antenna complex is made up of divinyl chlorophyll a and b, prochlorophyte chlorophyll-binding protein (Pcb), as well as accessory pigment [60, 64]. Collectively, these pigments increase blue light absorption that is the dominant wavelength in deep waters, restricting Prochlorococcus to warmer, oligotrophic oceans [65]. Since Prochlorococcus is reported to be a leading example of a naturally 'streamlined' genome [59, 66], this suggests that these genomes may require less engineering to efficiently produce high yields of FFA. Moreover, Prochlorococcus can be inexpensively cultivated using seawater [67].

Reference strains of the order Chroococcales, including Synechococcus PCC 7002, Synechocystis PCC 6803 and S. elongatus PCC 7942, were engineered, and demonstrate the production and secretion of FFA, which provides proof-of-concept. However, none of the predicted top ranked strains of the order Chroococcales has been shown to produce FFA. Nonetheless, Synechococcus UTEX 2973 (which was not included in this analyses because its genome sequence was not available at the time of this study), has recently been reported to be a fast growing chassis strain for biosynthesis using light and carbon dioxide, growing two times faster than S. elongatus PCC 7942 [68]. This finding demonstrates that there are possibly more suitable chassis strains that have not been investigated. Moreover, the Chisholm group [69] have reported that the Prochlorococcus strain MIT9313 produces lipid-containing vesicles that are released into the surrounding seawater. These released lipid-containing vesicles maybe collected without disturbing the growth of the Prochlorococcus, as opposed to other cyanobacteria or algae that require destroying one batch of cells and starting with a new batch, to retrieve lipids. FFASC ranked the Prochlorococcus strain MIT9313 at position 41, suggesting that if the MIT9313 mechanism is a Prochlorococcal trait, there are several other possible vesicle-releasing Prochlorococcus strains that may be a better chassis for FFA production. Moreover, the fact that the candidate chassis strains are clustered primarily in orders Synechococcus and Prochlorococcus, is a welcomed surprise that could constitute an additional criterion for positive prediction.

Compilation of protein groups that characterize FFA production and excretion

The PubMed database was queried using the query: "biofuel production" OR "free fatty acid production" on 2015/06/30, resulting in 1392 PubMed abstracts retrieved. We conducted a literature search to a compile list of proteins relevant for FFA production from organisms that have been genetically engineered for FFA/biofuel production, as well as proteins required for fatty acid synthesis. In total, we identified 64 such proteins in various organisms including Escherichia coli, cyanobacteria, algae, diatoms, plants, and yeast (Additional file 1: Table S1 and S2). These 64 proteins can be classified into 49 OGs, with 43 from KEGG and six from the eggNOG (evolutionary genealogy of genes: Non-supervised Orthologous Groups) [70] database. The 43 KEGG orthology KO identifiers were associated to these proteins using the KOALA (KEGG Orthology And Links Annotation) [71] tool. For the remaining six OGs with no associated KO identifiers, we used eggNOG (Additional file 1: Table S2) to associate OGs to the remaining proteins from the group of 64. All protein sequences included in the 49 OGs were extracted from the UniProt [72] database.

Based on the effects that the presence or absence of relevant genes have, a set of rules is derived to quantify these effects (see Criteria Generation section).

Compilation of control and target datasets

Sequence homology and domain search

Protein sequences included in the 49 OGs were mapped to 125 cyanobacterial proteomes using a protein homology search, with the local installation of BLASTp [80, 81], and with an e-value threshold of 0.0001.

We identified 81 conserved protein domain families in the 64 originally identified proteins, using the Pfam database and HMMER[82] with the cut-off gathering threshold (Additional file 1: Table S6). The hidden Markov model (HMM) profiles of these domain families were retrieved from the Pfam database.

The homologous protein sequences identified in the 125 cyanobacterial were further screened with the 81 HMM profiles using a locally installed HMMER [83] program with the trusted cutoff score as a threshold. In the analysis, only homologous protein sequences that have all domains of the associated protein from the group of 64 proteins are used (refer to Fig. 7).

Criteria generation

In order to provide an integral score of the potential for a species to produce and excrete FFA, we need to quantify the effects of presence or absence of genes that encode for relevant proteins. In our case these will be proteins from different OGs. We consider this quantification as criteria, and we derive one criterion for each OG.

The number of BLASTp hits of all proteins from an OG to the proteome of a species represents an OG hit number (hitN). hitNs are used to define criterion for the OGs. In determining hitNs, only proteins matched by BLASTp that have all domains of the source protein were used. One can conveniently describe species and OGs in terms of hitNs as follows. Suppose that n is the number of species and m is the number of OGs. We can create an n × m matrix C. In our case C is 125 × 49 (see Additional file 1: Table S7). The element (i,j) of C represents hitN of j-th OG in i-th species.

The quantification rules are defined as follows. Proteins from nOGs receive the values equal to “–hitN” that correspond to the considered species and the OG. Proteins from pOGs receive the values of “hitN” that correspond to the considered species and the OG. If, however, a pOG has “hitN = 0”, then we assign to it a value of “-1” as a penalty. Proteins from rOGs receive the values of “hitN” that correspond to the species and the OG (Fig. 8).

where c(i,j) is an element of C. While this is not the only possible way to calculate this score we find it simple and suitable. Note that in (1) we assume that all criteria have the same weight equal to 1.

FFASC method

1. 1.

   Ranking Algorithm (Algorithm 1)

    

In order to determine scores for each of the species so as to be able to rank them, we will determine the C matrix and use it as the input to the algorithm. This algorithm evaluates each of the considered species and generates scores according to (1). Then, the species are ranked, with the higher score being better. The top rank is 1. In this manner we are able to rank the considered species for their FFA production potential based on the scores derived. A pseudo code for the algorithm (Ranking algorithm) is presented in Fig. 9.

1. 2.

   Optimization

    

In Algorithm 1 we assume that all criteria considered have the same level of influence to the potential of an organism for FFA production as expressed through Equation (1). However, it is reasonable to expect that different criteria have different levels of effects and thus they should have different weights. Because we have no data to determine precisely what values of these weights should be, we used an optimization approach in order to estimate suitable values of these weights. The general ‘constraint’ is that good producers of FA should be ranked higher and well separated from the poor ones. Thus, for the optimization process we selected a positive reference strain, Synechocystis sp. PCC 6803 and a negative reference cyanobacteria strain, A. platensis NIES.39. The goal of optimization was to make the score difference between these two selected species as big as possible, while having the positive reference strain ranked above the negative reference strain. Optimization was preformed using the pattern search solver (PSS) of the global optimization toolbox in MATLAB. For the PSS, a generalized pattern search algorithm was used with default values. The optimized solutions for the weights found by the optimizer were between 0.001 and 1, where p (ranking effect coefficient) is equal to 0.010241 at 1744 iterations, with the objective function value at the solution equal to 0.0232 (convergence level). The proposed objective function to achieve our goal is based on maximizing the difference in scores for the two species used; Synechocystis sp. PCC 6803 and A. platensis NIES.39 as defined below:

$$ \begin{array}{l}\underset{w}{ \max\ }{w}^T\ast \kern0.5em \left|{x}_1-{x}_2\right|+p \ast \kern0.5em rank\kern1.75em where\kern1.25em 1\ge {w}_j\ge 0.001,\ \\ {}\kern16em {\displaystyle \sum }{w}_j=12, \\ {}\kern16em {w}^T\kern0.5em \ast \kern0.5em \left|{x}_1-{x}_2\right|>\kern0.75em 0.001\end{array} $$

Here, x ₁ and x ₂ are data vectors describing Synechocystis sp. PCC 6803 and A. platensis NIES.39, respectively, obtained as rows of C; T denotes the transposition; |()| denotes the absolute value of (); w is a weight vector with values indicating the contribution of features as suggested by PSS; p is a coefficient to introduce a ranking effect on the optimization; rank is the difference in ranking between the Synechocystis sp. PCC 6803 and A. platensis NIES.39. In this optimization, an optimized set of weights are bounded and constrained as described above. Finally, having optimized the weights, we ranked 125 cyanobacteria, with the scores determined as

$$ score\kern0.5em =\kern0.75em {w}^Tx $$

(2)

Here, w is a column vector of dimension 49. Note that this procedure can be applied to the newly sequenced cyanobacteria species (or other species) added to the set we considered. The pseudocode of Algorithm 2 that describes ranking based on score determined by (2) is presented in Fig. 10.

Based on these optimized weights of different criteria, we propose a list of chassis candidate cyanobacteria strains, where the final ranking reflects the potential of the chassis strain to produce FFA (Additional file 1: Table S8).

Heatmap generation

The matrix C was modified following the MATLAB syntax to

C+(-0.5 + rand(size(C))*10^-10)

by adding a small level of noise to avoid numerical problems with singular matrices.

Generating data for comparison used in FFASC and Model SEED

The EC numbers corresponding to the 49 OGs were used for comparison with Model SEED. In addition, we submitted 25 cyanobacteria (which include the 20 top-ranked cyanobacterial strains by FFASC and the five control reference strains) to the Model SEED resource (using default values) and obtained the SEED metabolic models and corresponding genome annotations. Similarly, we had binary (presence/absence) output from our FFASC method. We compared the identified EC numbers of 41 OGs in both models and generated the comparison data for Model SEED and FFASC with binary values (0/1) (Additional file 1: Table S5).

We subtracted data for Model SEED from data for FFASC row-wise and obtained values ranging from -25 to 25, where values less than zero indicate the fulfillment of criteria in some strains as required by FFASC only, while values more than zero indicates the fulfillment of criteria in some strains as required by Model SEED only, while zero indicates that the same criteria were required by both FFASC and Model SEED.

Phylogenetic analyses

In order to see if the ranking obtained as described above reflects any phylogenetic similarities, we performed phylogenetic analyses of cyanobacteria. We used 16S rRNA sequences for the 124 cyanobacterial strains retrieved from INDIGO [79]. Synechococcus sp. CB 0101 was not included in this analysis as its 16S rRNA was not available. We also included 16S rRNA of the outgroup (Chlorobium tepidumdum, Rhodobacter sphaeroides and Chloroflexus aurantiacus). The 16S rRNA sequences for the 124 strains and outgroup were aligned using MAFFT (Multiple Alignment using Fast Fourier Transform) [84] with default parameters on the T-REX Web Server [85]. A maximum likelihood tree [86] was then generated based on the aligned 16S rRNA sequences using RAxML (Randomized Axelerated Maximum Likelihood), with default parameters and 1000 bootstrap runs for the GTRCAT substitution model [87]. The maximum likelihood tree was visualized using FigTree [88] and edited to improve visualization using Inkscape 0.91 [89].

K-means clustering

To further substantiate the results obtained by applying FFASC, K-means clustering was preformed on the 125 species using all 49 OGs. The K-means procedure in the Package ‘stats’ of R (R 3.1.2) [56] was used. To determine the proper number of k clusters, we established 1/the point at which artificial fusions are omitted, that is, when diatoms fall into a separate cluster (determined to be where k = 6) and 2/the point at which the natural clusters are divided in an artificial way, that is, when diatoms start to separate into individual clusters (determined as k = 8). Thus, based on the properties of the dataset, natural clustering was found to range from cluster 6 to 8 (Additional file 1: Table S9). Further analysis was restricted to natural clusters 6 to 8. To determine the optimal number of k clusters from this range, we used the largest average silhouette width as the measure of ‘natural’ clustering and calculating the CH index.