s2fieldsofstudy: Biology Medicine

text:
MCM6 promotes metastasis of hepatocellular carcinoma via MEK/ERK pathway and serves as a novel serum biomarker for early recurrence

Background The high incidence of recurrence and metastasis of hepatocellular carcinoma (HCC) necessitate the discovery of new predictive biomarkers of invasion and prognosis. Minichromosome maintenance complex component 6 (MCM6), which has been reported to up-regulate in multiple malignancies, was considered to be a novel diagnoses biomarker in HCC. However, its functional contributions and prognostic value remain unclear. Methods The expression of MCM6 was analyzed in 70 HCC tissues and 5 HCC cell lines by immunohistochemistry and real-time RT-PCR. The roles of MCM6 in HCC cell proliferation, migration and invasion were explored by CCK8, Wound healing and Transwell assays, respectively. Western blotting and Immunofluorescence staining were conducted to detect the protein expressions of ERK signaling pathway and EMT-related markers. To verify the above findings in vivo, we established subcutaneous xenograft tumor and orthotopic xenograft tumor models in nude mice. Finally, Enzyme-linked immunosorbent assay was used to evaluate the serum MCM6 level. Results MCM6 was significantly up-regulated in HCC tissues. Increased MCM6 expression was associated with aggressive clinicopathological features and worse prognosis in HCC patients. These results were consistent with our analyses of The Cancer Genome Atlas database (TCGA). Furthermore, knockdown of MCM6 significantly decreased proliferative and migratory/invasive capability of HCC cells in vitro, as well as decreased tumor volume, weight and the number of pulmonary metastases in vivo. Mechanistic analyses indicated that MCM6 promoted EMT and activated MEK/ERK signaling. More importantly, serum MCM6 levels in HCC patients were significantly higher than those in cirrhosis and healthy controls (P < 0.0001), and allowed distinguishing early recurrence with high accuracy (AUC = 0.773). Conclusions Our findings indicate that MCM6 predicts poor prognosis and promotes metastasis in HCC. Postoperative serum MCM6 level could be valuable to detect preclinical early recurrence, indicative of a need for more careful surveillance and aggressive therapeutic intervention. Electronic supplementary material The online version of this article (10.1186/s13046-017-0669-z) contains supplementary material, which is available to authorized users.

Background
Hepatocellular carcinoma (HCC) is one of the most prevalent malignant tumors worldwide. Although surgical techniques and adjuvant therapy have improved, long-term survival remains low due to tumor recurrence and metastasis [1]. Early recurrence is defined as HCC that recurs within 1 year after resection, and it is most closely related to cancer metastasis spread and is the leading cause of early death in patients with HCC after liver resection [2]. Therefore, it is critical to understand the molecular mechanism associated with HCC invasion and metastasis and identify the predictive biomarkers of HCC recurrence and prognosis to help monitor disease progression and guide diagnosis and treatment.
The Mini-chromosome maintenance protein (MCM) family is a recently reported gene family that was identified in the control of eukaryotic genome replication. In the DNA replication licensing process, the MCM2-7 complex primes chromatin for DNA replication by binding origins of DNA replication during the late M to early G1 phase of the cell cycle [3]. Also, it is involved in the formation of replication forks and in the recruitment of other DNA-replication-related proteins. The MCM complex is a replicative helicase that is essential for DNA replication initiation and elongation in eukaryotic cells [4].
The MCM protein family plays an important role in genome duplication of proliferating cells. It had been reported that abnormal expression of MCMs contributes to tumorigenesis. For example, a significant rise in the MCM5 level was observed in cervical epithelium squamous cancer [5]. MCM7 expression may predict poor postoperative prognosis for HCC and promote cancer progression through cyclin-D1-dependent signaling [6,7]. Altered Mcm5, Mcm2 and Mcm7 were found in malignant cells in the lung, kidney and prostate [8][9][10][11]. Furthermore, a rise of MCM6 has been shown to correspond with a high risk of recurrence in meningioma [12], suggesting that the MCM proteins are promising prognostic markers for monitoring various cancers.
MCM6 is one of six members of the Mini-chromosome maintenance family. Previous studies verified the important role of MCM6 in cell proliferation and indicated its potential role in prompting tumor progression. MCM6 is up-regulated in various types of tumors [13][14][15]. It has been reported that a high serum MCM6 level is a diagnostic biomarker for HCC [16]. However, little is known about its relation to cancer survival or its prognostic significance in HCC. The previous studies give rise to the hypothesis that MCM6 plays an important role in HCC. In the present study, we performed a series of experiments to explore its molecular function in HCC cells and evaluate the prognostic value of both the tissue and plasma MCM6 expression level.

Tissue and serum specimens
Tumor tissue specimens, including tumor and matched adjacent non-tumor tissues, were obtained from 70 patients who had undergone curative liver resection. Detailed clinical pathological parameters were listed in Table 1. Thirty  Table S1.

Immunohistochemistry assay
Immunohistochemical staining of tissue was performed according to the manufacturer's instructions (DAKO, Glostrup, Denmark) on formalin-fixed, paraffin-embedded tissue sections that had been cut to 4-μm thickness. The specimens were dried at 62°C for 2 h and then dewaxed in xylene and rehydrated using graded alcohols, followed by incubation in 3% hydrogen peroxide for 15 min to exhaust the endogenous peroxidase activity. The antigens were then retrieved in Tris-EDTA (pH 9.0)/ 0.01 M sodium citrate buffer (pH 6.0) using a microwave oven for 25 min, followed by blocking with 10% goat serum for 30 min to prevent non-specific staining. Primary antibodies were incubated overnight at 4°C in a humidified chamber, followed by HRP-labeled anti-mouse/rabbit secondary antibody (DAKO) incubation for 1 h at room temperature. Antibody binding was detected by DAB, and the reaction was stopped by immersion in distilled water once the brown color appeared. Finally, specimens were counter-stained with hematoxylin for 3 min. Two different pathologists who specialized in liver cancer evaluated the IHC results without knowledge of the clinical data. Both the extent and intensity of immunostaining were taken into consideration when analyzing the data. The intensity of staining was scored from 0 to 3, and the extent of staining was scored from 0% to 100%. Final quantitation of staining was obtained by multiplying the two scores. The protein expression was classified as high expression if the score was higher than 1.5 and low expression if the score was 1.

Immunofluorescence staining
Cells grown on coverslips were rinsed with phosphatebuffered saline (PBS) and fixed with cold 4% paraformaldehyde for 5 min at RT. Subsequently, the cells were blocked with Triton X-100 at a concentration of 0.2% for 30 min. Cells were then blocked for 1 h with 5% BSA and washed for 30 min, followed by incubation with primary monoclonal antibodies against MCM6 (1:500; Abcam, Cambridge, UK). Vimentin (1:100, ZSGB-BIO, China) and E-cadherin (1:50, ZSGB-BIO, China) overnight at 4°C. The next day, the coverslips were incubated for 1 h in a dark room with Alexa Fluor 564 goat anti-rabbit IgG and Alexa Fluor 488 goat anti-mouse IgG (1:100 dilution; Bioworld Technology, Inc). Furthermore, the coverslips were stained with DAPI for 5 min at 4°C. Finally, an LSM80 fluorescent microscope (ZESS, Germany) was used to observe the expression in cells.

Western blot
Total proteins were extracted with RIPA lysis buffer supplemented with protease inhibitors, separated by SDS-PAGE and transferred to a PVDF membrane (Millipore, Billerica, MA). The membrane was blocked with 5% non-fat milk at room temperature for 1 h and incubated with the appropriate antibody. Protein bands were visualized by the chemiluminescent HRP detection system

Cell viability assay
Cell viability was analyzed using the Cell Counting Kit-8 (CCK-8) (Dojindo Molecular Technologies Inc. Shanghai, China). Cells at a density of 2 × 10 3 /well were seeded into 96-well plates and cultured in 100 μL of DMEM containing 10% FBS for 4 days. Ten microliters of CCK8 solution was added to each plate, and the cells were incubated for 3 h at 37°C. The absorbance value (OD) of each well was measured at 450 nm.

EDU incorporation assay
For the EdU incorporation assay, proliferating cells were examined using the Cell-Light EdU Apollo 567 in vitro Imaging Kit (RiboBio, Guangzhou, China) according to the manufacturer's protocol. Cells at a density of 5 × 10 3 /well were seeded into 96-well plates and incubated with 10 mM EdU for 2 h, followed by fixation with 4% paraformaldehyde, permeabilization in 0.2% Triton X-100 and staining with Apollo fluorescent dyes. Next, 50 μL/well DAPI was used to stain the cell nuclei for 10 min. The number of EdU-positive cells was counted under a fluorescent microscope in five random fields.

Wound-healing assay
For the wound-healing assay, cells were grown to confluence in a six-well plate. Artificial wound tracks were created by scraping the confluent cell monolayers with a pipette tip. The cells were fed with serum-free medium.
The ability of the cells to migrate into the wound area was assessed every 24 h after scratching.

