Genome-wide survey of the phosphofructokinase family in cassava and functional characterization in response to oxygen-deficient stress

Background Glycolytic pathway is common in all plant organs, especially in oxygen-deficient tissues. Phosphofructokinase (PFK) is a rate-limiting enzyme in the glycolytic pathway and catalyses the phosphorylation of fructose-6-phosphate to fructose-1,6-bisphosphate. Cassava (M. esculenta) root is a huge storage organ with low amount of oxygen. However, less is known about the functions of PFK from M. esculenta (MePFK). We conducted a systematic analysis of MePFK genes to explore the function of the MePFK gene family under hypoxic stress. Results We identified 13 MePFK genes and characterised their sequence structure. The phylogenetic tree divided the 13 genes into two groups: nine were MePFKs and four were pyrophosphate-fructose-6-phosphate phosphotransferase (MePFPs). We confirmed by green fluorescent protein fusion protein expression that MePFK03 and MePFPA1 were localised in the chloroplast and cytoplasm, respectively. The expression profiles of the 13 MePFKs detected by quantitative reverse transcription polymerase chain reaction revealed that MePFK02, MePFK03, MePFPA1, MePFPB1 displayed higher expression in leaves, root and flower. The expression of MePFK03, MePFPA1 and MePFPB1 in tuber root increased gradually with plant growth. We confirmed that hypoxia occurred in the cassava root, and the concentration of oxygen was sharply decreasing from the outside to the inside root. The expression of MePFK03, MePFPA1 and MePFPB1 decreased with the decrease in the oxygen concentration in cassava root. Waterlogging stress treatment showed that the transcript level of PPi-dependent MePFP and MeSuSy were up-regulated remarkably and PPi-dependent glycolysis bypass was promoted. Conclusion A systematic survey of phylogenetic relation, molecular characterisation, chromosomal and subcellular localisation and cis-element prediction of MePFKs were performed in cassava. The expression profiles of MePFKs in different development stages, organs and under waterlogging stress showed that MePFPA1 plays an important role during the growth and development of cassava. Combined with the transcriptional level of MeSuSy, we found that pyrophosphate (PPi)-dependent glycolysis bypass was promoted when cassava was under waterlogging stress. The results would provide insights for further studying the function of MePFKs under hypoxic stress. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-03139-7.


