Process Engineering of Biopharmaceutical Production in Moss Bioreactors via Model-Based Description and Evaluation of Phytohormone Impact

The moss Physcomitrella is an interesting production host for recombinant biopharmaceuticals. Here we produced MFHR1, a synthetic complement regulator which has been proposed for the treatment of diseases associated to the complement system as part of human innate immunity. We studied the impact of different operation modes for the production process in 5 L stirred-tank photobioreactors. The total amount of recombinant protein was doubled by using fed-batch or batch compared to semi-continuous operation, although the maximum specific productivity (mg MFHR1/g FW) increased just by 35%. We proposed an unstructured kinetic model which fits accurately with the experimental data in batch and semi-continuous operation under autotrophic conditions with 2% CO2 enrichment. The model is able to predict recombinant protein production, nitrate uptake and biomass growth, which is useful for process control and optimization. We investigated strategies to further increase MFHR1 production. While mixotrophic and heterotrophic conditions decreased the MFHR1-specific productivity compared to autotrophic conditions, addition of the phytohormone auxin (NAA, 10 µM) to the medium enhanced it by 470% in shaken flasks and up to 230% and 260%, in batch and fed-batch bioreactors, respectively. Supporting this finding, the auxin-synthesis inhibitor L-kynurenine (100 µM) decreased MFHR1 production significantly by 110% and 580% at day 7 and 18, respectively. Expression analysis revealed that the MFHR1 transgene, driven by the Physcomitrella actin5 (PpAct5) promoter, was upregulated 16 h after NAA addition and remained enhanced over the whole process, whereas the auxin-responsive gene PpIAA1A was upregulated within the first 2 hours, indicating that the effect of auxin on PpAct5 promoter-driven expression is indirect. Furthermore, the day of NAA supplementation was crucial, leading to an up to 8-fold increase of MFHR1-specific productivity (0.82 mg MFHR1/g fresh weight, 150 mg accumulated over 7 days) compared to the productivity reported previously. Our findings are likely to be applicable to other plant-based expression systems to increase biopharmaceutical production and yields.


