Comparative Transcriptomics Analyses of a Vero Cell Line in Suspension Versus Adherent Culture Conditions


 The Vero cell line is the most used continuous cell line for viral vaccine manufacturing. Its anchorage-dependent use renders scaling-up challenging and operations very labor intensive which affects cost effectiveness. Thus, efforts to adapt Vero cells to suspension cultures have been invested but hurdles such as the long doubling time and low cell viability remain to be addressed. In this study, building on the recently published Vero cell line annotated genome, a functional genomics analysis of the Vero cells adapted to suspension is performed to better understand the genetic and phenotypic switches at play during the adaptation of Vero cells from anchorage-dependent to suspension cultures. Results show a downregulation of the epithelial to mesenchymal transition (EMT) pathway, highlighting the dissociation between the adaptation to suspension process and EMT. Surprisingly, an upregulation of cell adhesion components is observed, notably the CDH18 gene, the cytoskeleton pathway, and the extracellular pathway. Moreover, a downregulation of the glycolytic pathway are balanced by an upregulation of the asparagine metabolism pathway, promoting cell adaptation to nutrient deprivation. A downregulation of the adherens junctions and the folate pathways alongside with the FYN gene are possible explanations behind the currently observed low cell viability and long doubling time.


Introduction
Derived from a female Chlorocebus sabaeus (African Green Monkey) kidney, the Vero cell line is susceptible to infection by a wide range of viruses. Consequently, the Vero cell line served as a platform for the development and production of approved vaccines against dengue fever, in uenza, Japanese encephalitis, polio, rabies, rotavirus, smallpox and more recently, Ebola (using a recombinant Vesicular Stomatitis Virus) and Covid-19 (Sinopharm, Sinovac, CoronaVac) [1][2][3], thus representing the most widely used continuous cell line for the production of viral vaccines over more than 40 years of manufacturing experience [4].
Most cell-based vaccine production platforms, including Vero cells, are adherent cells cultured in an anchorage-dependent manner thus hindering scale-up attempts due to the time and cost challenges posed by the development of large-scale anchorage-dependent cell culture platforms [5][6][7]. Therefore, in an effort to circumvent those challenges, stable cell lines adapted to suspension culture have been proposed over the years [8]. In the case of Vero cells in particular, signi cant efforts have been dedicated to developing micro-carriers cultures [9], agitation [10], and medium engineering [11] to achieve scalable processes and industrialization. Despite several reports stating that Vero cells adapted to suspension showed a higher production rate for viruses such as measles virus, rabies virus and vesicular stomatitis virus (VSV) [10,11], some issues such as the low cell viability and long doubling time of these cells were reported underlining the needs for more effective adaptation of Vero cells to suspension by exploring new avenues such as genetic engineering [10,12].
Recent advances in functional genomics and gene editing paved the way for a better understanding and high-throughput engineering of vaccine production cell lines, thus providing new possibilities for cell line development and vaccine bioprocessing intensi cation. With the recent publication of an annotated assembly of the Vero cell genome [13], we propose a transcriptomics analysis of the differences between adherent and suspension Vero cells developed by Shen et al. [11]. This process will help highlight key differentially expressed genes and their impact on the cell's phenotype and their metabolic pathways for a better understanding of the process of adaptation to suspension, thus providing insight for new strategies to successfully adapt Vero cells to suspension cultures and enable streamlined scale-up and industrialization.

Differential expression analysis and Pathway enrichment analysis
Following the DESeq2 differential expression (DE) analysis and an applied p-value cut-off of 10 -4 , among the 6627 DE genes that have been identi ed as highly signi cant (p-value < 10 -4 ), 1753 genes were identi ed as highly differentially expressed, with a |log2 fold change| > 2 ( Figure 1, Supplementary table S1).
The upregulated DE genes are highly enriched in key pathways such as vacuole, regulated exocytosis, plasma membrane bounded cell projection morphogenesis, actin lament-based process and regulation of cell adhesion with a p-value < 10 -10 ( Figure 2). On the other hand, downregulated DE genes are highly enriched in adherens junction, MYC targets, RNA processing and cell cycle related pathways with a pvalue < 10 -10 ( Figure 3).