Background
Hypoxia is common in plant tissues. Oxygen falls to very low level in metabolically active or large sink tissues. Low internal oxygen levels have been found in many plants, such as growing tuber [1], developing seeds [2], fruits [3], root [4] and phloem tissue [5]. Low oxygen concentration in the phloem would interfere with sugar retention [6]. The size of Arabidopsis seed shows a notable reduction when the external oxygen concentration drops below 15%, and seed production is remarkably inhibited when the oxygen concentration drops to 2% [7].
Glycolysis is an important metabolic pathway in all organisms. Glycolysis provides energy for cell metabolism and help cells adapt to abiotic stresses, such as hypoxia, cold and drought. In classical plant glycolytic pathway, phosphofructokinase (PFK) is the main rate-limiting enzyme and regulatory point. PFK catalyses the phosphorylation of D-fructose-6-phosphate (F-6-P) to D-fructose-1,6-bisphophate (F-1,6-BP). Plants have two forms of PFK based on phosphoryl donors, namely, the ATPdependent PFK (EC 2.7.1.11) and the pyrophosphate (PPi)-dependent pyrophosphate-fructose-6-phosphate phosphotransferase (PFP, EC 1.7.1.90) [8]. PFK catalyses the irreversible phosphorylation of F-6-P to F-1,6-BP, and PFP catalyses the reversible phosphorylation of F-6-P to F-1,6-BP. PFK remains very complex in higher plants. Latzko and Kelly found the PFK isoenzymes in spinach as early as 1977 [9]. Since then, PFK from various plant sources have been studied. PFKs are found in plastids and cytosol in many plants, such as tomato fruit [10], R. communis seeds [11] and potato [12]. Winkler firstly purified PFK protein from spinach leaves and identified the existing PFK and PFP in the spinach genomic database [8]. PFK from growing potato contains four polypeptides with molecular weight between 46,300 and 53,300. PFK's isolation during the purification process is unstable; thus, study on PFK biochemistry is relatively few. Plant PFP has been characterised in much detail in many plants, such as S. tuberosum [13], R. communis [14] and S. officinarum [15]. PFP is common in all plant tissues, and its enzyme activity is usually equal to or greater than that of PFK [16]. PFP consists of two different subunits, namely, alpha and beta, which form a heterotetramer [12]. In potato, PFP is composed of two polypeptides with molecular weight between 58,000 and 55,700.
PFP is a self-adaptive enzyme, and its activity and subunit composition depend on environmental, developmental, tissue-specific and species conditions. PFP has multiple functions, including glycolysis, gluconeogenesis, regulation of PPi concentration in sucrose metabolism and plant adaptability to environmental changes. In Arabidopsis plants with PFP RNA interference, PFP activity is reduced and the growth of plants is severely inhibited [17]. The result shows that AtPFP is involved in carbohydrate metabolism. The transcript levels of several PFKs and PFPs from O. sativa (OsPFKs and OsPFPs) roots are increased during anoxia [18]. This phenomenon shows that PFP may be involved in glycolysis metabolism to help plants adapt to poor environment, such as drought, cold or low oxygen stress. In antisense potato plants of PFP, PFP expression is decreased remarkably with more than 90% reduction in PFP activity in growing tuber and stored tubers, but no visible phenotype is observed [19]. The decrease in PFP expression resulted in the decreased concentrations of phosphoenolpyruvate and glycerate-3-phosphate. Fructose-6-phosphate,2-kinase was stimulated to increase fructose-2,6-bisphosphate concentration, and ATP-dependent PFK was stimulated to compensate the decrease in PFP protein. The activity is increased remarkably in germinated coleoptiles and the suspension cell of rice under hypoxic stress, but the activity of PFK does not change [20]. Although PFP has received much attention, the function of PFP has not been well elaborated.
Progress in sequencing technology has enabled the whole-genome sequencing of many plant species. PFK gene families have been characterised in some plants, such as Arabidopsis [21], spinach [8], Saccharum [22], rice [20] and Chinese white pear [23]. Eleven members of the PFK gene family are present in A. thaliana: four members belong to AtPFP, and seven belong to AtPFK. Fifteen PFK genes are present in rice; five belong to OsPFP, and 10 belong to OsPFK. Fourteen members are present in white pear (Pyrus bretschneideri) PFK family: 10 belong to PbPFK, and four belong to PbPFP. Overall, the PFK family has more than 10 members, and PFK is the centre in plant growth and development and has diverse functions.
Cassava (M. esculenta Crantz) is an important food crop apart from rice and maize in Africa and Asia. Cassava has tuber root with starch that provides dietary carbohydrate and is used to produce industrial starch and bioethanol [24]. Tuber root, which is a huge storage organ, is predicted to be hypoxic; thus, glycolysis will be prevalent in cassava tuber root, especially the pyrophosphate-dependent glycolysis pathway, which will save energy for the survival and metabolism of cassava. Although the PFK family has received much attention, the PFK family in cassava has not been reported yet. In this study, 13 PFK genes were identified from the cassava genome and divided into two groups on the basis of the conserved PFK domain. Their sequence structures and subcellular localisation were investigated. The transcript levels of the PFK family under normal condition and hypoxic stress were explored. Overall, this study identified several tissue-specific and hypoxic stressresponsive PFK genes, and these genes will be helpful for the study of hypoxic stress in cassava plant.