Colony formation assay
Cells were plated in 6-well culture plates at 100 cells/ well. After incubation for 2 weeks at 37°C, the cells were washed twice with PBS and stained with 0.1% crystal violet solution. The number of colonies containing ≥50 cells was counted under a microscope. The colony formation efficiency was calculated as (number of colonies/number of cells inoculated) × 100%.

Cell migration and invasive assays
Cell motility was assessed by cell migration and invasion assays using transwell chambers with or without Matrigel (BD, Biosciences, CA). Approximately 6 × 10 4 cells in medium without FBS were seeded on transwell chambers with or without Matrigel and incubated at 37°C for 15 h. Medium containing 10% FBS was put in the lower chamber. The invasive cells attached to the lower surface of the membrane insert were fixed, stained using Giemsa (Jiancheng, Jiangsu, China) and quantified.

Enzyme-linked immunosorbent assay (ELISA)
The serologic level of MCM6 was detected using the human MCM6 ELISA kit (Cusabio, Wuhan, CN) according to the manufacturer's instructions. One hundred microliters of serum samples and standard proteins was added to the appropriate wells and incubated at 37°C for 2 h. After being washed three times with washing buffer, HRP-conjugated secondary antibody (100 μL/well) was added to the cells, which were incubated at 37°C for 1 h. The plate was washed three times, and the cells were then incubated with DAB solution for 15-30 min at 37°C. Then, 50 μL of Stop Solution was added to each well, the optical density at 450 nm was detected, and the concentration of MCM6 was calculated from the standard curve.

Animal studies
Both subcutaneous and orthotopic models were used for animal studies. A total of 2 × 10 6 SMMC7721 cells transfected with MCM6 shRNA, lentiviral vectors and negative control (scramble) vector in 0.1 ml PBS medium was used in both the subcutaneous model and orthotopic model. The cells were injected into the left leg of 4-6-week-old male BALB/c nude mice. The subcutaneous tumor size was calculated and recorded every 3 days with a Vernier caliper using the following equation: tumor volume (mm 3 ) = (L × W 2 )/2, where L = tumor long axis and W = short axis. The measurements were repeated three times. The mice were maintained in a barrier facility on HEPA-filtered racks and fed an autoclaved rodent diet. After 27 days, the mice were killed and the tumor tissues were surgically excised, weighed and stained with hematoxylin and eosin (H&E). For the orthotopic model, the cell suspension was injected into the left hepatic lobe of 4-week-old nude mice with a micro-syringe. Transplanted cells were allowed to grow for up to 6 weeks, when the mice were sacrificed. The livers and lungs were dissected, fixed, and paraffinembedded. All animal handling and procedures were approved by the Animal Care and Use Committee of Southern Medical University.

Statistical analysis
All statistical analyses were performed with SPSS statistical software (version 21.0; SPSS, IBM, Armonk, NY). Survival curves were constructed using the Kaplan-Meier method and analyzed by the log-rank test. Significant prognostic factors found by univariate analysis were entered into a multivariate analysis using the Cox proportional hazard regression model. The Pearson's chi-square test was used to analyze the association of MCM6 expression with various clinicopathologic characteristics. A paired Student's t test or χ2 tests was used to compare the values between subgroups. Receiver operating characteristic (ROC) curve analysis was used to determine the predictive value of the parameters, and the differences in the area under the curve (AUC) were detected. R project for Statistical Computing (R version 3.4.3) was used to calculate Harrell C-index. Data were expressed as the mean ± SD of at least three sample replicates. A value of P < 0.05 was statistically significant, and single, double, and triple asterisks indicate statistical significance of * P < 0.05, ** P < 0.01 and *** P < 0.001.

MCM6 expression is significantly up-regulated in HCC tissues and correlates with poor prognosis and high recurrence rates in HCC patients
To evaluate the expression of MCM6 in HCC tissues, we first conducted immunohistochemistry assay in 70 HCC specimens (tumor and matched adjacent nontumorous tissues) and 30 normal hepatic tissue specimens. We observed that MCM6 was primarily localized in the nucleus. The results revealed that MCM6 expression in HCC tissues was markedly higher than the adjacent non-tumor liver tissue and normal liver tissue (Fig. 1a). We detected high expression of MCM6 in 50/ 70 (71.4%) HCC tissues, compared with only 8/70 (11.4%) adjacent non-tumor tissues and no staining in normal tissues (P < 0.01; Fig. 1b).
Second, to verify the functions of MCM6, we tested the correlation of MCM6 expression status in 70 HCC specimens with 12 widely recognized clinicopathological parameters. Statistical analysis indicated that high MCM6 expression was associated with hepatitis B infection, liver cirrhosis, multiple tumor nodules, high serum AFP and early recurrence (Table 1). Furthermore, univariate analysis revealed that tumor size, GGT level, vascular invasion, early recurrence and MCM6 expression were unfavorable predictors for TTR and OS (Table 2). To assess the prognostic significance of MCM6, we do the Kaplan-Meier survival analysis. The 1-, 3-and 5year recurrence rates in the MCM6-positive group were significantly higher than in the MCM6-negative group (58.0%, 80.1% and 81.5% versus 19.2%, 55.1% and 61.9%, respectively; Fig. 1c). Similarly, the 1-, 3-and 5-year OS rates in the MCM6-positive group were significantly lower than in the MCM6-negative group (61.2%, 38.3% and 32.4% versus 95.1%, 70.7% and 45.1%, respectively; Fig. 1d). Further multivariate analysis demonstrated that MCM6 was an independent prognostic factor for both the OS (P < 0.05) and TTR (P < 0.05) of HCC patients (Table 3).
In addition, to confirmed the prognostic value of MCM6, we applied MCM6 to 291 HCC samples with RNA-sequencing data in the Cancer Genome Atlas project (TCGA). Based on the relative expression of MCM6, 124 and 167 patients were classified into high-and low-expression groups, respectively. We performed systematic Kaplan-Meier survival and clinicopathological features analysis of the whole mRNA transcriptomics and found that high MCM6 mRNA expression was significantly correlated with early recurrence (recurrence rates 80.2% vs. 53.2%) (Fig. 1e) and predicted poor survival and high recurrence risk ( Fig. 1f-g).