Metabolic pathway analysis
The extraction of metabolic genes and their scoring based on their relation with established metabolic pathways revealed an upregulation of the metabolic pathway of gluconeogenesis, galactose, pyrimidine, glycine, threonine, but also the serine pathway which plays a key role in unrestrained cell cycle progression [14], the asparagine metabolism pathway that is linked to cells adaptation to nutrient deprivation and/or hypoxia [15] and the mitochondrial one-carbon metabolism pathway which is implicated in rapid cell proliferation [16]. On the other hand, the adaptation to suspension led to a downregulation of key metabolic pathways such as proline, folate, aspartate, lipids and glycolytic pathways ( Figure 4, Supplementary table S2). Notably, de ciencies in the folate metabolism pathway were linked to growth limitations of BHK-21 in suspension culture [17].

Gene set enrichment and network topology analysis
Gene set enrichment analysis (GSEA) with the hallmark gene set showed a signi cant upregulation of the interferon gamma response pathway and a down regulation of MYC targets (variant 1 and 2), E2F targets pathways but most importantly the epithelial mesenchymal transition pathway ( Figure 5). Only pathways with an FDR(False discovery rate)< 0.05 were considered. In order to identify key genes involved in protein-protein interaction (PPI) networks, network topology analysis(NTA) was done for both upregulated genes (Table 1) and down-regulated genes ( Table 2). The signaling, response to chemical and anatomical structure pathway associated with CRYAB, RHOU, ESR1, and C3 are upregulated alongside the cell adhesion pathway which is associated with CDH18. On the other hand, pathways associated with FYN and SMAD9 (cell differentiation, system development, cellular response to growth factor) and the positive regulation of cell migration pathway associated with TRIP6 are downregulated.