Results
Identification and phylogenetic analysis of MePFK genes in cassava Thirteen MePFK genes in the M. esculenta genome were identified using the online BLAST programme of JGI cassava genome data using the known AtPFK gene as reference. All of these genes contain conserved domain PF00365, which is the basic characteristic of the PFK family. A neighbour-joining (NJ) phylogenetic tree was drawn based on the multiple alignments of the MePFK amino acid sequences and other PFK sequences from Arabidopsis, rice, castorbeen, tomato and potato to investigate the evolutionary relationships between MePFK protein and other PFKs from other species (Fig. 1). Thirteen MePFK proteins were divided into two groups, namely, the MePFK and MePFP subfamilies. The MePFK subfamily was divided into three subgroups, namely, PFK-A, PFK-B and PFK-C. The MePFP subfamily was divided into two subgroups, namely, PFP-α and PFP-β. The 13 predicted MePFK proteins ranged from 318 amino acids (MePFPB2) to 617 amino acids (MePFKA1) ( Table 1). The length of proteins varied distinctly; thus, different PFKs have different biological functions.
Exon-intron structure and motifs of MePFK were highly conserved We compared the exon-intron organisation in the coding sequences (CDSs) of MePFK genes to obtain a further insight into the structural diversity of the PFK genes. As shown in Fig. 2A, the members of MePFK with closely genetic relationship showed similar exon-intron structure. Most of the MePFK genes in cassava have more than 10 exons, except for MePFK08, MeFPK09 and MePFPB2. MePFK08 and MePFK09 belong to PFK-B and only contain four and two exons, respectively. Fig. 1 Phylogenetic relationship of PFK genes from cassava and four other species. The NJ tree was constructed with PFK proteins from cassava, Arabidopsis, potato, rice and castor bean using the Muscle programme of MEGA7 software with 1000 bootstraps. PFK protein was classified into two groups and five subgroups. The five subgroups are indicated with different background colours MePFK08 and MePFK09 have longer introns and exons than the other members. This special structure may be endowed with a special function. A total of 13-14 exons were found in the PFK-A and PFK-C subgroups. In the PFP subfamily, 19 exons were observed in MePFPA1 and MePFPA2, 16 exons in MePFPB1 and 8 exons in MePFPB2. The length of MePFPB1 transcript was 2209 bp, and that of MePFPB2 transcript was 2209 bp. By comparing the nucleotide sequences of MePFPB1 and MePFPB2, it was found that they had high homology (89%, Fig. S1). The homologous region is located at base 153 to base 1111 of 5′ in MePFPB1. We inferred that MePFPB2 is a mutation of MePFPB1, which stops transcription in advance.
We also predicted the conserved domains of all cassava MePFK protein sequences.   reflected that the PFK and PFP subfamilies have similar functions, and each has its own division of labour.

Chromosomal and subcellular localisation of PFK in cassava
The genomic distribution of PFK genes on the chromosomes of cassava was identified. A total of 13 MePFK genes were distributed throughout 10 of 18 chromosomes. Most of the 10 chromosomes contained one PFK gene, with chromosomes 9 and 16 containing two and three genes, respectively (Fig. 3).
According to the online subcellular localisation software, five of the seven MePFKs were predicted to be localised at the chloroplast, two were predicted to be localised at the cytoplasm, and four PFPs were predicted to be localised at the cytoplasm (Table 1). Two PFKs (MePFK03 and MePFPA1), which belong to different groups and predicted to be located in the chloroplast and cytoplasm, respectively, were selected to construct MePFK03-GFP and MePFPA1-GFP fusion proteins to confirm the aforementioned result. Fluorescent signal results showed that MePFK03 was localised in the chloroplast and MePFPB1 was localised in the cytoplasm (Fig. 4). This finding is consistent with the predicted results.

Cis-element prediction of cassava MePFK gene promoters
The cis-element within these promoters were identified using the online software new PLACE to better understand the functions of the MePFK genes. In addition to some common cis-elements, we found some special elements in MePFKs (Table 2). Amongst these elements, 17 elements were very typical, amongst which six belonged to organ-specific expression, seven were hormonerelated, and the remainder were environment related. As for the organ-specific cis-elements, all the 13 PFK promoters contained many mesophyll DE expression modules; the maximum was 26, the least was 2, and most of them had more than 15. OSEROOTNODULE, which is a motif of specific activated elements in root modules, was widely distributed in 12 promoters; the number was up to 74 in MePFK07 promoter but 0 in MePFK08 promoter. Many ROOTMOTIFTAPOX (rootspecific expression element) were distributed widely in 12 promoters but was 0 in the MePFK06 promoter. Thus, most of were expressed in cassava root. POL-LEN1LELAT52 (pollen-specific activation) was also widely distributed in 12 promoters but 0 in MePFPB2. A total of 28 RY-ELEMENT (specific in seed storage protein genes) were found in MePFK08 promoter but few or none in the other promoters. Therefore, MePFK08 had protein-related function in cassava. The cis-element prediction of these promoters indicated that MePFK was distributed in leaf, root and flower. Thus, the functions of the PFK lie in the leaf, root and flower. LTRE1HVBLT49 (low-temperature-responsive element) was distributed in 7 of 13 promoters. ANAERO1CON-SENSUS/ANAERO3CONSENSUS (anaerobic-responsive element) was widely distributed in all 13 promoters, which reflected that MePFK was involved in hypoxic stress.