INTRODUCTION
The production of biopharmaceuticals, which include enzymes, vaccines, antibodies, growth factors, and hormones, is a complex task where small differences in production conditions can influence product quality, efficacy, safety and yield. As most commercial biopharmaceutical proteins are complex glycoproteins which cannot be produced in microorganisms, more than 50% are produced in mammalian cells (Tripathi and Shrivastava, 2019). Plant-based biopharmaceuticals have gained increasing attention as an alternative to mammalian cell systems in the last decade; however, only βglucocerebrosidase to treat Gaucher disease (Elelyso ® ) produced in carrot cell suspensions is currently available on the market (van Dussen et al., 2013).
One advantage of plant-based expression systems is biosafety due to lack of contamination with animal-borne viruses. Furthermore, there are technical advantages associated with up-stream processing: For instance, a rapid scale-up of production using Nicotiana benthamiana and a transient expression system is feasible within 2 months (Shoji et al., 2012;Capell et al., 2020), which would allow a quick response to address a public health crisis. During the current COVID-19 pandemic, it became evident that diagnostic reagents and vaccine production capacity was not sufficient to meet demand, and transient expression in plants could compensate this (Tusé et al., 2020). The plant-derived virus-like particles vaccine candidate for COVID-19 by Medicago completed phase I clinical trials with promising results and the phase II/III trials are ongoing Ward et al., 2021). There are many plant-based vaccines in phase I clinical trials (Shim et al., 2019), and some of them completed phases II and III such as avian (monovalent) and seasonal (quadrivalent) influenza vaccines produced in N. benthamiana (Ward et al., 2020).
Plant cell-suspension cultures in bioreactors have advantages over whole-plant systems, because biopharmaceuticals can be produced under controlled and reproducible conditions, while complying with Good Manufacturing Practice (GMP) requirements. Furthermore, downstream processing (DSP) may be simplified, thus reducing production time and costs (Huang and McDonald, 2009). Among other features, the high rate of homologous recombination observed in the nuclear DNA of somatic cells positioned the moss Physcomitrium patens (Physcomitrella) as an attractive platform for the production of recombinant therapeutic proteins (Decker et al., 2014;Reski et al., 2018;Bohlender et al., 2020).
Physcomitrella can be cultivated in suspensions in a differentiated, genetically stable, filamentous stage, the protonema tissue. This haploid tissue consists of two cell types: chloronema and caulonema. Chloronemal cells are characterized by cross-walls perpendicular to the growth axis of the filament and are rich in chloroplasts, while caulonemal cells are narrower with oblique cross-walls and fewer chloroplasts (Reski, 1998). Both cell types expand by tip growth with different rates, chloronemal cells divide every 24 h, caulonemal cells every 7 h (Jang and Dolan, 2011). The transition from chloronema to caulonema can be triggered among others by glucose or the phytohormone auxin (Thelander et al., 2005). Glucose is coupled to developmental progression of Physcomitrella via cyclin D (Lorenz et al., 2003), and auxin is a signaling molecule involved in many developmental processes in plants including mosses (Decker et al., 2006;Guillory and Bonhomme, 2021;Suzuki et al., 2021).
Protonemal tissue grows under autotrophic conditions in inorganic media. Furthermore, the bioprocess has been scaled up to 500 L in wave bag photo-bioreactors and stirred tank photobioreactors complying with GMP conditions (Reski et al., 2015;Decker and Reski, 2020). The first recombinant pharmaceutical protein produced in moss which completed clinical trial phase I is α-galactosidase (Repleva aGal; eleva GmbH) to treat Fabry disease (Hennermann et al., 2019). In addition, α-glucosidase to treat Pompe disease (Repleva GAA, eleva GmbH) completed preclinical studies (Hintze et al., 2020). Both products showed superior characteristics compared to the mammalian cell-based products and are catalogued as biobetters (Hennermann et al., 2019;Hintze et al., 2020).
Around 70% of biopharmaceuticals are produced in Chinese hamster ovary (CHO) cells, with high productivities ranging between 0.1 and 10 g/L (Tihanyi and Nyitray, 2020). One of the bottlenecks of plant-based production is the low yield of the recombinant protein (Schillberg and Finnern, 2021). Albeit the cell-specific protein productivity of plant cells is comparable to CHO cells, the larger size of plant cells results in much lower volumetric productivity (Havenith et al., 2014). Protein yields generally range from 0.01 to 200 mg/L (Xu et al., 2011). For example, antibodies produced in tobacco leaves by transient expression yielded up to 2 mg/g FW (Zischewski et al., 2016) and a peptide vaccine produced in chloroplasts from tobacco yielded up to 7 mg/g FW (Molina et al., 2004). However, most complex biopharmaceuticals cannot be targeted to the chloroplasts but should be targeted to the secretory pathway, where they undergo posttranslational modifications. Therefore, many products could not yield more than 100 μg/g FW yet (Schillberg et al., 2019;Schillberg and Finnern, 2021). Productivities up to 300 mg/L of an IgG monoclonal antibody have been reported in the moss Physcomitrella (Niederkrüger et al., 2014) or up to 100 µg/g FW moss of a synthetic complement regulator (Top et al., 2019).
Many efforts have been undertaken to increase the yield in different stages of the bioprocess, such as optimization of the gene sequence, search for suitable promoters and terminators, design of synthetic 5′and 3′unstranslated regions (UTRs), insertion of multiple copies of the transgene into the genome, targeting strategies, co-expression with protease inhibitors, or suppression of protease gene expression, among others (Rozov and Deineko, 2019). Culture conditions also influence protein stability and productivity, which can be enhanced by optimizing physical parameters such as light, pH, temperature, agitation speed, aeration or nutritional requirements, bioreactor design and operating mode (batch, fed-batch, semi-continuous, continuous, perfusion). For instance, protein productivity increased between 1.2 and 25-fold in plant cell suspension cultures by changing bioreactor operating modes (Huang et al., 2010;Huang and McDonald, 2012).
The need for optimization and better process control in the biopharmaceutical industry has increased during the last years due to the rise in demand for these drugs and more competitive markets (López-Meza et al., 2016). Mathematical models in biotechnology are a powerful tool for rational and straightforward process development thus saving time and resources. Among these are phenomenological models which derived mainly from microbial systems, and need understanding of the biochemical process and relationship between macroscopic variables such as biomass, substrate and product (Marconi et al., 2014a). These models are also widely applied to optimize processes for biofuels and high-value compounds produced by microalgae, such as carotenoids, phycobilins and fatty acids, e.g. eicosapentenoic acid and docosahexaenoic acid (DHA) (Lee et al., 2015). Dynamic models to predict and optimize recombinant protein production are less exploited in molecular farming (biopharmaceutical production), probably because transient platforms often lack reproducibility between batches, due to high variability in expression levels (Marconi et al., 2014a). A model based on a design-of-experiments approach was implemented to solve the uncertainty and predict the amount of recombinant protein produced in these platforms by identifying the main factors involved in recombinant protein productivity (Buyel and Fischer, 2012). Plant cell-suspension cultures can also benefit from mathematical modelling to optimize the process with less experimental effort, however, just a few studies have applied them (Marconi et al., 2014a;Marconi et al., 2014b). Nevertheless, so far there are no dynamic models for recombinant protein production in the moss bioreactor.
Here, we aimed at increasing the recombinant protein yield in moss bioreactors. As a case study, we produced the complex glycoprotein MFHR1, a synthetic complement regulator (Michelfelder et al., 2018), in Physcomitrella. The human complement system constitutes a crucial part of innate immunity and protects the body from invading pathogens (Merle et al., 2015). MFHR1 is a fusion protein, which combines the regulatory activity of human complement factor H (FH) and FH-related protein 1. Alongside moss-produced FH, it is considered a potential therapeutic agent against complementassociated diseases such as C3 glomerulopathies (C3G), atypical hemolytic uremic syndrome (aHUS), or viral infections where complement over-activation has been associated to the pathogenesis, such as COVID-19 or severe dengue fever (Michelfelder et al., 2017;Michelfelder et al., 2018;Top et al., 2019;Ruiz-Molina et al., 2021). In our present study, we produced two MFHR1 variants, differing in one amino acid at position 62 of FH equivalent to position 193 of MFHR1 (V62I) (Ruiz- Molina et al., 2021).
Due to a high variation in MFHR1 yields between transgenic moss lines, we evaluated if product yield positively correlates with mRNA expression levels and transgene copies integrated into the genome. Our findings suggest that multiple transgene copies are desirable to maximize protein yields in Physcomitrella. We explored the effect of process operating mode (batch, fed-batch, and semicontinuous) on MFHR1 production as a means for targeted changes of environmental conditions for the cells. Based on these experiments we propose a kinetic model of the moss bioreactor which can accurately predict recombinant protein production, biomass growth and nitrate uptake under batch and semi-continuous operation, and can be used to maximize protein and biomass productivity. This model contributes to a deeper understanding of a plant-based system and dynamic changes of the variables. Moreover, we explored the effect of the crucial plant growth regulator auxin on recombinant protein yield and its influence on actin gene expression. The kinetic model and the positive influence of auxin on biopharmaceutical production might be applied to other plant-based systems aiming at increasing biopharmaceutical yields.

Protonema Suspension Cultures
Previously established transgenic moss lines P1 and N-179 producing MFHR1 V62 and MFHR1 I62 variants, respectively, were used (Top et al., 2019;Ruiz-Molina et al., 2021). These variants are characterized by valine or isoleucine in position 193, corresponding to position 62 of FH, which lead to a difference in their activity (Ruiz- Molina et al., 2021). For simplification and better understanding of the results, P1 is called line A and N-179 line B.
Suspension cultures were initiated with a cell density around 100 mg dry weight (DW)/L in 100 ml shaken flasks (25 ml working volume, silicone sponge), agitated at 120 rpm. Moss lines were grown in Knop medium supplemented with microelements (Knop ME) according to Heck et al. (2021) 2-(N-morpholino)ethanesulfonic acid (MES, 0.1%) was added and the pH was set to 5.5. To determine the correlation between copy Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2022 | Volume 10 | Article 837965 number of integrated expression cassettes, mRNA expression and protein production, transgenic moss lines were propagated in suspension cultures during 1 year. Cell suspension cultures initiated with the same cell density were harvested at day 18.