Knockdown of MCM6 suppresses HCC cell migration, invasion and proliferation in vitro
In vitro experiments suggested that MCM6 expression in high-metastasis potential cell lines, such as HCC-LM3 and SMMC7721, was higher than in the lowmetastasis potential cell lines PLC/PRF5 and HepG2 (Fig. 2a). These data indicated that MCM6 might regulate the invasion and metastasis of HCC cells. This speculation was confirmed by cytological experiments. To understand the function of MCM6 in HCC cells, we designed three short hairpin RNAs (shRNA1, shRNA2, and shRNA3) to silence MCM6 expression. More than fourfold decrease in MCM6 expression was observed in SMMC7721 and HCC-LM3 cells treated with shRNA compared with the control group by qRT-PCR (Student's t-test, with P < 0.05 for both) and western blotting; shRNA3 was the most effective shRNA and was chosen for further study (Fig. 2b-c). The wound-healing and Matrigel-coated (for invasion) or -uncoated (for migration) Transwell assays showed that HCC invasion and metastasis were effectively inhibited by MCM6 knockdown. The results showed that SMMC7721-scramble and HCC-LM3-scramble cells had a faster wound closure rate and more invasive properties than shMCM6 cells (Fig. 2d-e). These results uncovered the role of MCM6 in promoting HCC progression in vitro.

MEK/ERK signaling activation is critical for MCM6-induced EMT
Invasion and metastasis are responsible for the unsatisfactory long-term survival. Emerging evidence suggested that EMT endows epithelial cells with their migratory and invasive capacity and is involved in many cancer metastasis processes, including HCC [17,18]. To explore the relationship between MCM6 and the epithelialmesenchymal transition (EMT), we performed real-time PCR to assess EMT markers in MCM6 knockdown cells. We found that the suppression of MCM6 resulted in higher E-cadherin and ZO-1 levels, accompanied by decreased levels of the mesenchymal marker vimentin, but did not significantly influence the other EMT markers (Fig. 3a). Western blot and immunofluorescence staining also confirmed that MCM6 expression inversely correlated with E-cadherin but positively correlated with vimentin expression (Fig. 3b and c).
Next, we sought to determine the signaling mechanisms involved in MCM6-mediated EMT. The extracellular signal-regulated kinase 1/2 (ERK, also known as p42/44 mitogen-activated protein kinase MAPK) signaling pathway is one of the most intensively studied signaling mechanisms that regulates migration, invasion, and EMT in cancer [19][20][21]. We found that MCM6 affected activation of ERK signaling. We analyzed the activation status of ERK and the upstream gene MEK1/2 in cells with stable MCM6 knockdown. The data showed that MCM6 increased the levels of phosphorylated MEK1/2 (p-MEK1/2) and ERK1/2 (p-ERK1/2), but the total MEK and ERK had similar expression levels in all experimental conditions (Fig. 3d). Meanwhile, we used immunofluorescence assays to detect the localization and expression of ERK protein.
It is widely known that ERK shuttling between the nucleus and cytoplasm provides a switch-like transition effect in the activation of the ERK signal, and active nuclear ERK1/ 2 levels were assessed by immunofluorescence. Consistent with the western blot data, the confocal photos revealed that the activation of ERK1/2 was strongly suppressed when MCM6 was knocked down (Fig. 3e). Taken together, these data suggest that MCM6 promotes EMT, at least in part, through the upregulation of MEK/ERK signaling pathways.

MCM6 augments tumor-initiating and metastatic potential in vivo
To further investigate the effects of MCM6 in vivo, we established subcutaneous xenograft tumor and orthotopic xenograft tumor models in nude mice. After 4 weeks, SMMC7721-Scramble cell-derived tumors at the subcutaneous implantation sites were larger and grew more rapidly than SMMC7721-shMCM6 cell-derived tumors. (Fig. 4a). Consistently, the liver orthotopic xenograft tumor also showed that MCM6 knockdown had imposed significant regulatory effects on HCC metastasis. The intrahepatic and pulmonary metastasis rates in mice   with tumors derived from SMMC7721-shMCM6 cells were significantly lower than the corresponding control group (Fig. 4b and c). Moreover, the effect of MCM6 on Ecadherin and vimentin was also observed by IHC in the xenograft specimens (Fig. 4d). We found that MCM6 expression inversely correlated with E-cadherin but positively correlated with vimentin and Ki-67 expression. Taking these results together, our studies demonstrated that MCM6 could promote HCC growth and metastasis in vivo.

Elevated serum MCM6 lever is correlated with increased risk of early recurrence after curative hepatectomy
Previous studies demonstrated that the MCM6 levels in plasma had diagnostic value for HCC [16]. Additionally, our earlier results revealed that the elevated expression of MCM6 in tissues and serum is correlated with early recurrence (Table 1) (Additional file 1: Table S2). Therefore, to further explore its potential value in early recurrence non-invasive serologic detection, we examined MCM6 expression in patient serum. Firstly, when discerning HCC patients from cirrhosis and healthy controls, serum MCM6 had significant sensitivity (P < 0.0001, Fig. 5a). Moreover, the data obviously revealed that the serum MCM6 concentration from HCC patients who received curative hepatectomy dropped sharply compared with the pre-treatment level (2039.17 ± 206.02 pg/mL preoperative versus 1077.51 ± 234.31 pg/ mL postoperative, P < 0.05). In addition, in the nonearly recurrence subgroup, plasma MCM6 levels declined more dramatically after operation and some were even reduced to the normal range over a certain time (Fig. 5b). To analyze the value of serum MCM6 in predicting the prognosis of HCC, we prepared the receiver operating characteristic (ROC) curve for serum MCM6 expression according to early recurrence. The cut off value was determined using receiver operating characteristic (ROC) curves and the Youden's Index. (Youden's index = sensitivity + specificity-1). According to Youden's index, optimal sensitivity (78.6%) and specificity (76.5%) were seen at a cut-off of ≥941.635 pg/ml to distinguish the non-early recurrence group from the early recurrence group. To estimate the value of the dynamic changes in serum MCM6 in the early recurrence of HCC after liver resection, we analyzed the differences in the recurrence rate of the two groups. Early-recurrence was observed in 8 of 10 (80%) patients with high postoperative MCM6 levels (≥ 941.635 pg/ml), whereas only 6 of 21 (28%) patients with low postoperative MCM6 levels (< 941.635 pg/ml) had early-recurrence (Fig. 5c). Additionally, a high postoperative MCM6 level correlated significantly with cumulative recurrence after curative hepatectomy (Fig. 5d).
Next, the predictive performance was compared between serum preoperative MCM6 and postoperative MCM6 (Fig. 5e). Our data showed that after curative hepatectomy, the area under the ROC curve (AUROC) for MCM6 was 0.773 (95% confidence interval [CI], 0.594-0.953). The postoperative serum MCM6 was Fig. 4 MCM6 augments tumor proliferation and metastatic potential in vivo. a Xenografts in nude mice were established by subcutaneous injection of SMMC7721-shMCM6 and control cells. b and c Representative pictures and external whole-body fluorescence images of intrahepatic and lung metastasis, which were obtained 6 weeks after orthotopic liver injection, respectively. *P < 0.05 based on Student's t test. Scale bar: 5 mm (top), 100 μm (bottom). d Representative IHC images of MCM6, as well as Ki67 and EMT-associated genes in subcutaneous tumors of mice injected with SMMC7721-shMCM6 and control cells. Scale bars: 50 μm more informative than the preoperative level in predicting early recurrence in HCC. Furthermore, the ROC analysis also revealed that postoperative MCM6 is superior to postoperative AFP in differentiating early recurrence. (Fig. 5f ). The AUC of postoperative MCM6 was 0.773 (95% confidence interval [CI], 0.594-0.953), which was larger than that of AFP at 0.681 (95% confidence interval [CI], 0.474-0.887). Additionally, the sensitivities of the postoperative MCM6 and AFP were 78.6% and 71.4%, respectively. Notably, the sensitivity of the combination with AFP and MCM6 was increased to 80.1% (Table 4). Harrell C-index could also proved that postoperative serum MCM6 not only significantly correlated with tumor early recurrence but also displayed improved indices of prognostic performance. Harrell C-index = 0.7249 95%CI (0.6143-0.8356). This result showed that the prediction model is accurate and the credibility is high (a C-index of 1 indicates 100% predictive accuracy). To summarize, in patients with resectable HCC, high postoperative serum MCM6 individuals had a higher tendency of early recurrence. In particular, patients with dual positivity of postoperative serum MCM6 and AFP Fig. 5 Elevated perioperative plasma MCM6 correlates significantly with increased risk of early recurrence after curative hepatectomy. a Compared with healthy control individuals and cirrhotic patients, serum MCM6 was significantly higher in patients with HCC (P < 0.0001). b In the non-early recurrence subgroup, plasma MCM6 levels decreased dramatically after operation and some were even reduced to the normal range in a certain time. c The high serum MCM6 subgroup was more likely to suffer from early recurrence than the low MCM6 subgroup after surgical resection (80% versus 28%). d A high serum MCM6 level correlated significantly with cumulative recurrence after curative hepatectomy (P < 0.0001). e-f ROC curve by serum biomarkers indicated that postoperative serum MCM6 has a significant predictive ability of HCC early recurrence: AUC postoperative MCM6 > AUC preoperative MCM6 (0.773 versus 0.513), and AUC postoperative MCM6 > AUC postoperative AFP (0.773 versus 0.681). ROC, receiver operating characteristic. AUC, area under the curve should be supplied with more intensive care and close follow-up after they undergo tumor resection.

Discussion
Invasion and metastasis are responsible for the unsatisfactory long-term survival of HCC patients. An early diagnosis and aggressive intervention may offer HCC patients a significant survival benefit. Thus, the identification of new predictive biomarkers of HCC invasion and prognosis is critical. In this study, we provided evidence that MCM6 had the potential to be a novel prognostic biomarker for HCC patients. IHC assays suggested that MCM6 was highly expressed in HCC tissues. Clinical significance analysis indicated that MCM6 was associated with tumor number, liver cirrhosis and early recurrence. Also, HCC patients with increased MCM6 expression had worse over survival and higher cumulative recurrence rates. All of these results suggested that MCM6 might play an important role in HCC development. In addition, depletion of endogenous MCM6 expression in LM3 and SMMC7721 cells suppressed cell growth, migration, and invasion. Mechanical analysis indicated that MCM6 promoted EMT and activated the MEK/ERK pathway. Moreover, we found that high serum MCM6 levels significantly correlated with high early recurrence risk after hepatectomy. All of these results strongly suggest that MCM6 has the potential to be a useful prognostic biomarker for HCC.
Accumulated data show that EMT results in increasing cell migration and invasion in several cancers [22,23]. Functional experiments revealed that HCC invasion and metastasis were effectively inhibited by MCM6 knockdown. More significantly, we found that MCM6 knockdown increased epithelial marker expression and decreased mesenchymal marker expression. Additionally, the relationship of MCM6 and EMT markers was confirmed in xenografted tissue specimens. Therefore, our results illustrated that MCM6 promoted the EMT process. Previous reports showed that activation of the MEK/ERK signaling pathway contributed to cell growth, invasion, and EMT [24,25]. We showed that phosphorylation of ERK was inhibited by MCM6 silencing. This result suggests that MCM6 might affect EMT via regulation of the activation of the ERK signaling pathway. Strategies to target this pathway might be developed to inhibit HCC metastasis, although more mechanistic studies are needed.
Our study found that MCM6 was significantly upregulated in HCC tissues and predicted poor prognosis and high recurrence risk. Data from the TCGA confirmed this result. Since liquid biopsy has been widely applied in diagnostic situations and the detection of many cancers [26,27], serum biomarkers are also recommended to be used in surveillance for HCC early recurrence. Interestingly, our data indicated that the serum MCM6 protein level can predict early recurrence of HCC patients who accepted radical resection. Plasma analyses suggested that postoperative plasma MCM6 levels of early-recurrence patients were significantly higher than those of non-early-recurrence patients. HCC tissues are not always readily obtained during routine follow-up visits, so MCM6 has meaningful potential to be translated and applied to serum detection although further verification based on large sample tests is needed.