Preliminary investigation on the oxygen concentration in cassava root
According to existing literature [1][2][3][4][5], the active metabolism tissues of plant storage organs are often in the state of low oxygen. Therefore, cassava tuber root, which is a large storage organ, is supposed to be in low oxygen state. According to the transverse section of the root tuber, the middle segment of the tuber root in 1-yearold roots of different cassava varieties (Arg7 and SC124) was divided into three regions (Fig. 5A). Region 1 refers to the periderm and includes the sclerenchyma, parenchyma and phloem; region 2 refers to the parenchyma;  Table 3. For Arg 7, the oxygen concentration under the periderm was 9-11.2%, and the concentration decreased to 2.58-3.96% in the centre. For SC124, the oxygen concentration under the periderm was 6.95-8.66%, and the concentration decreased to 3.7-5.2% in the centre, which showed a typical oxygen level profile throughout a growing tuber root. Overall, the oxygen concentration in root tuber decreased sharply from the outside to the inside. Compared with the oxygen concentration in the air (21%), hypoxia occurred in the inner cassava root.

Expression profiles of MePFK in cassava
The expression level of MePFK genes from different tissues, including leaves, petiole, stem, tuber root cortex, tuber root stele, fibre root and flower, were identified in SC124 to find some clues on the functions of PFK during the growth and development of cassava. As shown in Fig. 6A, MePFPA1 and MePFPB1 showed higher expression in all tissues than other genes. MePFPA1 and MePFPB1 were highly expressed in flower, leaf and fibrous root, which are metabolically active tissues. These genes also showed moderate expression in tuber root cortex and stele. MePFKs showed lower expression in all cassava tissues than MePFPs. Amongst MePFKs, MePFK02 and MePFK03 showed visible expression. The expression levels of PFK genes at three development periods of tuber root (90, 150 and 240 days after planting [DAP]) were also studied. As shown in Fig. 6B, MePFK02, MePFK04, MePFK06, MePFK08, MePFK09, MePFPA2 and MePFPB2 were almost not expressed in the root block. Thus, we only showed the higher expressed members of MePFKs. The expression levels of MePFPA1, MePFPB1 and MePFK03 increased gradually with the development of root tubers. This finding reflects that the functions of the three members were related to the development of tuber root.
MePFPA1, MePFPB1 and MePFK03 had relatively high expression in cassava root. Thus, they were selected to identify the expression pattern in cassava roots from different depths (Fig. 5B). The results showed that    Number of every selected cis-element in the MePFKs gene promoters. The yellow part indicated Tissue specific-relative elements; the orange part indicated hormone-relative elements; and the blue part indicated stress-relative elements members, and the expression of MeSuSy2, MeSuSy5, MeSuSy6 and MeSuSy7 is very low [25]. In this study, the expression profile in Fig. 7A indicates that the expression of MeSuSy3, MeSuSy4 and MeSuSy6 in roots was remarkably increased under waterlogging stress compared with normal control. Among all the MeSuSy genes, only MeSuSy5 was highly expressed in leaves. The expression of MeSuSy5 was down-regulated in cassava leaves under waterlogging. The expression profile of MePFKs under waterlogging stress (Fig. 7B) shows that the expression of MePFPA1 in waterlogged leaves was down-regulated markedly, but that in waterlogged roots was up-regulated considerably. We can speculate that MePFPA1 acts more in waterlogged roots and may be the main force that deals with hypoxic stress. The expression of MePFPB1 in waterlogged leaves was up-regulated at 24 h of waterlogging and down-regulated at 72 h of waterlogging. Therefore, MePFPB1 may act well at the early stage of waterlogging stress. The expression of MePFK03 was lower in leaves than in roots; the expression of MePFK03 in roots under waterlogging stress was up-regulated considerably but was nearly unchanged in waterlogged leaves. MePFK03 was responsive to waterlogged roots.