Discussion
Over the years signi cant efforts have been made to adapt Vero cells to suspension in order to engineer high-throughput and scalable vaccine production platforms, however, those efforts were limited by hurdles such as low cell viability and long doubling time. In order to design more e cient strategies for successful adaptation to suspension cultures maintaining acceptable cell viability and doubling time, it is necessary to better understand the genetic and phenotypic changes triggered by the adaptation to the suspension state. Thus, we propose in this paper a comparative functional genomics analysis of adherent and suspension Vero cells from the gene to the protein interaction network level.
Indeed, correlations between differential expression analysis ( Figure 2, 3), metabolic pathway enrichment analyses (Figure 4), gene set enrichment analyses ( Figure 5) and network topology analyses (Tables 1, 2) highlight key events such as the upregulation of immune response, exocytosis and vacuole pathways which are involved in the storage of waste and other exocytosis molecules.
Furthermore, several events of checks-and-balances were observed across the different analysis methods used (Table 3): lipids metabolism and lipid pathways are also upregulated via pathway enrichment while the metabolic analysis revealed a downregulation of CPT1 and palmitate associated metabolic reactions which regulate fatty acid oxidation in mitochondria and whose upregulation impairs glucose homeostasis [18]. In parallel, the gluconeogenesis metabolic pathway is upregulated while the glycolytic pathway is downregulated thus hindering ATP generation. That resulting stress is met with not only a downregulation of the proline metabolic pathway which constitute a checkpoint that is reported to promote proline accumulation during stress [19], but also an upregulation of asparagine metabolism which promotes the cell adaptation to nutrient depravation and/or hypoxia [15] thus promoting cell survival.
Surprisingly, as observed in HEK293 cells adapted to suspension [20], cellular component organization pathways associated with cell adhesion such as actin lament-based process, regulation of cell adhesion, apical part of cell, plasma membrane bounded cell projection morphogenesis and extracellular matrix are upregulated via pathway enrichment analysis, while NTA showed an upregulation of the cell adhesion PPI network with CDH18 as its central gene and which can be considered as possible target for engineering in order to improve the cells adaptation to suspension. This upregulation of cell adhesion related genes could be due to the cells' attempt to restore the attachment to culture surfaces and surrounding cells which could explain the aggregates that are often observed in Vero cells suspension cultures and the cell rings that form on the suspension culture dishes.
On the other hand, the adherens junction pathway, which regulates cell-cell adhesion and is essential for viability (via the control of cell proliferation, polarity, shape, motility and survival) [21], is downregulated alongside the aspartate metabolic pathway, the mitochondrial 1-carbon metabolic pathway, the MYC and E2F targets pathway, which could explain the low cell density observed during Vero cells adaptation to suspension.
Moreover, the pathways related to cell division, mitotic cell cycle phase transition are downregulated via pathway enrichment, which is con rmed by the downregulation of the folate metabolic pathway which leads to cell cycle arrest at G0/G1 [22] and the downregulation of FYN which is central to the networks associated with the control of cell growth, thus, providing some insight in the origin of the long doubling time observed in Vero cells adapted to suspension culture. Nonetheless, Vero cells attempt to balance that effect on doubling time via the upregulation of the glycine, threonine and serine pathway which associated with an unrestrained cell cycle progression [14], and the upregulation of RHOU, ESR1 which are central to the anatomical structure morphogenesis pathway and more precisely cell proliferation and antiapoptotic regulations.
Lastly the epithelial to mesenchymal transition (EMT) pathway is downregulated in Vero cells adapted to suspension, thus highlighting the fact that the adaptation to suspension is not associated with EMT as previously shown with HEK293 cells adapted to suspension [20].
To conclude, we present in this paper key genes pathways at play during the adaptation of Vero cells to suspension and their complex checks-and-balances which could assist in the successful adaptation of Vero cells to suspension. Indeed, those key genes, notably associated with cell adhesion or hindering cell viability or doubling time could be potentially targeted via gene editing as new strategies for the adaptation to suspension, but also the observed competitions between the regulation of competing pathways can be studied more in detail via targeted perturbations using gene editing tools such as CRISPR [23].

Cell lines and culture media
The adherent Vero WHO cell line studied in this work was at passage 138 (Neovacs). This cell line was itself derived from a vial of Vero ATCC CCL-81 which was send to WHO at passage 124 for analysis and establishment of the Vero WHO master cell bank approved for vaccine production. The cells were grown in static culture at 37°C and 5% CO 2 in a humidi ed incubator (Infors HT, Switzerland). Cells were passaged twice weekly using TrypLE Express (Thermo Fisher Scienti c) as dissociation reagent. A serumfree adapted sub-cell line grown in OptiPRO medium (Thermo Fisher Scienti c) supplemented with 4 mM GlutaMAX (Thermo Fisher Scienti c) was cryopreserved at a passage number of 151 in OptiPRO medium supplemented with 4 mM GlutaMAX and 10% DMSO (Sigma, USA).
The suspension Vero cell line was provided by the National Research Council (NRC) of Canada [11]. The cells were maintained at 37°C, 135 rpm and 5% CO 2 in a humidi ed Multitronorbital shaker (Infors HT, Switzerland) and were cultivated in 20 mL working volume of either IHM03 medium, provided by the NRC, or in MDXK medium (Xell AG, Germany), supplemented with 4 mM GlutaMAX (Thermo Fisher Scienti c, USA) [4].
For genome analysis, Vero WHO cells at passage 153 were washed in PBS (Wisent, Canada), harvested using TrypLE Express and centrifuged at 300 × g for 5 minutes. Cell pellets of around 6 million cells were quickly frozen in a mixture of dry ice/ethanol and stored at −80°C until further analysis. Adherent and suspension samples were generated in triplicates.