Conclusions
MCM6 functions as a tumor promotor by activating MEK/ERK signaling and subsequently regulates metastasis and EMT in HCC. In addition, MCM6 can be a potential indicator to predict poor prognosis and early recurrence after receiving curative hepatectomy. These efforts to provide personalized surveillance are expected to improve the overall clinical management of HCC.

Additional file
Additional file 1: Table S1. The primers used in this study. Table S2.


extfieldsofstudy: Medicine Chemistry Biology

text:
Thermodynamic constraints on the regulation of metabolic fluxes

Nutrition and metabolism are fundamental to cellular function. Metabolic activity (i.e. rates of flow, most commonly referred to as flux) is constrained by thermodynamics and regulated by the activity of enzymes. The general principles that relate biological and physical variables to metabolic control are incompletely understood. Using metabolic control analysis and computer simulations in several models of simplified metabolic pathways, we derive analytical expressions that define relationships between thermodynamics, enzyme activity, and flux control. The relationships are further analyzed in a mathematical model of glycolysis as an example of a complex biochemical pathway. We show that metabolic pathways that are very far from equilibrium are controlled by the activity of upstream enzymes. However, in general, regulation of metabolic fluxes by an enzyme has a more adaptable pattern, which relies more on distribution of free energy among reaction steps in the pathway than on the thermodynamic properties of the given enzyme. These findings show how the control of metabolic pathways is shaped by thermodynamic constraints of the given pathway.

Metabolism enables the utilization of nutritional resources to provide free energy, material, and cellular communication for functions of cells (1). Some metabolic pathways such as central carbon metabolism, which processes macronutrients or the major caloric sources in diet (proteins, fats, and carbohydrates), have been largely defined for over 50 years (2), and even genome-scale reconstructions of metabolic networks that include thousands of metabolites and reactions in different intracellular compartments are available in many unicellular organisms and metazoans (3)(4)(5). Despite the complexity of metabolic networks in living organisms, metabolic processes still follow the laws of thermodynamics. It is thus important to understand the principles that link fundamental thermodynamic quantities to the control of metabolic pathway activity (i.e. the rates of flow of materials through the network, most commonly referred to as fluxes).
Analytical frameworks have been developed to understand steady-state and dynamic behaviors of metabolic networks. The most widely applied method in computational modeling of metabolism is flux balance analysis, which assumes that the network is in steady state and that configurations of metabolic fluxes are determined by optimizing an objective function such as growth rate (6). This approach does not require that enzyme properties such as expression levels be known, but the utility of the method depends on the objective function and other assumptions such as nutrient uptake rates. Moreover, flux balance analysis does not provide information about pathway regulation and the control of metabolic flux, which requires additional knowledge such as enzyme activities or metabolite abundances, which are routinely measurable (7,8).
Another framework for understanding metabolism is metabolic control analysis (MCA), 3 originally developed in the 1970s once the biochemistry for many of the key pathways in metabolism such as glycolysis and the tricarboxylic acid cycle was established (9 -15). MCA quantitatively measures how the flow through a metabolic pathway responds to changes in parameters such as the abundance of an enzyme or availability of a nutrient. MCA defines the sensitivity of a metabolic flux to a perturbation in a given metabolic reaction (i.e. the flux control coefficients (FCCs)) and also provides a series of rigorously derived relationships between metabolic fluxes, metabolite concentrations, and enzyme activities. This framework can also be applied in computer simulations (16 -19) or in experimentation (20,21) when some but not all of the relevant variables are measured (22)(23)(24)(25).
Metabolism is subject to the laws of thermodynamics. These laws place constraints on the dynamics of metabolic reactions (26 -32) and are known to affect the control of fluxes in linear pathways (9,31,33,34). Moreover, the deviation from equilibrium at an individual reaction step has been applied as the criterion to identify rate-limiting steps in metabolic pathways (2,29). However, this rule of thumb has been challenged by MCA, at least for linear pathways (9), but a quantitative evaluation of the relationship between thermodynamics and flux control in metabolic pathways with different topologies such as the control at branching points is still lacking.
In this study, we use a set of models with representative topological structures observed in metabolism to investigate the quantitative relationships between thermodynamics and regulation of metabolic fluxes. Notably, all models are based on This work was supported by National Institutes of Health Grants R01 CA193256 and R00CA168997 and American Cancer Society Grant TBE-TBE434120 (to J. W. L.). The authors declare that they have no conflicts of interest with the contents of this article. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This article contains supporting information and Figs. S1 and S2. 1 To whom correspondence may be addressed. E-mail: ziwei.dai@duke.edu. 2 To whom correspondence may be addressed. E-mail: jason.locasale@ duke.edu. some simplifying assumptions, but some (and hopefully these) are useful. To our surprise, we find that, in both linear and branched pathways, the regulation of pathway fluxes by individual enzymes is in general loosely constrained by the deviation from thermodynamic equilibrium, or the thermodynamic driving force, of the whole pathway. Only pathways very far from equilibrium have their flux regulation strictly constrained by the thermodynamic driving force, in which all fluxes are almost fully controlled by the upstream enzymes. The relationships are further studied in a mathematical model of glycolysis as an example of a more complicated metabolic model. These results unravel simple principles of how metabolic pathways are regulated by the interaction between thermodynamics and enzyme activity.

Metabolic flux and thermodynamics
We first consider unimolecular, first-order kinetics of reversible metabolic reactions for simplicity. In this framework, each reaction has one substrate (S) and one product (P) (Fig. 1). It is noteworthy that this simplified case approximates the more complicated Michaelis-Menten mechanism when the abundance of substrate is far below the K m . The forward reaction rate v ϩ ϭ k ϩ S and backward reaction rate v Ϫ ϭ k Ϫ P are linear in substrate and product concentrations. Because the rate constants of the forward and backward reactions are coupled by the equilibrium constant K:K ϭ k ϩ /k Ϫ , the net flux carried by this reaction is as follows: v ϭ v ϩ Ϫ v Ϫ ϭ k ϩ (S Ϫ P/K). Let k ϭ k ϩ for simplicity, and we have the following.
The equilibrium constant K can be further connected to the standard Gibbs free energy of this reaction.
Finally, the Gibbs free energy is determined by the standard Gibbs free energy in combination with concentrations of the substrate and product.
We note that Equation 1.1 can also be rewritten to explicitly incorporate the reaction free energy change. Here, we let g ϭ ⌬G/RT for simplicity. Thus, we have the following.
According to the connectivity theorem in MCA (13), flux control coefficients of reactions directly associated with a metabolite are coupled through local elasticity coefficients of these reactions with respect to the metabolite. For a reaction and a metabolite directly associated with this reaction (either as substrate or as product), the elasticity coefficient is defined as the partial derivative of the reaction rate with respect to the concentration of the metabolite on the logarithmic scale. The elasticity coefficients can be computed from Equations 1.1 and 1.4.

Linear pathways
Based on the definitions under "Metabolic flux and thermodynamics," we now consider a linear pathway consisting of n reactions ( Fig. 2A). Each reaction has unimolecular, linear kinetics, as described previously. The first substrate S in is converted to the end product S out by this reaction chain in n steps with n Ϫ 1 intermediary metabolites S 1 , . . . , S n Ϫ 1 . The ith reaction in this chain has rate constant k i and equilibrium constant K i . The concentrations of S in and S out are treated as fixed parameters, because these two values, together with rate constants and equilibrium constants of the enzymes, determine concentrations of all other metabolites.
Previous work has shown that in linear pathways with either first-order or zero-order (i.e. substrate-saturated) kinetics, the pattern of flux control can be completely determined if free energies of all reactions involved in the pathway are known (34). Here, we use these results as necessary background to study the relationship between thermodynamics and flux control in metabolic networks with different topological structures. Briefly, from MCA (13), there are analytical relationships between thermodynamic properties and the flux control coefficients, which quantify the relative importance of each enzyme in regulating the flux through the pathway.
To apply the summation theorem and connectivity theorems in the theory of metabolic control analysis, we assume that the A, forward and backward fluxes of a reversible reaction. S is the concentration of substrate, P is the concentration of product, v ϩ is the rate of the forward reaction, v Ϫ is the rate of the backward reaction, k ϩ is the rate constant of the forward reaction, and k Ϫ is the rate constant of the backward reaction. B, net flux and Gibbs free energy change of a reversible reaction. v is the net reaction rate, k is the rate constant, K is the equilibrium constant, ⌬G is the reaction Gibbs free energy change, R is the universal gas constant, and T is the temperature. Other variables are the same as in A.

Thermodynamics in control of metabolism
pathway is in steady state, in which all reactions have identical net flux. In other words, all intermediary metabolites have balanced fluxes feeding and consuming them. Thus, we can write n Ϫ 1 equations for the steady-state constraints.
Let J denote the net flux through this pathway. The flux control coefficient, C vi J , quantifies the sensitivity of the pathway flux J to perturbation in activity (i.e. rate constant) of the ith reaction. It is defined as the ratio of relative change in the flux J to relative change in the rate of the ith reaction when an arbitrary parameter p that only affects the rate of the ith reaction and no other reactions has a small change.
If p ϭ k i , we have the following.
The flux control coefficients can be uniquely determined by solving equations derived from the summation and connectivity theorems. For linear pathways, the summation and connectivity theorems do not explicitly include the pathway flux J.