Phylogenetic analysis of the MePFK family in cassava and four other species
We identified 13 members of the MePFK family in the whole genome of cassava. We submitted all MePFKs with four other PFK families from Arabidopsis, castor bean, rice and potato for phylogenetic tree analysis. The results showed that all PFK genes could be classified into PFK and PFP subfamilies. PFKs could be divided into PFK_A, PFK_B and PFK_C. PFPs could be divided into PFP_α and PFP_β. Amongst the 13 MePFK genes, nine belong to PFKs (five PFK_A, two PFK_B and two PFK_ C) and four belong to PFPs (two PFP_α and two PFP_β). PFP is a heterotetramer with two different α and β subunits. Similar to Arabidopsis, cassava has two genes that encode α subunits, which are the regulatory subunits of PFP, and two genes that encode β subunits, which are the catalytically active subunits of PFP. Similar results were found in Saccharum [22], white pear [23] and rice [20].  The middle segment of the tuber root was divided into three regions: region 1 refers to the periderm, which includes the sclerenchyma, parenchyma and phloem; region 2 refers to the storage parenchyma; and region 3 refers to the tuber centre, which includes xylem vessels and fibres. Values with different letters within one dividual are significantly different at p < 0.05

Expression analysis of MePFKs
The expression analysis of 13 MePFKs showed that MePFPA1 and MePFPB1 were highly expressed in seven tissue parts, including leaf, stem, petiole, fibre root, root cortex, root stele and flower. Therefore, MePFPs play important roles in the whole growth and development of cassava. PFP is involved in plant glycolysis [26][27][28]. PFP provides energy for plant morphogenesis and biochemical reaction through the glycolysis pathway. The expression level of MePFK genes in cassava were lower than that of MePFP genes. Among 9 MePFK genes, MePFK02 and MePFK03 showed relatively high expression. The highest expression for MePFK02 was in the functional leaf, and that for MePFK03 was in root tuber. Therefore, MePFK02 may be involved in glycolysis in leaves and MePFK03 may be involved in glycolysis in the tuber root expansion stage.

Hypoxia occurred in the cassava tuber root
Hypoxia is a common phenomenon in plant tissues. We monitored the oxygen concentration in the roots of two cassava varieties (SC124 and Arg7). The results showed that the oxygen concentration decreased from the outside to the inside until the centre at a drop of below 5%. Similar phenomena were found in potato [29], castor [30], banana [31] and carrot [32]. The mechanism that can make plants conduct a series of biochemical reactions under hypoxia without causing internal anoxia is still unclear.