Description of the Mathematical Model
The model developed is built up from a physiological level represented by formal kinetics for the specific growth rates. Material balances are employed on the macroscopic level to describe the dynamics of biomass (c x ), nitrate (c N ) and recombinant protein (c p ) concentrations. r X , r N , r p define the specific growth rate, the specific substrate uptake rate and the specific recombinant protein production rate, respectively (d −1 ).
For light limitation: r N 1 y X,N .r X (5) Equations 1-5 describe the specific growth and substrate uptake rates. One important factor in autotrophic organisms is the light. Limitation of nutrients and light have been widely reported in microalgae (Béchet et al., 2013;Lee et al., 2015). To construct this kinetic model, just nitrate and light were taken into account since these are the first to be limited in Physcomitrella cultures under our conditions. This co-limitation was described using the so-called minimum law or threshold model, in which the growth is affected by the most limited substrate Eq. 4 (Lee et al., 2015). The most common models for growth and substrate consumption correspond to Monod type (Monod, 1949), which introduced the concept of a growth-limiting substrate. Further development has been introduced by Jhon Pirt interpreting substrate uptake as enzymatic step combined with Michaelis-Menten kinetics and in a second step a yield with which the cell can support growth based on the substrate taken up. In this way also an equation for the substrate uptake rate (r N ) is delivered (Eqs 1, 2 for nitrate limitation and Eqs 3-5 for light limitation), where y X,N is the biomass yield coefficient (g biomass/g nitrate consumed), and k i , and k N are the saturation constants of the limiting substrates, light and nitrate, respectively (g/L). Nitrate used for non-growth associated purposes (maintenance term) was neglected. Monod-like equations have been used to model growth kinetics in autotrophic processes, where light, the energy source, is considered a substrate (Béchet et al., 2013). The growth kinetics of Physcomitrella under different light intensities suggested Han Model behavior (Cerff and Posten, 2012b;Schediwy et al., 2019); therefore, we incorporated Monod's equation in the threshold model to explain the substrate consumption and specific growth rate (Eq. (3)).
Modelling of microalgae processes generally uses the Beer-Lambert law to predict the light distribution in the bioreactor (Béchet et al., 2013). This equation describes an exponential decay of light intensity from the external surface, which is partially compensated by the geometric intensification against the axis of the cylinder. This results in a reasonable constant light intensity in the reactor. Beer-Lambert assumes that the cells do not scatter the light, which does not reflect reality. Although more complex models exist, they are difficult to implement and are computationally intensive (Dauchet et al., 2013), however, the Beer-Lambert law could be accurate enough to describe light distribution in the bioreactor (Luo and Al-Dahhan, 2012). Therefore, we implemented this equation in our model, assuming that the reactor surface is homogenously illuminated, following the methodology described by Evers (1991)  The nitrogen source is essential for protein production, therefore we used Eq. 5 to simulate recombinant protein production originated from Monod model. r p,max represents the maximum specific recombinant protein production rate (d −1 ), k p is nitrate saturation constant for protein of interest synthesis (g/L). The degradation rate (r p,d ) (d −1 ) was also taken into account, depending on biomass and product concentration (Eq. 6).
The Eqs 7-10 describe a batch, fed-batch or semi-continuous process, where D is the dilution rate, and dV/dt is the dynamics of volume in the bioreactor. D is zero in a batch run and dV/dt is zero in semi-continuous and batch operations. c N,f is the substrate concentration in the feeding stream in a fed-batch or semicontinuous operation.

Parameter Estimation
Ordinary differential equations (ODE) were solved by the numeric method Runge-Kutta, using the solver solve_IVP implemented in SciPy (Python) (Virtanen et al., 2020).
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2022 | Volume 10 | Article 837965 Estimation of parameters was performed by minimizing the distance between experimental data (Y ij ) and the ODE solution (function (X i , P k )). The minimization of the objective or error function (least square minimum) (Eq. (12)) was performed by using differential evolution, a global optimization method implemented in SciPy (Python).

Bioreactor Cultures
The cultures were scaled up to stirred tank photobioreactor (5 L) with Knop ME medium using the following conditions: aeration: 0.3 vvm (2% CO 2 ), agitation with pitched 3 blade impeller at 500 rpm under continuous light. The light was increased from 160 to 350 μmol/m 2 s after 2 days. The illumination using LEDs were evenly distributed on the surface of the bioreactor, the pH was controlled at 4.5 and temperature was kept at 22°C. Light intensity (µmol m −2 s −1 ) was measured and calibrated with a light meter (LI-250 A, LI-COR) at the inner wall of the bioreactor before sterilization. For fed-batch operation the process was carried out in batch with 4 L volume for 7 days and the feeding was started at day 7. The feeding flow (F) with 5× concentrated Knop ME medium was set to 200 ml/d.
For the semi-continuous operation, the process was carried out in batch with 5 L volume for 5 or 6 days, respectively, and the semi-continuous phase started at day 6 or 7. For this, 2 L or 1 L suspension were harvested and replaced with fresh Knop ME medium daily, which correspond to a dilution rate of 0.4 or 0.2 d −1 , respectively.

Auxin Supplementation and Auxin Biosynthesis Inhibition
For the supplementation with the auxin 1-Naphtalene acetic acid (NAA, Sigma), 100 mM stock solution was prepared in 0.5 M KOH. The auxin biosynthesis inhibitor, L-kynurenine, (Sigma) was prepared in DMSO.

Determination of Biomass dry Weight, Nitrate Concentration and Protein Accumulation
For dry weight (DW) measurement, 10-50 ml of tissue suspension were vacuum-filtered and dried for 2 h at 105°C. For experiments at shaken-flask scale, growth index was calculated as [(maximum biomass-initial biomass)/initial biomass]. Nitrate concentration in the culture supernatant was measured in the reflectometer (RQflex ® ) using the Reflectoquant ® Nitrate Test following the manufacturer's instructions. During bioreactor runs, nitrate was also measured at day 0, because it differs slightly between batches, due to experimental error.