Summation theorem
Connectivity theorem According to Equation 1.5, we have the following. A, diagram of a linear pathway and related parameters. S in is the input substrate, S out is the final product, S i is the ith intermediary metabolite, k i is the rate constant of the ith reaction, K i is the equilibrium constant of the ith reaction, and J is the pathway flux. B, scatter plots comparing flux control coefficients and reaction free energy changes in a linear pathway with randomly sampled parameters. The pathway includes 10 reaction steps. C vi J is the ith flux control coefficient, and g i ϭ ⌬G i /RT quantifies the ith free energy change. C, scatter plots comparing flux control coefficients and the thermodynamic driving force (max{g i }) in a linear pathway with randomly sampled parameters. Variables are the same as in B. D, relationship between free energy of a reaction and the elasticity coefficient of its velocity toward the substrate. E, relationship between free energy of a reaction and the elasticity coefficient of its velocity toward the product. F, relationship between the free energies of two sequential reactions and the ratio of the flux control coefficients associated with the two reactions.

Thermodynamics in control of metabolism (Eq. 2.6)
Flux control coefficients can be solved by combining Equations 2.1, 2.4, 2.5, and 2.6 (see supporting information). They can be written as either explicitly including only the reaction free energy terms or explicitly including only the rate constants and equilibrium constants.
According to Equations 2.7 and 2.8, the flux control coefficients are completely determined by kinetic and thermodynamic parameters of reactions in the pathway. Therefore, alterations in concentrations of the input substrate and the final product do not influence the flux control coefficients.
To illustrate the quantitative relationships defined in Equations 2.7 and 2.8, we consider a linear pathway consisting of 10 reactions; randomly sample 20,000 combinations of the parameters S in , S out , {k i }, and {K i }; and compute the corresponding flux control coefficients. We select the size 10 because it approximates the typical length of linear metabolic pathways. All parameters and substrate concentrations are sampled from a log-normal distribution.
We first correlate the flux control coefficient with reaction free energy change for each individual reaction (Fig. 2B). We find that the relationship between reaction free energy changes and flux control coefficients calls into question the longstanding hypothesis that reactions with the most negative free energy changes, such as the reactions catalyzed by hexokinase, phosphofructokinase, and pyruvate kinase in glycolysis, serve as rate-limiting steps of a pathway (35-38). Although previous work in metabolic control analysis has provided the theoretical framework to study the role of thermodynamics in the regulation of metabolic fluxes (34), this hypothesis is still widely accepted (39 -41). However, according to our analysis, for all reactions except the first reaction, their flux control coefficients correlate poorly with the reaction free energy changes, suggesting that a large absolute value of free energy change is dispensable for a rate-limiting step, and regulation of metabolic fluxes by an individual reaction is determined by both thermodynamics property and position along the pathway.
We next investigate how the flux control coefficients are associated with global thermodynamic properties of the entire pathway. We quantify the deviation from the thermodynamic equilibrium by the reaction free energy change that is closest to zero (i.e. the maximal value of free energy change among all reactions). When this quantity equals zero, the pathway has zero net flux (i.e. thermodynamic equilibrium). Moreover, when this quantity is close to zero, at least part of the pathway is not efficiently driven by thermodynamics, resulting in inefficient usage of enzymes with free energy changes close to zero (34). Thus, it is termed the thermodynamic driving force (TDF).
We correlate the thermodynamic driving force with flux control coefficient of each reaction step (Fig. 2C). From these simulations, a pattern emerges, in which the thermodynamic driving force determines either the upper bound or lower bound of the flux control coefficients; for the flux control coefficient of the first reaction (C vi J ), its lower bound depends on the thermodynamic driving force, whereas for the rest of the reactions, the upper bounds of their flux control coefficients depend on the thermodynamic driving force.
Analytical relationships that constrain flux control coefficients by thermodynamic driving force can also be derived from Equations 2.7 and 2.10 (see supporting information).
From Equations 2.11 and 2.12, we can derive corollary relationships for the flux control coefficients when the system is very far from equilibrium (i.e. thermodynamic driving force approaches minus infinity).
(Eq. 2.13) Equation 2.13 suggests that in a linear pathway far from thermodynamic equilibrium, the metabolic flux through the pathway is fully controlled by the enzyme catalyzing the first step. In other situations when the thermodynamic driving force has a smaller absolute value, the flux control is more evenly distributed among all enzymes in the pathway.
An intuitive interpretation of the relationship between the thermodynamic driving force and flux control coefficients can be given based on the fact that for a reaction with firstorder kinetics, the elasticity coefficients toward the sub-Thermodynamics in control of metabolism strate and product can be completely determined if the free energy of this reaction is known (Fig. 2 (D and E) and Equation 2.6). According to the connectivity theorem (Equation 2.5), in a linear pathway, the ratio of the ith flux control coefficient to the i ϩ 1th flux control coefficient can be written as a function of the free energies of the ith and i ϩ 1th reactions (Fig. 2F).
Thus, the thermodynamic driving force constrains the distribution of flux control among enzymes by constraining the elasticity coefficients that determine flux control coefficients through the connectivity theorem. When the thermodynamic driving force has a large absolute value (i.e. TDF approaches minus infinity; bottom left corner in Fig. 2F), for any reaction in the pathway, its elasticity coefficient with respect to the product is very close to zero, which means that altering the concentration of the product has negligible influence on the rate of this reaction. In this situation, perturbation of any downstream enzyme has a minimal effect on the first reaction, whose rate is equal to the overall flux, because the effect must be propagated through the metabolite S 1 , which serves as the product of the first reaction. On the other hand, when the thermodynamic driving force is close to zero (top right corner in Fig. 2F), the ratio between two sequential flux control coefficients is more flexible, resulting in what we term as a more adaptable pattern of flux control in these pathways.
To summarize, for a linear metabolic pathway with first-order kinetics, we have derived analytical relationships between the flux control coefficients and thermodynamic properties of both individual reactions and the entire pathway. We have shown that the flux control coefficients correlate poorly with free energy changes of the corresponding reactions but that the thermodynamic driving force of the whole pathway places bounds on how much any enzyme can control pathway flux. However, when the pathway is very far away from thermodynamic equilibrium, the metabolic flux through this pathway is completely controlled by the first reaction. We next investigate whether these principles are conserved in metabolic networks with more complex topologies.

Branch point with two upstream fluxes
We next consider a metabolic network with two converging fluxes, J 1 and J 2 , at a branch point (Fig. 3A). Through the two fluxes, the metabolite at the branch point, S BP , is produced by two input substrates, S in,1 and S in,2 . The final product S out is produced from S BP with the flux J 1 ϩ J 2 at the steady state. The three reactions included in this network have rate constants {k 1 , k 2 , k 3 } and equilibrium constants {K 1 , K 2 , K 3 }. At the steady state, the fluxes J 1 and J 2 and the concentration of S BP can be solved from the rate constants, equilibrium constants, and concentrations of the input substrates and end product based on the kinetic rules. We can then determine the flux control coefficients from the summation and the connectivity theorems. Here, we have three fluxes and three reactions in the network, giving nine flux control coefficients. It is worth noting that Equation 3.5 is not sufficient for determining all values of the flux control coefficients because it provides six linear equations for nine unknowns. To uniquely calculate every element of C v J , three additional equations from the connectivity theorem are needed.

Connectivity theorem
In Equation 3.7, ⑀ represents the matrix consisting of the normalized elasticity coefficients.  Thermodynamics in control of metabolism (Eq. 3.9) Furthermore, the reaction free energy changes can also be written as functions of the kinetic parameters and substrate concentrations.
(Eq. 3.10) We then combine Equations 3.9 and 3.10 to express the flux control coefficients in terms of the reaction parameters and free energy terms. To simplify the expressions, we let h i ϭ e gi . Then we have the following.

(Eq. 3.11)
To illustrate the relationships between flux control coefficients and thermodynamic properties of the network, according to Equations 3.9-3.11, we simulate the flux control coefficients and reaction free energy changes based on 20,000 combinations of parameters randomly sampled from a log-normal distribution as described previously in Equation 2.9. Combinations of parameters resulting in negative fluxes are dis-carded. The distributions of flux control coefficients (Fig. 3B) show that C v2 J1 , C v1 J2 , C v3 J1 , and C v3 J2 can have absolute values much larger than 1 under certain parameter combinations, indicating that the two upstream fluxes can be dramatically altered in response to perturbation of activities of enzymes in other branches. Based on the simulation, we also compare flux control coefficients with the reaction free energy changes (Fig. 3C) and the thermodynamic driving force (Fig. 3D). Consistent with the case of the linear pathway, most of the flux control coefficients correlate poorly with free energy changes of the individual reactions catalyzed by the corresponding enzyme (Fig. 3C) but exhibit a dependence on the thermodynamic driving force (Fig. 3D). Moreover, the quantitative relationships characterizing how flux control coefficients are constrained by the thermodynamic driving force can also be analytically derived (see the supporting information).