Differential gene expression analysis
Total RNA sequencing (TrueSeq) was performed using Illumina NovaSeq6000 Sprime v1.5, PE100. Following standard quality control, the reads were rst aligned to the recently published Vero cell genome [13] using STAR [24] and the resulting BAM les were sorted by name using SAMtools [25] before read count. Transcripts were quanti ed using featureCounts [26]. Differential expression analysis of the raw read counts was done using DESeq2 [27] and quality control graphs were produced using DESeq2 and R package. The resulting differentially expressed gene list was ltered with a p-value cut-off of 0.0001.

Downstream analysis of differentially expressed genes
The differentially expressed (DE) genes were ranked based on their log2 fold change. Upregulated and downregulated DE genes between adherent and suspension cell lines were used for pathway enrichment analysis via the Metascape webtool and plotted using default parameters [28].
In order to identify the metabolic deregulations that distinguish adherent cells from suspension cells, metabolic genes from our list of differentially expressed genes were extracted and their Reaction Activity Score (RAS) were computed by solving Gene-Protein-Reaction (GPR) association rules based on the HMRcore metabolic network model via MaREA (Metabolic Reaction Enrichment Analysis) [29]. A statistical comparison of the RAS computed for adherent samples and suspension samples was done via the Kolmogorov-Smirnov test and a metabolic map was generated.
WebGestalt (WEB-based GEne SeT AnaLysis Toolkit) [30] was used for gene set enrichment analysis with the hallmark gene set collection. To nd differentially expressed pathways of genes between adherent and suspension cell lines, gene sets were ltered and the top 20 gene sets with an adjusted p-value lower than 0.05 were considered as signi cantly changed.
The upregulated part of the gene list generated by DESeq2 was ltered to consider genes with a |log2 fold change| > 2 for Network Topology Analysis (NTA) based on the Network Retrieval & Prioritization construction method [30] by rst using random walk analysis to calculate random walk probability for the input gene IDs (seeds), then identifying the relationships among the seeds in the selected network to return a retrieval sub-network where the top 20 genes with the top random walk probability are highlighted. Indeed, assuming a tight connection between mechanistically important genes and a random distribution of other genes on the network, the Network Topology-based Analysis (NTA) uses random walk-based network propagation by identifying those genes which are potentially biologically signi cant. Our input gene IDs (upregulated genes previously ltered) were used as seeds and, based on its overall proximity (quanti ed by the random walk similarity) to the input seeds, each gene in the protein-protein interaction (PPI) network was attributed a score. Then the statistical signi cance of those scores was calculated via two p-values: a global p-value which signi cance is the result of a non-random association between the gene in the PPI network and the input seeds; and a local p-value which signi cance ensures that the gene did not acquire a signi cant association with the input seeds simply because of network topology.
Finally, enrichment analysis of the retrieved sub-networks was done using the PPI BIOGRID   DESeq2 differential expression data quality control: Volcano plot of DE genes with applied p-value and log2 fold change cut-offs. FDR: False Discovery Rate; LogFC: log2 fold change.

Figure 2
Network of terms enriched by upregulated DE genes: (a) colored by cluster ID, where nodes that share the same cluster ID are typically close to each other; (b) colored by p-value, where terms containing more genes tend to have a more signi cant p-value. Each node represents one enriched term and edges link similar terms (the edge thickness is proportional to the similarity between terms).

Figure 3
Network of terms enriched by downregulated DE genes: (a) colored by cluster ID, where nodes that share the same cluster ID are typically close to each other; (b) colored by p-value, where terms containing more genes tend to have a more signi cant p-value. Each node represents one enriched term and edges link similar terms (the edge thickness is proportional to the similarity between terms).

Figure 4
Adherent and suspension comparative metabolic map. Blue and red arrows refer respectively to downregulated and upregulated reactions. Dashed gray arrows refer to non-signi cant dysregulations according to Kolmogorov-Smirnov test with p-value 0.01. Solid gray arrows refer to reactions with a variation lower than 20%.