Nucleic Acids Extraction and Quantitative PCR
Auxin-treated samples were harvested from the bioreactor. The suspension was vacuum-filtered and 30 mg FW plant material frozen in liquid nitrogen. The frozen tissue was homogenized and disrupted using two beads (glass and metal) in a tissue-lyzer for 1 min at 30 Hz. RNA was extracted using RNeasy plant mini kit (Qiagen). After DNAse I treatment, cDNA synthesis was carried out from 1 µg RNA using Multiscribe RT and TaqMan Reverse Transcription reagents (Thermo Fisher Scientific) and a control of complete DNA digestion was included without addition of Multiscribe RT. Quantitative RT-PCR was performed with 10 ng/ μL cDNA, 400 nM of each primer and 1× SensiFast SYBR ® No-ROX (Bioline) in the LightCycler 480 System (Roche) under the following conditions: 95°C for 2 min (one cycle); 95°C for 5 s, 60°C for 10 s, and 72°C for 15 s (45 cycles). A thermal ramping stage was included to obtain melting curves.
Primers were designed using the Universal Probe Library Assay Design Center (Roche) and checked for specificity against the Physcomitrella transcriptome in the Phytozome database (Goodstein et al., 2012) release v13. Every primer pair was tested for specificity and efficiency using serial dilutions DNA. Primers with an efficiency lower than 1.9 were discarded. Expression levels of the transgene (MFHR1), PpIAA1A (Pp3c8_14720V3.1), and actin 7 (PpAct7a) (Pp3c3_33410V3.1), were analyzed using the primers listed in Supplementary  Table S1.
Normalized expression was calculated according to the ΔΔC T method, using the Light Cycler 480 software (Roche). Genes coding for elongation factor 1-α (EF1-α) (Pp3c2_10310V3.1) and 60S ribosomal protein L21 (Pp3c13_2360.V3.1) were used as housekeeping references to normalize the data. Relative quantification was performed using a time point before addition of NAA in the bioreactor as calibrator. Every sample was tested in triplicates.
The number of integrations of the MFHR1 transgene in the genome was estimated by qPCR as described above. DNA was Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2022 | Volume 10 | Article 837965 isolated using GeneJet Plant Genomic DNA Purification Mini Kit (Thermo Scientific) according to the manufacturer's instructions. The single copy gene PpCLF (Pp3c22_22940V3.1) was used as a reference to normalize the data using two primer pairs listed in Supplementary Table S1. The number of transgene copies was determined using two controls as calibrators, a moss line with a single integration of the hygromycin resistance gene (hpt), and the parental line (Δxt/ft) which has a single integration of the 5′HR region included in the construct (Ruiz- Molina et al., 2021). The MFHR1 transgene was amplified with the same primers used for the qRT-PCR. All primers are listed in Supplementary  Table S1.

Statistics
Analyses and figures were performed with the GraphPad Prism software version 8.0 for Windows (GraphPad software, San Diego, California, United States). Statistical significance was evaluated by ANOVA followed by Bonferroni or Tukey test (p < 0.05).

RESULTS AND DISCUSSION
A High Copy Number of Integrated Expression Cassettes Is Important for High MFHR1 Productivity First, we studied the correlation between copy number of integrated expression cassettes, mRNA expression and protein production for 14 MFHR1-producing moss lines. We analyzed the transgenic plants after a long period of cultivation (1 year after transfection) to avoid episomal expression (Murén et al., 2009). DNA-copy number, mRNA levels, and recombinant protein yields were analyzed by qPCR, qRT-PCR and ELISA, respectively. It is important to mention that under our conditions, MFHR1 accumulated in the cellular fraction, and was very low or not detected in the medium, therefore the MFHR1 concentrations described were measured from cellular extracts in the whole study. Our results reveal a positive correlation across the three molecular levels DNA-copy number, mRNA level and protein yield (Figure 1). According to previous studies protein and mRNA levels do not follow a normal distribution; therefore, the Spearman correlation was considered more suitable than the Pearson correlation (Maier et al., 2009). Accordingly, in our study mRNA expression and copy number of integrated expression cassettes could explain around 90% of the MFHR1 level variation (Figure 1). Many regulatory mechanisms could explain the remaining 10% such as DNA positional effects, methylation, miRNA and protein stability (Khraiwesh et al., 2010;Vogel and Marcotte, 2012;Myhre et al., 2013;Payne, 2015).
Although the correlation between protein and mRNA is weak or insignificant when groups of proteins or the proteome are analyzed (R 2 0.3-0.89) (Maier et al., 2009;Schwanhäusser et al., 2011), a better correlation was reported for individual proteins which are essential for growth, such as those implied in energy metabolism (Nie et al., 2006). This is congruent with our results, since MFHR1 expression is driven by a native actin gene promoter (PpAct5) (Top et al., 2019). Actins, essential components of the cytoskeleton, are ubiquitous proteins in eukaryotic cells and play important roles in cell division and extension, among others (Szymanski and Staiger, 2018;Wu and Bezanilla, 2018).
For further analyses, we selected the lines with the highest productivity for each MFHR1 variant, referred from now on as line A (producing MFHR1 V62 ) and line B (producing MFHR1 I62 ). At shaken-flask scale, line A produces approx. 20-30 µg MFHR1/ g fresh weight (FW) and is further referred to as the highproducer line, while line B shows a lower productivity of up to 10 µg MFHR1/g FW. Recombinant protein production was scaled up to 5 L stirred tank bioreactors, and batch, fed-batch and semi-continuous operation modes were tested in order to increase protein productivity and yield. For the sake of clarity, in this study we referred to specific productivity as the amount of recombinant protein per gram fresh weight biomass, recombinant protein specific production rate as g recombinant protein/(g biomass d), while yield is the mass of accumulated recombinant protein during the whole process.