(Eq. 3.12)
From Equation 3.12, we can derive the limits of all the flux control coefficients except for C v1 J1 ϩJ2 and C v2 J1 ϩJ2 when the thermodynamic driving force goes to minus infinity. Thus, when the pathway is far from equilibrium, the two upstream fluxes J 1 and J 2 are fully regulated by activities of the Thermodynamics in control of metabolism enzymes that generate them. Moreover, the downstream flux J 1 ϩ J 2 is also fully controlled by the two upstream enzymes. The inequalities in Equation 3.12 also indicate that the feasible regions of the flux control coefficients C v2 J1 , C v1 J2 , C v3 J1 , and C v3 J2 become semi-infinite (i.e. only bounded in one direction) when the thermodynamic driving force approaches zero (i.e. there exists at least one near-equilibrium reaction in this pathway).
To summarize, in a metabolic pathway with a branch point involving two converging fluxes, regulation of the fluxes is largely dependent on how close to thermodynamic equilibrium the pathway in entirety is. When the pathway is far away from equilibrium, all three fluxes are completely controlled by the two upstream enzymes. On the other hand, when there exists at least one nearequilibrium reaction in the pathway, there is more flexibility in the regulation of fluxes, and, unlike in the case of linear pathways, the two upstream fluxes are able to be dramatically altered even by enzymes that are not directly involved in the reaction.

Branch point with two downstream fluxes
We next consider a pathway with one branch point in which one upstream flux, J 1 , diverges to two downstream fluxes, J 2 and J 3 , and show that the principles of metabolic flux regulation and thermodynamics also hold in this case. This network includes one input substrate, S in , and two final products, S out,1 and S out,2 , which are all connected to the branch point metabolite S BP . All kinetic rules and assumptions are as described above. Thus, at the steady state, concentration of S BP can be solved from the reaction parameters and concentrations of other metabolites in the network.  Here, we also have in total nine flux control coefficients that can be solved based on the summation theorem and connectivity theorem. To investigate the relationships between flux control coefficients and thermodynamic variables, we randomly sample 20,000 combinations of reaction parameters and substrate concentrations and calculate flux control coefficients and free energy changes corresponding to the randomly sampled parameters from Equations 4.3 and 4.6. Distributions of the sampled flux control coefficients suggest that the two downstream fluxes are able to be regulated in response to variations in the activity of enzymes not directly carrying them (i.e. C v1 J2 , C v3 J2 , C v1 J3 , and C v2 J3 Thermodynamics in control of metabolism can have absolute values larger than 1; Fig. 4B). Again, we observe a poor correlation between flux control coefficients and free energy change of the corresponding reaction steps (Fig. 4C). Among the nine flux control coefficients, only C v1 J2 ϩ J3 exhibits a clear dependence on the free energy change g 1 .
Despite the poor correlation between flux control coefficients and free energy changes of individual reaction steps, seven of the nine flux control coefficients are strictly con-strained by the thermodynamic driving force (Fig. 4D). This corroborates the findings in other types of network topologies that the regulation of metabolic fluxes is constrained by the deviation of the whole pathway from thermodynamic equilibrium. We also analytically derive the quantitative relationships between the upper and lower bounds of the flux control coefficients and the thermodynamic driving force of the pathway (see the supporting information).  Thus, based on a combination of analysis and simulation, we have demonstrated that, at branch points, the regulation of fluxes by enzyme activities is constrained by the thermodynamic driving force but has little correlation with free energy changes of the individual enzymes. These findings, together with similar results in other network structures, suggest that the influences of thermodynamics on regulation of metabolic fluxes by enzyme activities primarily occur at the pathway level rather than at the individual reaction level. These principles are applicable to simple metabolic network models regardless of the structure of the metabolic network, but they also rely on the assumption that all enzymes follow first-order kinetics. We next examine whether the conclusions are conserved in metabolic networks with Michaelis-Menten kinetics.

Thermodynamics and flux control in Michaelis-Menten models
We now consider zero-order kinetic models as an extreme case of substrate-saturable enzymes and repeat the theoretical analysis and simulation in these models to characterize the relationships between thermodynamics and flux control. Briefly, although the expressions of steady-state concentrations of the intermediary metabolites and steady-state fluxes are much more complicated with this assumption, the flux control coefficients can still be written in the form of simple functions of the reaction parameters and free energy changes. However, the thermodynamic driving force no longer constrains the flux control coefficients in this case (Figs. S1 and S2). Most of the flux control coefficients in all three network structures have the form of a zero-order homogeneous function of the thermodynamic terms h i . Consequently, moving all reactions away from or toward the thermodynamic equilibrium without influencing the ratios between the h i values does not change the resulting flux control coefficients.
Because flux control is no longer constrained by the thermodynamic driving force in zero-order kinetic models, we hypothesize that in Michaelis-Menten models, which lie between the extreme cases of first-order and zero-order kinetic models, the saturation of enzymes by metabolites determines whether the flux control coefficients are still constrained by the thermodynamic driving force as they are in first-order kinetic models. In Michaelis-Menten models (28), the rate of a reaction depends on the catalytic constant (i.e. turnover number) k and the equilibrium constant K of the enzyme, the concentrations of its substrate S and product P, and the Michaelis-Menten constants K S and K P for the substrate and product, respectively.
(Eq. 5.1) We note that for Michaelis-Menten models, the elasticity coefficients toward the substrate and product concentrations differ from those in first-order kinetic models only by one term quantifying saturation of the enzyme by the given metabolite (i.e. substrate or product).

Thermodynamics in control of metabolism
For Michaelis-Menten models, because the rate law is more complicated than that of zero-order or first-order kinetic models, we use sets of parameters randomly generated, as previously done for first-order kinetic models to simulate the thermodynamic variables and flux control coefficients, and then compare them with those in first-order kinetic models. To quantify the saturation of enzymes, we introduce the saturation term for a reaction (Fig. 5A). This saturation term considers the total saturation of an enzyme by the substrate and the product. For a pathway, the overall level of enzyme saturation is quantified by the maximal saturation term (Fig. 5A).
To evaluate whether the thermodynamic variables and flux control coefficients in a Michaelis-Menten model are consistent with the relationships we have derived for first-order kinetic models, we compute the deviation of the flux control coefficients in the Michaelis-Menten model from those allowed in the first-order kinetic model with the same topology and thermodynamic driving force (Fig. 5B).
We first examine how the maximal saturation term affects the fraction of Michaelis-Menten models consistent with firstorder kinetic models. In all three pathway structures, we observe a gradual increase in the fraction of Michaelis-Menten models consistent with the first-order kinetic models (Fig. 5C) as the maximal saturation term decreases, suggesting that the flux control coefficients are largely constrained by the thermodynamic driving force when the enzymes in the pathway are not highly saturated. We next compare the maximal saturation term with the deviation from the results in the first-order kinetic model. We find that the deviation gradually decreases as the maximal saturation term decreases, corroborating the findings that the relationships between thermodynamics and flux control are conserved in Michaelis-Menten models if the enzymes are not highly saturated.
To summarize, here we show that the principles for thermodynamics and flux control in metabolic pathways with firstorder kinetics are still effective in Michaelis-Menten models if the enzymes are not highly saturated. It is now reasonable to ask whether the same rules apply to more complicated metabolic networks in which most enzymes not only follow the Michaelis-Menten kinetics, but also have their activities regulated by long-range interactions such as allosteric feedback.

Thermodynamics in control of metabolism Thermodynamics and flux control in a mathematical model of glycolysis
Thus, we considered a model of glycolysis to study the relationships between thermodynamics and flux control (42). We have previously constructed a mathematical model of glycolysis and used this model to study the regulation of the lactate production flux (17,19). The topological structure of this model can be approximated by a linear pathway consisting of 13 reversible reactions (Fig. 6A). Elements and interactions such as reactions with more than one substrate, allosteric feedbacks, conserved species (e.g. NAD ϩ /NADH, ATP/ADP/AMP), and irreversible reactions are also included.
We generated 50,000 sets of random parameters to simulate the glycolysis model and computed the flux control coefficients and thermodynamic driving force at the steady state for each set of parameters. We also computed the maximal saturation term and assessed how the saturation of enzymes influences the relationships between thermodynamics and flux control. We first compared the thermodynamic driving force and flux control coefficients in the glycolysis model with the values allowed in the model for a linear pathway with the same number of reactions (i.e. 13 reactions) according to Equations 2.11 and 2.12 (Fig. 6B). As expected, we find that for parameter sets yielding a high (i.e. close to 1; red dots in Fig. 6B) maximal saturation term, the glycolysis model is inconsistent with the relationships derived under the assumption of first-order kinetics. We then estimated the deviation from first-order kinetic models (Fig.  6C) and computed the fraction of models consistent with firstorder kinetic models (Fig. 6D). Consistent with the Michaelis-Menten model, we found that as the deviation decreases, the fraction of consistent models increases, suggesting that the relationships derived in first-order kinetic models with simple topological structures are largely conserved in more complicated metabolic networks.