Waterlogging stress in cassava variety SC124
We set up an experiment under extreme oxygen stress. Roots of cassava SC124 were completely immersed in water for 0, 24, 72 and 168 h. Then, leaves and tuber roots were collected for the transcription level analysis of the key genes of glycolysis. The immediate consequence of waterlogging was that oxygen was blocked and glycolysis flow was promoted under submerged condition. The expression level of MePFK03 in tuber root was up-regulated under waterlogging stress. This finding is similar to that for OsPFK05, which showed a moderate increase in stem and leaves upon anoxia [20]. MePFK03 and OsPFK05 belong to the PFK_A subgroup in the evolutionary tree (Fig. 1). This inducible expression of PFK genes was also found in Arabidopsis (At4g26270 and At4g32480) [21]. Therefore, we can speculate that the induction expression of PFKs is an important regulation for plant metabolism under oxygen deficiency. As for PFP, the expression of MePFPA1 in cassava root was upregulated substantially under waterlogging stress, but the transcripts of MePFPB1 showed up-regulated only at 24 h waterlogging. Similar trends were also observed in rice [20]. OsPFPA was more induced than OsPFPB in anoxic rice seedlings. PFPAs encode the regulatory subunit, and PFPBs encode the catalytic subunit. The adjustment of PFP activity is important under low oxygen stress, and PFPA plays a key regulatory role. The expression levels of different PFK isoforms are distinct in anoxic conditions, and this condition contributes to the balance of glycolysis capacity to cope with low oxygen stress. An oxygen-sensing system exists in plants under normal conditions, but hypoxia will lead to the decrease in adenylate energy level and respiration. Plants can optimise their metabolic pathways under hypoxia by saving ATP and improving the efficiency of oxygen utilisation. Two biochemical pathways from sucrose decomposition and glycolysis are present in plants, namely, the ATPdependent conventional glycolysis and the PPidependent glycolysis. From the point of view of energy consumption, the glycolysis alternative dependent on PPi saves more energy. The three enzymes involved are sucrose synthase (SUSY), PFP and pyruvate orthophosphate dikinase (PPDK), which correspond to the enzymes of conventional glycolysis: invertase (INV), PFK and PK. When hypoxia stress begins, glycolysis will shift from the ATP-dependent pathway to the PPi-dependent alternative pathway [33,34]. In this study, we investigated two genes of PPi-dependent glycolysis bypass. Firstly, sucrose decomposition can be catalysed by two enzymes, namely, the irreversible ATP-dependent INV or the reversible PPi-dependent SUSY. The second key step in glycolysis is the conversion of F-6-P to F-1,6-BP, which can also be catalysed by two enzymes, namely, the ATP-dependent PFK and PPi-dependent PFP. The transcription levels of MeSuSy and MePFP were upregulated substantially in cassava root under waterlogging stress (Fig. 8). The third key step is the conversion of phosphoenolpyruvate to pyruvate, which can be catalysed by the ATP-dependent PK or the PPi-dependent PPDK. The expression of MePPDK in leaves was upregulated after waterlogging treatment, but nearly no expression was observed in waterlogged root tubers [35]. The expression profiles of the above-mentioned enzymes indicate that the PPi-dependent glycolysis pathway was promoted in waterlogged root tubers, and this pathway may save energy and maintain the basic survival needs of plants.
In conclusion, 13 members were present in the MePFK family of cassava. These genes were distributed Error bars indicate the standard deviation calculated from three independent experiments. One and two asterisks (* and **)correspond to significant differences at p < 0.05 and p < 0.01, respectively in 10 of 18 chromosomes. Nine belonged to MePFKs, and four belonged MePFPs. MePFPs showed higher expression in the organs of cassava than MeFPKs. The expression of MePFK03, MePFPA1 and MePFPB1 decreased gradually with the depth increase of cassava tuber. We also found that the expression of MePFPA and MePFK03 increased remarkably in waterlogged roots. MePFPA1 plays an important role in cassava under waterlogging stress. This study will help suggest the candidate genes of how cassava adapts to hypoxic stress. The mechanism of PFP should be further determined by knocking out and overexpressing PFK in cassava root.

Plant materials and treatments
Materials for qRT-PCR. Cassava stalks (SC124 var.) were planted in Chengmai Experimental Field, Chinese Academy of Tropical Agricultural Sciences. The plants were grown for 180 days, and then different tissues from leaves, petiole, stem, root cortex, root stele, fibrous root and flower were collected and frozen in liquid nitrogen for RNA isolation and qRT-PCR. Tuber roots from different developmental stages at 90, 180 and 240 DAP were also collected and frozen in liquid nitrogen for RNA isolation and qRT-PCR.
Materials for waterlogging treatments. Cassava stalks (SC124 var.) were cut and placed in large pots. When the plants grew to 5 months old, the roots of cassava SC124 were completely immersed in water. The leaves and roots were collected at 0, 24, 72 and 168 h after waterlogging initiation and were frozen in liquid nitrogen for subsequent qRT-PCR.