Mathematical Modelling and Effect of the Bioreactor Operating Mode on Recombinant Protein Productivity
Bioreactor operation mode is one of the variables that needs to be optimized to increase product yield and to guarantee the reproducibility of the results. The choice of the operation mode in plants is host species-specific and depends on the FIGURE 1 | Positive correlation between three molecular levels: number of integrations of expression cassettes, mRNA expression level and recombinant protein production. The correlation was estimated with 14-15 MFHR1 producer lines, which include low and high producers. The experiment was performed 1 year after moss transfection. Linear regression in logarithmic scale is shown. Lines A and B selected for the next experiments are marked with arrows. mRNA expression levels were determined by qRT-PCR. Samples were normalized using two reference genes (coding for EF1-α and 60S ribosomal protein L21) and the line with the lowest expression level was used as calibrator. Number of integrations of expression cassettes was determined by qPCR. Samples were normalized with the single copy gene PpCLF (Pp3c22_22940V3.1).
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2022 | Volume 10 | Article 837965 expression system (constitutive or inducible promoter, intracellular or secreted protein) (Huang and McDonald, 2012). We evaluated different operation modes aiming at increasing the specific MFHR1 productivity to facilitate the purification of the protein from the cellular fraction. Furthermore, we used an unstructured kinetic model based on Eqs 1-11 described in the materials and methods section, which considers just the macroscopic variables, and tested its accuracy to predict recombinant protein production in moss suspension cultures in photobioreactors operated in batch, fed-batch and semi-continuous modes, respectively. The growth kinetic model, which describes the dynamics of biomass growth, MHR1 accumulation and nitrate consumption, was implemented with line B. The kinetic parameters were estimated using differential evolution algorithm implemented in SciPy (Python) ( Table 1).
The batch run took 8 days to reach stationary phase and MFHR1 ₋ specific productivity reached~65 ± 11.8 μg/g FW ( Table 2) at the end of the exponential growth phase and decreased steeply afterwards. Furthermore, growth cessation coincided with the complete consumption of nitrate, therefore we considered it in addition to light as a limiting substrate in the dynamic model.
Initially, the parental line (Δxt/ft) was cultivated in batch operating mode and the proposed model was fitted to experimental data to predict biomass growth and nitrate uptake rates. Kinetic parameters of the parental line ( Table 1, Supplementary Figure S2) allowed us to narrow the bounds to estimate the parameters of the moss transgenic line B. The maximum specific growth rate (r X,max ) of the parental line was 0.699 d −1 (doubling time approx. 1 day) and the biomass yield was 9.6 g biomass/g nitrate (Table 1). Our model fits accurately with the experimental data in batch operation until the beginning of the stationary phase, and it is able to predict concentration of product, nitrate and biomass (Figure 2A). The maximum specific growth rate of MFHR1-producing line B was almost identical to the parental line, but biomass yield was lower with 7.14 g biomass/g nitrate and maximal specific MFHR1 production rate (r p,max ) was 0.32 mg protein/g biomass d ( Table 1).
To achieve a higher biomass and subsequent higher MFHR1 yield, we fed our culture with nitrate (fed-batch). The kinetic parameters obtained in the batch mode were initially used to calculate the feeding flow rate (F) maximizing biomass productivity and F was set at 200 ml/d with a measured nitrate concentration of 2,500 mg NO 3 − /L. The bioreactor was initiated with 4 L, fed-batch operation started at day 7 and the process was ended at day 12. The experimental data obtained were used to estimate new parameters to fit the fed-batch mode accurately. The proposed model was able to predict the dynamic behavior of biomass and nitrate but not recombinant protein. We observed again a drop in MFHR1 production between day 9 and 10 before reaching stationary phase of growth although nitrate was not limiting in this case ( Figure 2B).
The discrepancy between model prediction and experimental data might be due to the age of the tissue and a change in the metabolic activities of the cells. After approx. nine days, the extracellular medium started turning slightly brownish, which might be attributed to accumulation of polyphenolic compounds and oxidative stress (Halliwell, 2003;Wilken and Nikolov, 2012). Although maximum biomass concentration increased from 3 g dry weight (DW)/L to 4 g DW/L compared to batch bioreactor, specific protein productivity and total product in the cellular fraction was similar (around 77.6 ± 9.1 µg MFHR1/g FW, Table 2). Kinetic parameters are not listed in Table 1, due to the discrepancy between model prediction and experimental data.
The semi-continuous bioreactor was run for 11 days under the same conditions used for the batch mode. From day 6-11, 2 L suspension culture were replaced daily with the same volume of fresh medium, which corresponds to a dilution rate (D) of 0.4 d −1 . However, as this D was slightly higher than the apparent specific growth rate (0.39 d −1 , calculated between day 6-11), washing out of the cells was likely to occur. Kinetic parameters are listed in Table 1. The proposed model was able to accurately predict biomass growth, nitrate uptake and protein production ( Figure 2C). The specific protein productivity was the lowest TABLE 1 | Estimated kinetic parameters for batch, fed-batch and semicontinuous model of an MFHR1 producer line (Line B, MFHR1 I62 ). r X,max : maximum specific growth rate. k N and k i : saturation constants of the limiting substrates, nitrate and light, respectively. y X,N : biomass yield coefficient. r p,max : maximum specific recombinant protein production rate. k p : nitrate saturation constant for recombinant protein synthesis. k pd : specific protein degradation rate. σ X : cell absorption cross section. r R : bioreactor cylinder radius. l 0 : incident light. NA: not-applicable.  compared with fed-batch and batch (around 50 µg MFHR1/g FW). In total, around 10 and 13 mg recombinant protein accumulated in the cellular fraction after 8 days in batch and fed-batch, respectively, while for the bioreactor in semicontinuous mode around 6 mg were produced within 11 days, including biomass removed during this operation ( Table 2). This operation could not be prolonged beyond 2 weeks, because the moss filaments formed big aggregates (pellets). This resulted in a reduced specific growth rate, making it difficult to reach a stable state (Supplementary Figure S3). Growth in pellets was also observed at the end of batch or fed-batch operation mode. Dense pellet formation could cause limitations in light, CO 2 and nutrient transfer, which can affect product formation. However, limitations in CO 2 could not be demonstrated in moss pellets (Cerff and Posten, 2012b). Formation of large entanglements and pellets can be avoided or delayed with low initial light intensities, however it reduces the maximum specific growth rate (Cerff and Posten, 2012b). The effect on recombinant protein production has to be further studied. The kinetic model was validated to assess its predictability under different conditions in the semi-continuous operation mode, using the same estimated kinetic parameters ( Table 1). For this the dilution rate and the starting time of the semicontinuous operation were modified to allow higher biomass density to accumulate and prevent washing out of the cells. From day 7 to day 10, 1 L of the suspension culture was replaced daily with the same volume of fresh medium (D = 0.2 d −1 ). Under these operation conditions, around 8.2 mg MFHR1 accumulated in the cellular fraction, 36% more than with a higher dilution rate (D = 0.4 d −1 ). The model fits well with experimental values for biomass, nitrate and recombinant protein concentration with a coefficient of determination (R 2 ) around 0.96, 0.99 and 0.88, respectively (Supplementary Figure S2B). According to the results, the model shows an acceptable predictability using several operating conditions with different variables.
Altogether, MFHR1 production is growth-associated, where both, a higher growth rate and a high cell density are crucial to increase the productivity. The specific MFHR1 productivity and total MFHR1 produced was affected by the change in medium conditions induced by the operating mode. The total amount of recombinant protein was doubled by using fed-batch or batch compared to semi-continuous operation (D = 0.4 d −1 ), although the maximum specific productivity increased by just 35%. Here, we propose an unstructured and nonsegregated kinetic model of the moss bioreactor to produce recombinant proteins, in this case MFHR1, which can be used as a starting point to optimize biomass and recombinant protein production under different operation modes or for the bioprocess design. The model fits accurately with the experimental data in batch and semi-continuous operation, and is able to predict recombinant protein production, nitrate uptake and biomass growth.