Discussion
In this study, we derive quantitative relationships between thermodynamics, enzyme activity, and regulation of metabolic fluxes. For a set of example pathways, we calculate the flux control coefficients as functions of the enzyme rate constants, reaction equilibrium constants, and Gibbs free energy. Based on numerical simulation and exact analytical results, we find that in all network topologies considered, the flux control coefficients are bounded by the thermodynamic driving force, as defined by the deviation of the entire pathway from thermodynamic equilibrium. Moreover, distinct patterns of flux regulation emerge when the thermodynamic driving force approaches two extreme values: if the thermodynamic driving force is very negative (i.e. all reactions are far away from equilibrium), enzymes catalyzing the upstream reactions (i.e. reac-

Thermodynamics in control of metabolism
tions directly consume the input substrates) exert full control of all fluxes; on the other hand, if the thermodynamic driving force is close to zero (i.e. near-equilibrium reactions exist in the pathway), there is more flexibility in the pattern of flux control, and fluxes in network topologies with branch points can be dramatically altered by enzymes not directly carrying the fluxes. These findings suggest that the coupling between thermodynamics and regulation of metabolic flux occurs at the systemic level and challenges the rule of thumb that the reaction with the most negative free energy change serves as the rate-limiting step.
It is also worth noting that all analysis that we have done to derive analytical expressions here relies on the assumption of first-order kinetics. This approximates the more generalized Michaelis-Menten mechanism when the abundance of substrate is far below the K m of the enzyme. Accordingly, in a study comparing substrate concentrations and K m values in different cell types, among all substrate-enzyme interactions investigated, around 70% exhibited higher substrate concentration than the K m value, thus requiring a Michaelis-Menten or substrate-saturated, zero-order kinetics (30). For Michaelis-Menten kinetics and a glycolysis model, we have shown, that the generalizations of our findings depend on the levels of substrate saturation. We also note that the tendency of an enzyme to be operating in the linear region or to be saturated by its substrate is highly dependent on the substrate used by this enzyme. For instance, kinases use ATP as the substrate whose physiological concentration is much higher than the typical K m value of enzymes and thus are more likely to be substrate-saturated. On the other hand, other enzymes, such as chromatinmodifying enzymes, use metabolites much less abundant as substrates (43). Nevertheless, because the relationships between thermodynamics and flux control are conserved under circumstances where the enzymes are moderately saturated, we expect that our findings are applicable to more complicated biochemical processes with Michaelis-Menten kinetics.
Finally, although the network structures that we study here are much simpler than real metabolic networks, the conclusions that we derive here can still be extended to more complicated network topologies. A metabolic network without cycles can always be simplified to a set of branch points connected by linear reaction chains with varying lengths, and the enzymes in the same linear reaction chain can be treated as an entirety to calculate the flux control coefficients with respect to simultaneous change in these enzymes (44). In this case, each linear reaction chain is simplified to a single reaction step, which enables us to apply the principles we found in the simple network structures with a branch point.
To summarize, we characterize the quantitative relationship between thermodynamics and regulation of metabolic fluxes by enzymes in a metabolic network in this study. We find that the global thermodynamic property of the network (i.e. thermodynamic driving force) constrains almost all flux control coefficients in both linear and branched pathways. Specifically, fluxes in metabolic networks far away from thermodynamic equilibrium are almost fully controlled by enzymes catalyzing the upstream reactions (i.e. reactions directly consuming the input substrates). On the other hand, near-equilibrium metabolic networks allow more flexibility in the pattern of regulation.
Semi-infinite feasible regions of flux control coefficients are only allowed in branched pathways with near-equilibrium reactions. These findings highlight the importance of global thermodynamic features in constraining the pattern of regulation of metabolism.