Sequence retrieval and phylogenetic tree construction
The whole protein sequence of the MePFK family was obtained from the cassava genome database ( h t t p s : / / p h y t o z o m e . j g i . d o e . g o v / p z / p o r t a l . html#!info?alias=Org_Mesculenta) [36]. We used the hidden Markov model profile of the PFK protein domains to search the database using the programme HMMER 3.0. The SMART online programme was used to assess the conserved domain of candidate PFKs (http://smart.embl-heidelberg.de/smart/set_mode. cgi? NORMAL = 1). Only proteins with PF00365 domain were regarded as PFK and reserved for further analysis. The multiple alignments of PFK proteins from five different species were made using the Muscle programme of MEGA7 software [37]. A bootstrap neighbour-joining (NJ) phylogenetic tree was established with 1000 bootstrap replicates. The phylogenetic tree of cassava was constructed to better exhibit PFK structure and conserved motifs.

Sequence analyses and protein properties of MePFK and cis-element prediction of MePFK promoters
The isoelectric points and molecular weights of PFK proteins were examined using the ExPASy proteomics server (https://web.expasy.org/protparam/) [38]. PFK gene structures were identified by an online Gene Structure Display Server [39]. The MEME programme (http:// meme-suite.org/tools/meme, Version 5.1.0) [40] was used to investigate the conserved motifs in PFK proteins. In the analysis, the maximum number of motifs was 15, and the optimum width of motifs was 6-50.
The promoter sequences (less than 2000 bp) of MePFKs were obtained from the cassava genomic sequence downloaded from the JGI database. These sequences were submitted to the new PLACE database (https://www.dna. affrc.go.jp/PLACE/?action=newplace) for the analysis of the cis-elements of MePFK promoters [41].

Chromosomal and subcellular localization of MePFK protein in cassava
The chromosomal locations of MePFK were determined on the basis of the chromosomal information derived from the JGI database. The position of MePFK was physically mapped to a chromosome using the Mapchart software according to the relative location of MePFK on the chromosome [42]. The subcellular localisation of the MePFK protein was predicted online using the Cell-PLoc 2.0 package [43]. We selected two genes for subcellular localisation by GFP fusion protein expression to confirm the predicted results. The CDS fragments of MePFK03 and MePFPA1 without stop codon were amplified using cDNA from Ku50 leaves by gene-specific primers (Table S1) designed according to the CDS from the JGI database. The fragments were cloned into the binary vector pCAMBIA1300R (modified based on pCAMBIA1300). Positive clones were confirmed by DNA sequencing, and the expression vectors of 35S:: MePFK-GFP and 35S::MePFPA1-GFP were generated.
The expression vectors were transformed into A. tumefaciens EHA105 and then infiltrated into the leaves of Nicotiana benthamiana for transient expression. After 2 days, the infiltrated leaves were placed on a confocal laser scanning microscope to observe the GFP fluorescent signals.

Oxygen concentration determination in cassava tubers
Oxygen concentration was analysed from two cassava varieties, namely, Arg7 and SC124. Three individual plants were collected from each variety, two tubers were sampled per plant, three regions were selected per tuber, and each region was tested 10 times. Internal oxygen concentration was measured 1-2 min later by inserting an O 2 microelectrode (< 1 mm tip diameter; Presens, Germany) into the tuber root tissue. The measured air oxygen concentration was measured 22%. Tissues were sampled at the same region for the transcript analysis of MePFKs.

Transcript analysis of MePFKs in cassava plants
RNA extraction and cDNA synthesis. Total RNA was extracted from cassava tissues through the sodium dodecyl sulfate method [44], and purified by NucleoTrap® mRNA kit (Macherey Nagel, Genmany) according to the manufacturer's instruction. First strand cDNAs were synthesized by SMARTScribe™ Reverse Transcriptase from 200 ng of total RNA using primers of CDS III and Oligo III according to product manual.
The expression levels of MePFKs from different cassava tissues and from roots at different developmental stages were investigated by qRT-PCR analysis. qRT-PCR was performed in the thermal cycler of a Rotor-Gene 6000 (Rcorbett, Australia) using a SYBR Premix Ex TaqTM fluorescence quantitative kit (TaKaRa, Japan). Tubulin-F: 5′GTGGAGGAACTGGTTCTGGA3′ and Tubulin-R: 5′TGCACTCATCTGCATTCTCC3′ were used as reference gene [45]. The gene-specific primers are shown in Table S1. The PCR programme proceeded as follows: 95°C for 30 s, 40 cycles of 95°C for 10 s, 59°C for 20 s and 72°C for 30 s. Each sample had three technical replicates. The relative expression level was determined by the 2 −Δ Δ Ct method [46].

Statistical analysis
Three biological replicates were used for each measurement. The data are represented as mean ±standard deviation. Significant differences between different samples were tested with the software IBM SPSS statistics 19.0 [47].