Effect of Auxin on the Production of MFHR1 at Shaken-Flask Scale
Auxin plays an important role in many plant processes as in protonema development, promoting the transition from FIGURE 2 | Mathematical modelling and simulation to describe the growth kinetics of the transgenic moss line B and the production of the recombinant protein MFHR1. Parameters were estimated using differential evolution (DE) implemented in Scipy in Python 3.6. The solid lines represent modelled data; experimental data are shown with symbols (A) Batch mode, (B) Fed-batch mode, (C) Semi-continuous mode. This operation was conducted between day 6 and 11. Data represent mean ± standard deviations (SD) from two measurements. The starting day of fed-batch and semi-continuous operations are marked by arrows.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2022 | Volume 10 | Article 837965 chloronema to caulonema cells (Decker et al., 2006). Exogenous auxin addition to the moss bioreactor triggered a temporal increase in MFHR1 production of 25% within 24 h (Top et al., 2019). Therefore, the effect of auxin on recombinant protein production was further studied in order to establish a strategy to increase protein yield. First we studied the effect of the synthetic auxin naphthalene acetic acid (NAA, 10 µM) on MFHR1-specific productivity in the highproducer line A compared to the effect of the auxin biosynthesis inhibitor L-kynurenine (L-Kyn, 10 and 100 µM) (He et al., 2011), and their respective solvent controls (0.00005% KOH and 0.01% DMSO), in shaken-flask scale for 7 and 18 days.
Addition of auxin (10 µM NAA) increased MFHR1 concentration by 470% at day 7 compared to the controls ( Figure 3A). In contrast, the specific protein production remained constant with 10 μM L-Kyn and decreased with 100 μM L-Kyn by 110% and 580% at day 7 and 18, respectively, while in the controls MFHR1 concentration doubled between day 7 and 18 ( Figure 3A).
The effects of auxin were also reflected by the morphological appearance of the tissues ( Figure 3B). Compared to the control, NAA-treated suspensions showed accelerated caulonema development at day 7. In contrast, the addition of L-Kyn (10 and 100 µM) slowed down caulonema formation. Furthermore, the higher concentration of L-Kyn led to the formation of round, thick-walled brachycytes ( Figure 3B), which are formed under stress or the addition of ABA (Decker et al., 2006). This observation supports a reduction of endogenous auxin levels by L-Kyn. Next, we studied the effect of different NAA concentrations (2-50 µM) at shaken-flask scale for 14 days. We did not find significant differences in specific recombinant protein productivity, although 50 µM NAA seemed to negatively affect protein production compared to 2-10 µM NAA (Supplementary Figure S4A). Likewise, we could not find significant differences in the growth index ((maximum biomassinitial biomass)/initial biomass) upon different auxin concentration treatments at shaken-flask scale (Supplementary Figure S4B).
Auxin responses normally follow a bell-shaped dose-response curve. At low concentrations and up to an optimum, auxin induces cell division and elongation; beyond this maximum, opposite effects were observed in Arabidopsis shoots and cell suspension cultures of tobacco BY2, when inhibition of growth and bundling of actin occurs (Huang et al., 2017). According to previous studies, usually for high NAA concentrations (>10 µM) the response is less pronounced (Waller et al., 2002;Huang et al., 2017). Therefore, we selected the concentration 10 µM for further studies in the bioreactor.
Our results confirm that exogenous auxin significantly enhances recombinant protein production in Physcomitrella (p < 0.0001) and the inhibition of the endogenous synthesis of auxin by L-Kyn decreases the productivity of the protein of interest. Two different mechanisms can play a role in the MFHR1 production response to NAA treatment: tissue differentiation to the faster-growing caulonema and an auxin signaling-dependent actin response, as the actin5 promoter is driving the expression of our recombinant product (Top et al., 2019). The transition from chloronema to caulonema can be mediated by auxin or the quality of energy supply (Thelander et al., 2005;Decker et al., 2006). Caulonemal cells exhibit a faster elongation rate than chloronema (~20 μm/h vs. 6 μm/h) (Menand et al., 2007). Optimal rates of tip growth are linked to this transition, which is characterized by non-stabilized actin cytoskeleton where actin genes are highly expressed. Furthermore, analysis of the transcriptome of caulonema and chloronema showed that all the processes related to tip growth are more active in caulonemal cells, such as cell-wall modification and regulation of cell size. An elongation of caulonemal cells can also result in enriched expression of certain actin genes (Ortiz-Ramírez et al., 2016). In addition, endoreduplication of the nuclear DNA was observed when caulonemal cells become older, and exogenous auxin led to a higher proportion of polyploid cells (Schween et al., 2003). Therefore, endopolyploidization can also contribute to the high specific recombinant protein productivity achieved upon auxin treatment.

Effect of Sugar on the Production of MFHR1 at Shaken-Flask Scale
The quality of energy supply can influence tissue differentiation in Physcomitrella, e.g., exogenous glucose favors caulonema development. However, these cells are different from those induced by auxin, since they are shorter and "heavily pigmented" (Thelander et al., 2005). Therefore, we evaluated the influence of mixotrophic and heterotrophic conditions on the productivity of MFHR1.
Transgenic moss line A was cultivated under autotrophic, mixotrophic (1% sucrose and 50-70 µmol/m 2 s light, photoperiod 16/8), and heterotrophic conditions (1% sucrose and darkness). Protonema suspension cultures were initially adapted to each condition during 2 weeks before the experiment. Interestingly, sucrose had a negative impact on recombinant protein production. The specific protein production at day 7 was similar between autotrophic and mixotrophic conditions, however, it decreased to zero at day 18 under the latter ( Figure 4A). Furthermore, recombinant protein concentration was almost negligible under heterotrophic conditions.
To analyze if the effect of sugar was dependent on the sugar type, we also tested 1% glucose for 15 days. Specific protein production decreased during the whole kinetics under mixotrophic as opposed to autotrophic conditions where protein specific productivity was~15-fold the productivity achieved with 1% glucose at day 15 ( Figure 4B). There was no significant difference in terms of growth between medium supplemented with 1% glucose and the standard medium during 15 days at shaken-flask scale and nitrate was not limiting (Supplementary Figure S5).