Experimental procedures
Analytical expressions of flux control coefficients were calculated using the function LinearSolve[] in Wolfram Mathematica version 11. For the computer simulations, 20,000 combinations of random parameters were generated using the function normrnd() in MATLAB R2016b for each model. Flux control coefficients and reaction free energies were then computed using the parameter sets that generate positive flux configurations. For the glycolysis model, a model file in SBML format was downloaded from the BioModels Database (https:// www.ebi.ac.uk/biomodels-main/) 4 (45) and translated to Cϩϩ codes using the SBML translator module in the Systems Biology Workbench (http://sbw.sourceforge.net/) 4 (46). Simulation of the model was done with the ODE solver gsl_odeiv2_ step_msbdf in the GNU Scientific Library (https://www. gnu.org/software/gsl/). 4 50,000 random sets of parameters were generated by sampling V max values of the enzymes with Latin hypercube sampling. The flux control coefficient of the glycolytic flux (defined as the rate of glucose consumption) with respect to a reaction step was estimated using the finite difference approximation, in which J glycolysis is the glycolytic flux and V max,i is the maximal velocity of the ith reaction. Free energy changes for the reactions were calculated from standard reaction free energies and steady-state concentrations of the metabolites. All source codes are available at the GitHub page of the Locasale laboratory (https://github.com/LocasaleLab/MCA_thermodynamics). 4


Prionospio from the coast of the Iberian Peninsula, with the description of two new species (Annelida, Spionidae)

Abstract Five species of Prionospio Malmgren, 1867, each with four pairs of branchiae, are studied from coast of the Iberian Peninsula. Two of these species, Prionospiocristaventralissp. n. and P.paraparisp. n., are new to science, whereas P.caspersi Laubier, 1962, P.fallax Söderström, 1920, and P.ehlersi Fauvel, 1928 have been previously recorded. Prionospiocristaventralissp. n. is characterized by having ventral crests present on chaetigers XI–XIX, dorsal crests and low crests on chaetigers X–XXXIV, triangular neuropodial postchaetal lamellae with pointed ventral edges on chaetiger II, oval neuropodial lamellae on chaetiger III, digitiform pinnules on the posterior face of the first and fourth pairs, and branchial pairs II and III are triangular. Prionospioparaparisp. n. is characterized by having rounded neuropodial postchaetal lamellae on chaetiger I, digitiform pinnules on the posterior face of the first and fourth pairs, branchial pairs II and III are cirriform, low dorsal crests on chaetigers VIII–IX, and oval neuropodial lamellae with enlarged dorsal edges on chaetiger III. A key is given to all Prionospio species with four pairs of branchiae known from the Iberian Peninsula coastline.

Introduction
Prionospio was established by Malmgren (1867) for P. steenstrupi Malmgren, 1867, a spionid species with branchiae on chaetigers II-V, the first and fourth pairs of which are pinnate, and the second and third pairs apinnate. With the discovery of additional species, the diagnosis of the genus was widened to include species with different branchial shapes and arrangements, with the chaetiger upon which the branchiae first arise being particularly relevant. Variability has thus increased, making Prionospio a very heterogeneous genus with about 100 species (Sigvaldadóttir 1998). As a result, several authors have suggested that some of these should be reclassified into new genera and/or subgenera (Foster 1971, Blake and Kudenov 1978, Maciolek 1985, Wilson 1990, Sigvaldadóttir 1998. So far Prionospio caspersi Laubier, 1962, P. fallax Södreström, 1920, and P. ehlersi Fauvel, 1928 have been recorded from the Iberian Peninsula: the former from Catalonia (Desbruyères et al. 1972) and from off Aveiro (Ravara and Moreira 2013); the second from several localities on the Portuguese coast (Quintino and Gentil 1987;Quintino et al. 1989;Dexter 1992;Pardal et al. 1992;Gil andSardá 1999, Ravara andMoreira 2013) and the third from the Bay of Biscay (Aguirrezabalaga and Cebeiro 2005), from Coruña (Amoureux 1972, López Jamar andGonzález 1986), from off Aveiro and off Porto (Amoureux 1974, Pardal et al. 1992, and from the South western continental shelf of Portugal (Gil 2011). In this study, we examined material deposited in the Museo Nacional de Ciencias Naturales de Madrid previously identified as P. caspersi and P. fallax.
Following the careful comparison between the redescription of P. fallax by Sigvaldadóttir and Mackie (1993), redescription of P. ehlersi by Mackie and Hartley (1990) as well as the original description of P. caspersi, we describe two new species: P. cristaventralis sp. n. and P. parapari sp. n. An identification key and a map with type localities are provided for all known Prionospio species with four branchial pairs from the coastline of the Iberian Peninsula (Fig. 1).

Materials and methods
The material examined belongs to the collections maintained at the Museo Nacional de Ciencias Naturales de Madrid, Spain (MNCNM) and the Alcalá de Henares University, Spain, as well as the personal collection of J Parapar at the University of La Coruña, Spain. Type specimens of the two newly described species are deposited in the MNCNM. For each species, the location where the specimens were collected is indicated in their respective sections.
All material was fixed in 10% formaldehyde in sea water and preserved and stored in 70% ethanol. In order to examine some morphological characters, specimens were dipped first into water and then into methyl green solution for staining. The color fades quickly when specimens are returned to the alcohol solution (Warren et al. 1994). Specimens were measured with a millimeter ruler: body widths refer to the maximum postbranchial body width (including parapodia but not chaetae) at about chaetiger VIII. The eye color mentioned in the descriptions of the species was based on preserved organisms.
Remarks. The specimens examined in this study agree with the original and subsequent descriptions of the species (Laubier 1962, Imajima 1990, Dagli and Çinar 2009. However, there is a slight difference between our specimens and Dagli and Çinar's specimens: the sabre chaetae were reported as present from chaetiger X by Dagli and Çinar (2009) whereas they first appear on chaetiger XI in our specimens. There are a few further differences between the specimens examined here and the description given by Imajima (1990). Imajima (1990) described the prostomium as being broadly flared anteriorly, with a slight medial indentation and unilimbate chaetae, whereas the specimens in this study have a narrower prostomium that is slightly truncate to convex anteriorly, with bilimbate chaetae. Due, to these morphological differences, the identity of the specimens recorded as P. caspersi from Japan by Imajima (1990)  Description. Incomplete specimens, 3.0-11.0 mm long, with 13-40 chaetigers, 0.6-1.0 mm wide. Posterior fragment 9.0 mm long for 26 chaetigers, 0.8 mm wide. Color in alcohol pale white. Some specimens with oocytes on chaetigers XXX-XXXIV.
Four pairs of branchiae present on chaetigers II-V, first pair longest and thickest (Fig. 2D); sometimes first and fourth pairs equal in length or fourth pair slightly longer than first. First pair with long, dense digitiform pinnules on lateral and posterior faces, continuing to tip; central stem of branchial pair 1 pinnate, cylindrical, very thick and with a blunt tip ( Fig. 2A, D). Pairs 2-4 apinnate; pairs 2 and 3 triangular, thick (Fig. 2D, E), slightly expanded distally, with rounded tips, densely ciliated laterally, shorter than notopodial lamellae and pinnate pair. Pair 4 cirriform (Fig. 2D), basally united by a short, moderate dorsal cord-shape (Fig. 2D); branchiae subequal to, or longer than, notopodial lamellae; some specimens with regenerating branchiae.
Remarks. Prionospio cf. ehlersi is very similar to the original and subsequent descriptions of P. ehlersi (Fauvel 1928, Mackie andHartley 1990), in that both describe the same prostomial shape, branchial arrangement, and hooded hook structure, and all have interparapodial pouches. However, specimens of this study differ from P. ehlersi in that the former have oval (rather than triangular to subtriangular) notopodial lamellae on chaetiger I, the neuropodial lamellae on chaetigers IV-V are oval (rather than rounded), and the second and third branchial pairs are triangular and thick (rather than expanded or swollen distally), a dorsal cord-shape (rather than low crest) on chaetiger V is present, and the sabre chaetae in P. cf. ehlersi are alimbate (rather than sheathed). Mackie and Hartley (1990) reported a variation in the shape of chaetiger I notopodial postchaetal lamellae, and so considered this unimportant. They also noted some parapodial variation around chaetiger XVIII-XX. They did not mention the oval anterior neuropodial postchaetal lamellae, however, Fauvel (1928, fig. 1b) shows a chaetiger IV with neuropodial lamella very similar to that in this study (Fig. 2J). We consider that these differences are important, but we consider it premature to erect a new species with these specimens, without being able to first compare them with Fauvel's syntype material. The syntypes are deposited at the Museum National d'Histoire Naturelle, Paris under MNHN A438, A449. However, the material cannot be sent abroad for re-examination.
Remarks. These specimens match the redescription given by Sigvaldadóttir and Mackie (1993), except that we found some specimens with eyes and others without eyes; one specimen had a single large brown eyespot. Possibly, the variation is due to the preservation of the specimens.
Pygidium missing. Remarks. Prionospio cristaventralis sp n. is similar to P. multicristata Hutchings & Rainer, 1979, P. nirripa Wilson, 1990, P. nielseni Hylleberg & Nateewathana, 1991, P. cornuta Hylleberg & Nateewathana, 1991, and P. paradisea Imajima, 1990 in that all show the same prostomial shape and neuropodial postchaetal lamellae on chaetigers II-III. However, P. cristaventralis sp. n. differs from these species due to the presence of ventral crests and dorsal crests limited to chaetigers X-XI, and low dorsal crests on chaetigers XII-XXII/XXXIV (end fragment) compared to low dorsal crests from chaetigers VII to XXV-XXX in P. multicristata, low dorsal crests from X-XIII to XX-VIII-XXXVII in P. nirripa, P. nielseni and P. cornuta and high crests from X to LX in P. paradisea. Prionospio cristaventralis sp. n. is also similar to P. pacifica Zhou & Li, 2009 in that both species have dorsal and ventral crests or folds. However, P. cristaventralis sp. n. can be distinguished from P. pacifica by having a prostomium that is rounded anteriorly (instead of being truncate), dorsal crests on chaetigers IV-V (instead of lacking such crests), and ventral crests on chaetigers XI-XIX (instead of only on chaetiger IX). The presence of ventral crests appears to be a unique feature of these two species.
Four pairs of branchiae present on chaetigers II-V (Fig. 4B, C); first pair longer than fourth pair, up to 1 time and a half longer than the fourth pair. Pairs 1 and 4 with long, slender, dense digitiform pinnules arranged along posterior face of the stems, pinnules thick, long, blunt in middle region of branchiae (Fig. 3B); branchiae with very long, naked, smooth, distal tips (Fig. 4B). Pairs 2 and 3 apinnate, cirriform, long, densely ciliated laterally, with pointed tips (Fig. 4E, F); subequal in length, up to 3 times shorter than pinnate pairs, but longer than notopodial lamellae (Fig. 4B, E, F).
Pygidium with one long median cirrus and two short, lateral lobes (Fig. 4X).
Remarks. Prionospio parapari sp. n. is very similar to P. fallax in having a prostomium that is truncated on the anterior margin, a neuropodial postchaetal lamellae on chaetiger II being the same shape, and a high dorsal crest on chaetiger 7. However, P. parapari sp. n. can be distinguished from P. fallax as redescribed by Sigvaldadóttir and Mackie (1993), by the former having rounded (rather than rectangular) postchaetal neuropodial lamellae on chaetiger I, the first pair of branchiae longer than fourth (rather than of equal size), and with dense digitiform (rather than sparse lateral) pinnules arranged along the posterior face of the stems on branchiae 1 and 4, and cirriform (rather than subtriangular) second and third branchial pairs. In addition, the low dorsal crests in P. parapari sp. n. are present on chaetigers VIII-IX while in P. fallax they are absent; the neuropodium on chaetiger III in P. parapari sp. n. is oval with an enlarged dorsal edge whilst in P. fallax it is subtriangular and dorsally pointed; and the sabre chaetae in P. parapari sp. n. are heavily granulated and bilimbate whereas in P. fallax they are distally granulated and alimbate. The pygidium lacks pigmentation in P. parapari.
Prionospio parapari sp. n. is also similar to P. komaeti Hylleberg & Nateewathana, 1991, P. depauperata Imajima, 1990, and P. oligopinnulata Delgado-Blas, 2015, in that all three species have a prostomium that is truncated on the anterior margin, the same shaped neuropodial postchaetal lamellae on chaetiger II, and a high membranous dorsal crest on chaetiger VII that decreases in height on chaetigers VIII-IX. However, P. parapari sp. n. differs from the first two species in that it has rounded (rather than square or lanceolate) notopodial and neuropodial postchaetal lamellae on chaetiger I, and an oval neuropodium on chaetiger III with the dorsal edge enlarged (rather than one that is square or triangular). In addition, the branchiae of P. parapari sp. n. have long, naked, smooth distal tips, whereas those of P. komaeti and P. depauperata have pinnules extending to the distal end, and P. parapari lacks low dorsal crests on chaetigers X-XI/XIII. Prionospio parapari sp. n. also differs from P. depauperata in that it has sabre chaetae without a distal filament, the notopodial hooded hooks start on chaetigers XIX-XXXVIII rather than XLV-XLVII, and the pygidium has two long lateral lobes rather than two short lateral cirri. Furthermore, P. parapari sp. n. is similar to P. oligopinnulata in that both species show the same pygidial structure, but differs in having cirriform rather than triangular second and third branchial pairs. In addition, the low dorsal crests on chaetigers X-XIV are absent in P. parapari sp. n., the neuropodium on chaetiger III is oval with the dorsal edge enlarged (rather than subtriangular and ventrally pointed), and the sabre chaetae are bilimbate (rather than alimbate). Prionospio parapari sp. n. is similar to P. rotunda Delgado-Blas, 2015 in that both species have large, rounded neuropodial postchaetal lamellae on chaetiger I, cirriform second and third branchial pairs, the first branchial pair always longer than fourth pair, and the same pygidium structure. However, P. parapari sp. n. differs from P. rotunda in having a bottle-shaped prostomium that is truncated on the anterior margin (rather than a pyriform and rounded one). In addition, the low dorsal crests in P. parapari sp. n. are present on chaetigers VIII-IX whereas in P. rotunda they are absent, and the neuropodium on chaetiger III is oval with an enlarged dorsal edge (rather than rounded). The differences between this new species and the other species examined are provided in the key.
Etymology. The species is named in honor of Dr Julio Parapar, in recognition of his major contribution to the study of polychaetes from Spanish coasts.
Type locality. Ría de Ferrol and the mouth of Piedras River, Huelva, Spain. Distribution. To date, this species has been recorded only on the Spanish Atlantic coast (Ria de Ferrol and the mouth of Piedras River, Huelva).