Effect of Auxin on the Production of MFHR1 in 5 L Bioreactor Under Different Operating Conditions
Subsequently, the effect of auxin supplementation was studied in 5 L stirred tank bioreactors operated in different conditions. Transgenic moss lines A and B were cultivated semicontinuously under the same conditions described above and 10 µM NAA were added daily from day 6 to keep the concentration stable ( Figure 5). To explore a putative auxin influence on actin gene expression, a kinetic expression profile of the MFHR1 transgene upon auxin treatment was generated by quantitative real-time PCR (qRT-PCR) from NAA-treated bioreactor samples at different time points. The expression of MFHR1 was compared to the auxin-responsive gene PpIAA1A and the actin gene PpAct7.
Aux/IAA proteins have an essential role in auxin response and gene regulation. Physcomitrella has three genes (PpIAA1A, PpIAA1B, PpIAA2) encoding these proteins, and a triple knockout of the genes rendered the plants completely insensitive to auxin (Lavy et al., 2016). Therefore, PpIAA1A (Pp3c8_14720V3.1), which is induced upon auxin treatment (Prigge et al., 2010), was chosen as a positive control during NAA treatment in the bioreactor.
As the expression of MFHR1 is driven by a PpAct5 promoter, we aimed to test the response of MFHR1 and actin gene expression, respectively, to auxin treatment. However, the Physcomitrella genome contains a duplication (Pp3c10_17070V3.1, PpAct5b) of PpAct5 (Pp3c10_17080V3.1, PpAct5a) with a different promoter sequence but almost identical coding sequences. Therefore, PpAct5 mRNA was not eligible for specific amplification. As the PpAct5a UTRs are present in our transgene construct, these sequences could not be used as a target for the analysis of PpAct5a expression in MFHR1 expressing lines either. A multigene family with different expression patterns and functions encodes actin proteins. P. patens has eight actin genes with highly conserved coding  regions. Although all actin genes are expressed in protonema, PpAct5a, PpAct7a (Pp3c3_33410V3.1), and PpAct7b (Pp3c3_33440V3.1) are highly expressed (Supplementary Figure S6), according to the expression data obtained from PEATmoss (Fernandez-Pozo et al., 2020). Moreover, according to a distance-based dendrogram (Supplementary Figure S7), PpAct5 and PpAct7 are closely related and the expression levels of PpAct5a and PpAct7a are similar and higher than PpAct7b (previously called PpAct3) (Weise et al., 2006). Therefore, the expression level of PpAct7a was evaluated in the MFHR1expressing moss line under auxin treatment.
PpAct7a expression was relatively stable during the bioreactor run suggesting that this gene is not regulated by auxin ( Figure 5A), as opposed to Arabidopsis thaliana Act7, which is the only Arabidopsis actin gene strongly responding to auxin (McDowell et al., 1996). However, the Aux/IAA gene and MFHR1 transgene were upregulated during the whole kinetics ( Figure 5A). Interestingly, the MFHR1 transgene was upregulated during the first 64 h in a similar way to the Aux/ IAA gene, whilst a further increase in transgene expression compared to Aux/IAA was observed after 72 h when the semicontinuous operation was stopped for 2 days to let biomass increase. Thus, MFHR1 RNA level-increase coincided with increased biomass and MFHR1 protein accumulation ( Figures 5A,B).
The transgene was upregulated within 16 h after NAA addition ( Figure 5A). To further unravel the kinetic response of MFHR1 to NAA, the expression level was analyzed within 2-23 h after NAA addition, using the MFHR1 high-producer line A ( Figures 5C,D). The expression level of the transgene was constant up to 8 h after NAA addition, and clearly upregulated after 23 h while the Aux/IAA gene was induced already within 2 h ( Figure 5C), which suggests that the effect of auxin on PpAct5 promoter-driven expression is indirect.
Although PpAct5 was slowly induced upon auxin treatment, our results can be explained by a link between auxin signaling and actin. The debundling of actin is necessary for an efficient cellular response to auxin and promotion of auxin transport (Nick et al., 2009;Guillory and Bonhomme, 2021). Both, the natural and the artificial auxins, indole-3-acetic acid (IAA) and NAA, respectively, induced debundling of actin filaments. In Arabidopsis, overexpression of the actin-binding domain of plant fimbrin (GFP-FABD2) decreased auxin transport significantly due to reduced actin dynamics (Zaban et al., 2013). Furthermore, the overexpression of the same domain led to desensitization of auxin response in tobacco BY2 cell suspension cultures. These observations reveal that the role of actin in auxin signaling is not just a structural effect, since actin is involved in polar auxin transport (Bennett et al., 2014;Huang et al., 2017).
We observed a high increase in MFHR1 production in the bioreactor with NAA after the semi-continuous operation was ceased and the bioreactor was run in a batch phase for 2 days before ending the whole process ( Figure 5A). This observation suggests that biomass concentration plays an important role in MFHR1 production, which concurs with our phenomenological model. Therefore, batch, fed-batch and semi-continuous bioreactors with lower dilution rates and auxin supplementation were evaluated. NAA (10 µM) was added at day 3 and 4 in the batch bioreactor process, reaching 823 µg MFHR1/g FW (150 mg MFHR1 accumulated over 7 days) and 200 µg MFHR1/g FW with line A and B, respectively ( Figure 6). Compared to batch bioreactor without exogenous NAA, this treatment led to an approx. 230% increment in MFHR1 productivity with line B ( Figure 6A). Product formation decreased at day 8, which could be attributed to nitrate limitation. With these conditions and auxin supplementation, we reached 8-fold the productivity achieved in an earlier study under semi-continuous operation mode for line A ( Figure 6B) (Top et al., 2019) and 2.7-fold the specific productivity achieved under the same conditions but with auxin supplementation since day 6 ( Figure 5D). During the bioreactor run with exogenous NAA (day 3), caulonemal cells started developing from day 2, even before auxin addition, and became more evident from day 4. Big pellets were evident from day 8 (Supplementary Figure S8). Auxins are mainly extracellular hormones (Reutter et al., 1998); endogenous auxin levels and cell density can influence the development of caulonema already at day 2. In Funaria Pp3c8_14720V1.1), PpAct7 (Pp3c3_33410V3.1) and MFHR1 determined by qRT-PCR of 10 µM NAA-treated bioreactor samples. Data were normalized using two reference genes (coding for EF1-α and 60S ribosomal protein L21) and a point immediately before the addition of NAA (indicated by the arrow) was used as calibrator. A comparison between the protein production measured by ELISA and the expression of mRNA by qRT-PCR is included (data obtained from the bioreactor shown in (A), NAA added at day 7). Mean ± standard deviations (SD) from three technical replicates are shown.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2022 | Volume 10 | Article 837965 13 hygrometrica, the effect of auxin is highly dependent on the cell density, and the manipulation of both can lead to cultures enriched predominantly with caulonema. At low cell densities (<200 mg/L) basal caulonema development depends on endogenous auxin levels, and below 50 mg/L approx. 10% of cells starts differentiation (Johri, 2020).
As the steep increase in MFHR1 accumulation was not immediate ( Figures 5A,B), different time points for the beginning of NAA treatment in fed batch and semicontinuous processes were evaluated using line B (Figure 7). NAA (10 µM) was added at day 3, 5 or 7, respectively. Fed-batch operation started at day 7 by adding 200 ml/d (5× Knop ME), with a measured nitrate concentration of 2,500 mg NO 3 − /L. Semicontinuous operation started at day 7, and the dilution rate was reduced from 0.4 to 0.2 d −1 (1 L culture was replaced daily instead of 2 L) to allow the biomass to increase (Figure 7). Under these conditions in semi-continuous mode, the highest specific protein productivity was between 134 and 178 µg MFHR1/g FW. However, washing out of the cells was not prevented by decreasing the dilution rate in the case of NAA supplementation. As semi-continuous operation was not the best mode to achieve a high MFHR1-specific productivity we abstained from evaluation of NAA supplementation at day 5.
In fed batch the time point of NAA addition was crucial. The highest specific recombinant protein productivity was obtained in fed-batch mode. When NAA was added at day 3 it reached up to 280 µg MFHR1/g FW, while only a third of this concentration was achieved when NAA was added later, at day 5 or 7 (Table 3; Figure 7A). The total amount of MFHR1 was doubled by the use of NAA from day 3 in fed-batch and the specific protein productivity increased by 260% and 230% compared to the fed-batch and batch operation without NAA, respectively (Tables 2, 3).
Parameters of semi-continuous operation supplemented with auxin also influenced the product yield, which increased by 250% by decreasing the dilution rate, and delaying the starting day of the semi-continuous operation to allow higher biomass accumulation ( Figures 5A, 7B). These results suggest that cell density or tissue age play a role in protein production. The daily medium exchange can dilute metabolites involved in signaling, such as cytokinin and auxin that accumulated in the medium or influence the differentiation of the tissue, with a negative impact on recombinant protein production. However, this operation mode also avoids the steep decrease in MFHR1 production observed under batch or fed-batch conditions.
A drop in MFHR1 concentration was observed in fed-batch in each evaluated condition approx. after 8 days, similar to the behavior without auxin ( Figure 7A), and is not related to nitrate limitation. A drop in protein accumulation might be due to a decrease in the expression of the transgene but also attributed to inefficient protein extraction, secretion of the protein to the extracellular medium or degradation. In order to further understand this decrease in protein concentration, we analyzed the transgene expression levels from fed-batch bioreactor samples treated with NAA at day 7. MFHR1 productivity quantified by ELISA was congruent with the transgene expression the first 3 days after NAA exposure ( Figure 7C). The change in expression levels remained relatively stable from day 1-5 after NAA exposure, and decreased 2 days later, while MFHR1 productivity decreased after 4 days ( Figure 7C). The expression patterns of the MFHR1 transgene, Aux/IAA gene and PpAct7 were similar to the results observed under semi-continuous operation ( Figures 5B, 7C). Our results explain only partially the kinetics of MFHR1 concentration, as a lower rate of production caused by a lower expression of the transgene cannot be solely responsible for a decrease in concentration. The mechanisms involved in a drop of MFHR1 remain unclear. As the culture supernatant turns slightly brownish at the time where protein concentration decreases, we speculate that accumulation of polyphenolic compounds is a crucial factor. Moreover, this phenomenon occurred faster upon auxin treatment than in the control and could be a stress signal leading to accumulation and oxidation of these compounds and impairment of growth rate. During the extraction of the protein, covalent bonding between protein and phenolic compounds can occur and alter protein functionality and physico-chemical properties (Menkhaus et al., 2004), which might affect the estimation of MFHR1 concentration by ELISA. Likewise, auxin transport and response is inhibited when there is an accumulation of flavonoids (McCurdy et al., 2001;Moody et al., 2021), which might explain the difficulty to increase protein production in the fed-batch mode with exogenous auxin even more.
In summary, we propose a phenomenological model to predict recombinant protein production in moss bioreactors operated in batch and semi-continuous operation, which can be used to optimize recombinant protein yield and biomass. This kinetic model can be applied to other plant cell suspension cultures to predict biopharmaceutical production, just by changing the terms of the energy source. Furthermore, exogenous auxin increased specific recombinant protein production by 470% in shaken flasks, and by up to 230 and 260% in moss bioreactors operated in batch and fed-batch mode, respectively. Under our conditions, the semicontinuous operation mode was inferior to batch or fed-batch for the production of MFHR1. A change in bioreactor operation mode along with NAA supplementation led to an increase of MFHR1specific productivity up to 8-fold (0.82 mg MFHR1/g FW) compared to the values reported for line A (Top et al., 2019). We conclude that PpAct 7 and PpAct 5 are not rapidly induced upon auxin treatment, which suggests an indirect response to auxin leading to upregulation of a transgene driven by the actin 5 promoter in Physcomitrella. The application day of auxin and the bioreactor operation mode are important factors to increase the recombinant protein specific productivity. In this case, MFHR1 production follows a mixed-growth associated pattern. Therefore, a high growth rate and a high biomass density are essential.
We suggest that the auxin effect on biopharmaceutical production in moss bioreactors can be extrapolated to other plant-cell suspension cultures. Auxin stimulates the mitotic activity in tobacco-cell suspension cultures (Huang et al., 2017), which is advantageous to producing therapeutic proteins, when their expression is driven by a constitutive promoter. Although the medium of plant-cell suspension cultures usually includes cytokinins and auxins, the optimization of auxin concentrations is worth being considered, e.g., medium optimization of tobacco-cell suspension cultures producing the antibody M12 revealed a strong influence of the auxins IBA and 2,4-D (Vasilev et al., 2013).
The effect of auxin was not included in the phenomenological model due to lack of understanding of the relationship between auxin and the macroscopic variables biomass, nitrate and protein concentration. However, the data generated in this study will be the basis for implementing a data-driven model involving the effect of auxin in fed-batch bioreactors based on machine-learning approaches to further optimize recombinant protein production.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
NR-M performed experiments and built the dynamic model. SS performed experiments at shaken-flask scale. NR-M, JP, CP, RR, and ED conceived the study and wrote the manuscript. All authors read and approved the final manuscript.

FUNDING
This work was supported by the Deutscher Akademischer Austauschdienst DAAD (to NR-M) and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC-2189 (CIBSS to RR). We acknowledge support by the Open Access Publication Fund of the University of Freiburg.