{'fraction_non_alphanumeric': 0.028354233714349128, 'fraction_numerical': 0.025024969481744534, 'mean_word_length': 5.124893797790994, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Defining lncRNAs Correlated with CHO Cell Growth and IgG Productivity by RNA-Seq

Summary How the long non-coding RNA (lncRNA) genome in recombinant protein producing Chinese hamster ovary (CHO) cell lines relates to phenotype is not well described. We therefore defined the CHO cell lncRNA transcriptome from cells grown in controlled miniature bioreactors under fed-batch conditions using RNA-Seq to identify lncRNAs and how the expression of these changes throughout growth and between IgG producers. We identify lncRNAs including Adapt15, linked to ER stress, GAS5, linked to mTOR signaling/growth arrest, and PVT1, linked to Myc expression, which are differentially regulated during fed-batch culture and whose expression correlates to productivity and growth. Changes in (non)-coding RNA expression between the seed train and the equivalent day of fed-batch culture are also reported and compared with existing datasets. Collectively, we present a comprehensive lncRNA CHO cell profiling and identify targets for engineering growth and productivity characteristics of CHO cells.


INTRODUCTION
Many recombinant protein biopharmaceuticals are expressed in mammalian expression systems due to the ability of such systems to correctly fold, assemble, and undertake ''human-like'' post-translational modifications and secrete the target protein out of the cell (Walsh, 2010). Of mammalian cell expression systems one predominates, with more than 60% of mammalian made biotherapeutic proteins produced from cultured Chinese hamster ovary (CHO) cells (Kunert and Reinhart, 2016;Leu et al., 2004;Mead et al., 2015;Povey et al., 2014;Walsh, 2010). Fed-batch culture is currently the most common bioprocess used for the industrial production of proteins in CHO cells, generating increased cell concentrations (and hence biomass) and sustained culture viability compared with batch culture, ultimately resulting in higher productivity and final product yields (Durocher and Butler, 2009;Pan et al., 2017;Wong et al., 2006). The introduction of small-scale parallel bioreactors allowing automated sampling and continuous control of fundamental culture parameters, including pH, stirring, and temperature has enhanced the ability to screen a wider range of culture parameters and cell lines, leading to improved upstream development timelines and experimental throughput (Bareither and Pollard, 2011).
The ambrâ15 cell culture system (Sartorius Stedim Biotech) has been shown to give similar cell growth and productivity data to those achieved in larger-scale stirred bioreactors, enabling more accurate predictions compared with shake flasks on the behavior of a cell line at larger scale (Alsayyari et al., 2018;Janakiraman et al., 2015;Nienow et al., 2013;Rouiller et al., 2016). This capacity to conduct small-scale experiments under controlled conditions, of a highly predictive nature at larger scale, allows the investigation of the behavior of different cell lines under alternative feeding regimes to determine how each respond. Indeed, recent reports state that the ambrâ15 small-scale automated and controlled bioreactor system provides an excellent scale-down model to facilitate studies on multiple cell lines under controlled industrially relevant conditions to identify robust targets linked to productivity for cell engineering and material and data for future regulatory submissions (Sandner et al., 2019).
Despite advancements in the ability of CHO cells to reach higher cell concentrations and generate increasing amounts of target biotherapeutic proteins, particularly monoclonal antibodies (mAbs), there remains a desire to further understand the limitations upon CHO cell phenotypes and to engineer cells for the production of more difficult-to-express products (Godfrey et al., 2017;Jossé et al., 2018;Mead et al., 2009). One approach that has been applied toward improving our understanding of the limitations on CHO cell growth and recombinant protein production is the field of transcriptomics (Tamo saitis and Smales, 2018). Transcriptomic studies in particular could benefit from generating material from controlled miniature bioreactors that predict behavior at larger scale, as an issue of such previous studies is reproducibility and robustness across different transcriptomic datasets, given the high heterogeneity of CHO cell lines and their intrinsic genetic instability Wurm, 2013;Wurm and Wurm, 2017). The availability of CHO cell and Chinese hamster genome sequences Xu et al., 2011) has greatly enabled omics-based studies (Faustrup Kildegaard et al., 2013), and since there has been an increasing number of publicly available databases for different CHO cell lines. However, the focus of these studies has been on either coding genes or microRNAs (miRNAs), with few studies investigating other classes of RNAs and their impact on CHO cell behavior (Singh et al., 2018;Tamo saitis and Smales, 2018).
Here we investigate the long non-coding transcriptome in CHO cells during fed-batch culture under controlled bioreactor conditions. Since the unraveling of multiple organisms genomes, particularly eukaryotic genomes, associated with the development of high-throughput sequencing technologies, new classes of noncoding RNA have been identified (Djebali et al., 2012). Among these, a class of transcripts known as long noncoding RNAs (lncRNAs) was identified. lncRNAs are defined as transcripts longer than 200 nucleotides that lack a significant open reading frame (ORF), are usually transcribed by RNA polymerase II, and are spliced with, or without, 3 0 polyadenylation (Kashi et al., 2016;Kung et al., 2013;Wilusz, 2016). In the nucleus, cis-acting lncRNAs regulate the chromatin state and transcription of nearby genes, whereas trans-acting lncRNAs can recruit RNA-binding proteins to form chromatin-modifying complexes and modulate splicing or organize functional nuclear domains (Kopp and Mendell, 2018). When transported to the cytoplasm, lncRNAs act at a post-transcriptional level by promoting specific mRNA translation or turnover and by competitively binding microRNAs (miRNAs), attenuating the repression of target genes (Geisler and Coller, 2013). The wide range of processes involving lncRNAs suggests that some of these may be potential cell engineering targets to rewire CHO cell phenotypes for enhanced cell growth and/or recombinant protein production and quality without placing a translational burden on the cell compared with overexpression of coding genes. However, the majority of our knowledge around lncRNAs comes from studies in model organisms related to disease and development (Perry and Ulitsky, 2016;Schmitt and Chang, 2016) with lncRNAs poorly annotated in the CHO cell genome and little known about their role in defining CHO cell phenotypes.
A recent report described the lncRNA landscape in CHO cells, showing regulated expression of thousands of lncRNAs under batch and fed-batch conditions over time (Vito and Smales, 2018). Others have demonstrated the potential power of lncRNA cell engineering to manipulate the cells' ability to produce target recombinant proteins with the engineering of a class of lncRNAs named SINEUPs enhancing the translation of target-specific mRNAs in various mammalian cell factories (Patrucco et al., 2015;Zucchelli et al., 2016). However, the limited number of studies and poor annotation of non-coding regions in the Chinese hamster genome means that transcriptomics across multiple cell lines associated with phenotypes of interest under industrially relevant and controlled conditions is required to identify lncRNAs whose manipulation may enhance mammalian cell factories' ability to generate secreted target products (Vishwanathan et al., 2014). Here we present a comprehensive coding and non-coding, particularly lncRNA, transcriptome analysis using RNA-Seq of five IgG1 producing CHO cell lines and one stable pool harboring the plasmid cassette without genes encoding for the IgG1, cultivated under fed-batch conditions in an ambrâ15 system to unveil lncRNA targets for cell engineering. The RNA-Seq datasets and analyses are made openly available to the community to promote further studies and comparisons, providing a detailed CHO cell lncRNA transcriptomic resource with a focus on particular genes and pathways.

RESULTS AND DISCUSSION
Following our previous landscape of lncRNA expression in CHO using microarray (Vito and Smales, 2018), we set out to provide an analysis of the lncRNA transcriptome using RNA-Seq in CHO cells producing three different model IgG1 monoclonal antibodies during fed-batch culture, defining those lncRNAs expressed in CHO cells and the flux of these during culture and between cell lines. To do this, we undertook RNA-Seq analysis on a panel of IgG-expressing CHO DHFR-cell lines (these cell lines being generated from a modified dihydrofolate reductase-deficient (DHFR-) CHO DG44 host cell line, (Urlaub et al., 1983)), sampling throughout fed-batch cultures in an ambrâ15 microbioreactor system generating profiles of the flux of coding RNAs and lncRNAs.

Analysis of Fed-Batch Culture Samples
The DAVI dataset included the 3068, 3080 and 3077 IgG1-producing cell lines and the null pool 3478, cultured for 12 days (Table 1). These cell lines maintained culture viabilities >80% throughout the 12 days of fed-batch culture, reaching a peak in viable cell concentration (VCD) at 14.48 x 10 6 viable cells/mL for cell line 3077, whereas 3080 had the lowest VCD throughout culture ( Figure 1A). 3080 reached the highest titer at 1.80 g/L and productivity (Qp) at 31.71 pg/cell/day (3080) associated with the highest amount of lactate accumulated during the first few days of culture, with a peak concentration of 2.2 g/L ( Figures 1B-1D). The glutamate, glutamine, and ammonia profiles were similar across all cell lines (Figure S1). Overall, the highest productivity was associated with the highest lactate accumulation causing the lowest VCDs for 3080, whereas the rest of the producers did not show this association compared with the null pool.
The JCE dataset included the 3068, 3080, and 4384 IgG1-producing cell lines and the 3936 IgG1-producing pool, cultured for 14 days (Table 1). Once again, 3080 was the cell line showing the lowest VCDs and the highest accumulation of lactate at approximately 3.5 g/L at Day 5, whereas culture viability profiles were broadly comparable across cell lines, being within the range of 60%-80% on the last day of culture (Figures 1E and 1G). When the Qp was calculated, the 3068 and 3080 cell lines had the highest Qps at 35.20 and 35.81 pg/cell/day respectively, whereas 3936 had the lowest Qp, as would be expected from a pool of clones ( Figure 1H). The glutamate, glutamine, and ammonia profiles were again broadly similar across cell lines ( Figure S1). In summary, 3080 repeated the trend observed in the DAVI dataset of lowest VCD associated with highest productivity and highest accumulation of lactate, but 3068 followed closely in productivity.
RNA Sequencing of ambrâ15-Generated Samples and Subsequent Analysis of the Data: The DAVI Experiment As described in the methods section, samples for RNA-Seq were collected in duplicate at Day 4 and Day 12 of fed-batch culture from the ambrâ15 reactors. Gene counts were calculated using featureCounts, then the differential expression analysis was undertaken with the R/Bioconductor package DESeq2 (Love et al., 2014). Clustering analysis of the RNA-Seq data revealed that the samples showed a consistent hierarchical clustering for the biological replicates and an evident separation between the two time points based on gene expression (Figures 2A and 2B). This separation was confirmed by principal component analysis (PCA); however, an additional layer of clustering emerged among the producers upon PCA with a difference in Qp, clearly showing clusters formed of 3077 and 3068 (23.05 and 26.53 pg/cell/day) that were distanced from cell line 3080 (31.71 pg/cell/day) and the null-pool 3478 ( Figure 2C). Overall, this preliminary cluster analysis based on total gene expression suggests a grouping of cell lines directly related to differences in Qp more than maximum viable cell concentrations for this dataset.
Differential transcript expression (DE) analysis was then conducted using the DESeq2 R software package, setting a fold-change (FC) threshold of 1.50. DE genes were considered significant if the adjusted p value for this FC threshold calculated using the Benjamini-Hochberg method was below 0.10. As expected, the longest distances observed in clustering and PCA analyses corresponded to a higher number of DE genes overall ( Figure 3). Consequently, samples collected at Day 12 showed the highest variability when compared between them or against Day 4 samples, with the exception of 3077 vs 3068, which indeed clustered together in PCA ( Figure 2). Within the identified DE genes, lncRNAs made up 10-30% of the total number identified ( Figure 3). A representative group of differentially expressed coding and non-coding RNAs were then selected (see section on identification of differentially expressed lncRNAs as potential cell engineering targets) for RT-qPCR validation, resulting in a positive correlation between the foldchanges measured by RNA-Seq and by RT-qPCR (Data S2).

RNA Sequencing of ambrâ15-Generated Samples and Subsequent Analysis of the Data: The JCE Experiment
In the JCE experiment, samples for RNA-Seq were collected in triplicate from the seed train (ST) flasks and at Day 4 and Day 7 of fed-batch culture. The samples showed a hierarchical clustering for each biological triplicate, but the separation into groups as observed in the DAVI experiment was not as evident ( Figures 2E and 2F). PCA revealed a similar pattern where only a clear separated cluster composed of the 4384 cell line samples was distinguishable from the other samples ( Figure 2G). The 4384 clones showed the highest VCD throughout culture ( Figure 1E) and the observed distance in the PCA clustering reflected in the high number of DE genes identified when comparing this cell line against the others, especially on Day 4 ( Figure 3). Although less predominant, the presence of sub-clusters grouped by time point with the same general trend can be observed among the rest of the samples ( Figure 2G). We then compared the gene expression profiles of the seed-train cultures that were used to start the fed-batch process, obtained from cells during logarithmic growth phase, with the Day 4 gene expression profiles of the fed-batch culture experiments. Although cells from the seed train and fed-batch Day 4 ambrâ15 bioreactor might be expected to be in a similar growth and metabolic state, we found significant differences in gene expression numbers at this early stage of culture, particularly for 3068 where the number of DE genes identified was 761 ( Figure 3). Overall, the hierarchical clustering, PCA, and DE analysis suggest the 4384 cell line has a distinct transcriptional landscape, whereas 3068, 3936, and 3080 have much closer gene expression profiles. Further, the seed train samples of each cell line show, to varying degrees, different gene expression profiles than that of cells taken from the fed-batch cultures in an equivalent growth phase.

Investigating Pathway Enrichment in DE Genes
KEGG pathway functional enrichment of the RNA-Seq datasets based on statistically significant differentially expressed genes showed two distinct patterns across the datasets. Firstly, a major theme of enrichment in the DAVI dataset was in the Replication and Repair area, where DE genes were found to be enriched in DNA replication, base excision repair, nucleotide excision repair, mismatch repair, homologous recombination, and Fanconi anemia pathways among the 3077, 3068, and 3478 cell lines when comparing Day 12 versus Day 4 expression profiles within the same cell line ( Figure 4). Interestingly, the only cell line in the DAVI dataset in which none of these pathways was enriched was the 3080 cell line. In the DAVI dataset this was the cell line with highest Qp associated with lowest VCD and a clear distinction from the others when the RNA-Seq data were analyzed by PCA. On the other hand, comparing the different cell lines to each other among Day 4 or Day 12 did not reveal any enriched pathways related to genome maintenance ( Figure S2). Thus, within a given cell line changes in genome maintenance pathways were observed over time between days 4 and 12 of culture; however, when same time points of different cell lines were compared this was not observed.
We then applied the same enrichment analysis to the JCE dataset, and surprisingly none of the pathways involved in genome maintenance were enriched within the DE genes between Day 4 and Day 7 ( Figure S3). An analysis of the seed train versus fed-batch DE genes for the 3080 and 3068 cell lines revealed enrichment of the DNA replication pathway and for the seed train versus 4384 cell line enrichment of the nucleotide excision repair, mismatch repair, and DNA replication pathways ( Figure S4). The pathway enrichment analysis also consistently revealed the enrichment of the PI3K-Akt signaling pathway, focal adhesion, and ECM-receptor interaction pathways in the DE genes across both the DAVI and JCE datasets. Overall, the enrichment suggests a prominent regulation of genome maintenance mechanisms is conserved across different cell lines at the passage from seed train to fed-batch culture and toward the end of culture, whereas the most differentially regulated pathways at the same stage in culture for the same cell lines are the PI3K-Akt signaling, focal adhesion, and ECM-receptor pathways.
Focusing on individual genes involved in the replication and repair domain, we identified exonuclease 1 (Exo1), Rad51, essential meiotic structure-specific endonuclease 1 (Eme1) and FA complementation group B (Fancb) in common among the top 30 most differentially expressed genes in both the DAVI and JCE datasets ( Figure S5). These genes are involved in a wide range of genome repair mechanisms from mismatch repair to homologous recombination and DNA double-strand break repair, suggesting a co-regulation of multiple facets of genome maintenance and the importance of high fidelity in these pathways to maintain cell integrity, viability, and growth as culture progresses.

Mapping of Long Non-coding RNA Expression During Fed-batch Culture
As outlined in the introduction section, a key aspect of this study was to generate a detailed description of lncRNA expression under controlled bioreactor conditions in different IgG-producing CHO cells and to identify lncRNAs whose manipulation may enhance the CHO cell factory's ability to generate secreted target products. In order to identify non-coding RNAs, all the significant differentially expressed non-coding genes were filtered based on NCBI annotation and are shown as a percentage of the total number of differentially expressed genes for each comparison in Figure 3 and in Table S1. The DAVI dataset showed a higher percentage of ncRNAs on average (22.2%) compared with JCE (14.0%) as a percentage of the total DE RNAs identified, most likely due to the higher coverage in the sequencing data for the DAVI dataset. The complete RNA-Seq dataset is provided reporting those lncRNAs identified as being expressed in CHO cells and hence providing a reference for the community to investigate individual lncRNAs in CHO cells (see GEO accession number GEO: GSE140671 and Data S1). The number of differentially expressed genes with an FC R 1.5, adj p value < 0.1 (DE genes were considered significant if the adjusted p value for this fold chnage (FC) threshold calculated using the Benjamini-Hochberg method was below 0.10). In blue is shown the DAVI dataset and in red is shown the JCE dataset while the size of the dot indicates the % of lncRNAs for each comparison. For each comparison, the exact number of genes is indicated on the right.

Identification of Differentially Expressed lncRNAs as Potential Cell Engineering Targets
We then looked to filter and refine the list of DE lncRNAs by counting the occurrence of each transcript in all the DE comparisons, assessing sequence conservation across mammalian species through the discontiguous megablast algorithm and secondary structure prediction based on the RFAM database (Kalvari et al., 2018). Each lncRNA identified via this strategy was then experimentally validated by RT-qPCR to confirm the differential expression between conditions in CHO cells (Data S2). The first lncRNA identified via this approach was Adapt15. LncRNA Adapt15 (also known as growth-arrested DNA-damage inducible gene 7, Gadd7) was discovered in hamster cells, with orthologs identified in the closely related C. longicaudatus and sequence conservation across rodents (Crawford et al., 1996). More recently, Adapt15 has been linked to oxidative lipotoxicity, with knockdown alleviating ER stress and cell death (Brookheart et al., 2009). In the DAVI dataset at Day 12, Adapt15 was consistently upregulated in the 3077, 3068, and 3080 cell lines compared with the null pool 3478 (Table 2). Interestingly, the lowest FC was measured in 3080, the clone with the highest Qp, suggesting Adapt15 expression could be a bottleneck for protein production when the ER is put under stress. In addition, Adapt15 is a known contributor to DNA damage, and 3080 was the only clone not showing DNA-damage-related pathway enrichment in the KEGG analysis previously described (Hollander et al., 1996). Clones with naturally lower levels of Adapt15 could be able to escape these detrimental effects on productivity and genome stability. Conversely, within the JCE dataset, which is focused only on producers and comparison of earlier stages of culture, there was no significant DE of Adapt15 identified, suggesting DE expression of this transcript is observed in producer cell lines compared with a non-producing control later in fed-batch culture.
A second lncRNA that showed significant DE was growth-arrest-specific transcript 5 (GAS5), a multiple small nucleolar RNA (snoRNA) host gene (Smith and Steitz, 1998) with a short ORF and a well-known tumor suppressor lncRNA in human cancer biology (Ma et al., 2016). As a result of cell growth arrest and mTOR pathway activity repression, translation of the GAS5 short ORF is blocked and the transcript accumulates, escaping the nonsense-mediated decay (NMD) pathway that depends on active translation (Ma et al., 2016;Tani et al., 2013). We identified the uncharacterized gene LOC103158913 as a ortholog of mouse GAS5 in CHO using BLAST and RFAM predictions (Table 2). We then investigated its correlation with mTOR expression by integrating the JCE and DAVI measurements for the 3068 and 3080 producers to obtain four time points ( Figure 5). We found a stable ratio between mTOR and GAS5 until Day 12, when mTOR expression became predominant, especially in the highest producer 3080. In addition, GAS5 was upregulated at Day 12 in 3068 compared with the null pool and not in the 3080, whereas mTOR was never   (Kang et al., 2017), the mouse homologue transcript with the corresponding E-value obtained using the disc megablast algorithm, and a summary of the biological function.
differentially regulated. Taken together, literature evidence and the presented data suggest the accumulation of GAS5 may be the result of slower translation efficiency caused by an attenuation of mTOR signaling in the later stages of culture, when recombinant protein production is at its peak, leading to lower productivity.
A third well-characterized lncRNA in other organisms that was identified in our system as differentially regulated was plasmacytoma variant translocation 1 (PVT1). PVT1 is a non-coding oncogene related to poor prognosis in different cancer types (Zhu et al., 2017) with reports of competing endogenous RNA (ceRNA) activity on the miR200 microRNA family and direct interaction with the Myc oncogene (Colombo et al., 2015;Tseng and Bagchi, 2015;Tseng et al., 2014). PVT1 has also been found to induce PI3K/AKT cascade activation, a pathway we found consistently dysregulated in KEGG functional enrichment (Sun et al., 2019). In a previous work (Vito and Smales, 2018), we identified the uncharacterized gene LOC103158906 as an ortholog of mouse PVT1 in CHO using BLAST and secondary structure prediction in RFAM. PVT1 is annotated in C_griseus_v1.0 as LOC103162981, giving the same BLAST alignment scores and RFAM. Myc proto-oncogene protein is predicted to be encoded in CHO as LOC100758352, containing Myc amino-terminal region (Myc_N), helix-loophelix DNA-binding domain (HLH), and a Myc leucine zipper domain (Myc-LZ). PVT1 expression was upregulated by more than two-fold at Day 12 in every producer compared with the null pool in the DAVI dataset and at both Day 4 and Day 7 in the JCE dataset for 3068 compared with 3936, the clone with the highest yield and the pool with the lowest yield, respectively. We then investigated the connection between PVT1 and Myc by measuring the ratio of expression in the JCE and DAVI datasets for 3068 and 3080 ( Figure 5). Although seed-train samples showed an initial discrepancy between the two clones, the ratio was comparable among later time points, a predominant Myc expression at Day 7, when the exponential growth phase is sustained by active proliferation, while PVT1 expression increases substantially at Day 12 ( Figure 5). When assessing the direct interaction probability of CHO PVT1 (XR_478426.2) and Myc (XP_003516054.2) using lncPro and RPISeq prediction tools, the output suggested a strong likelihood of interaction, scoring 85.1 and 0.9, respectively (Lu et al., 2013;Muppirala et al., 2011). Collectively these data suggest a relationship between PVT1, Myc, and proliferation in the early stages of culture when Myc activity is preponderant and a strong overexpression of PVT1 toward Day 12 that may promote IgG production.

Comparison of Data Presented Here with Existing Datasets
A recent meta-analysis publication compared transcriptomics studies in CHO cells, commenting on the difficulty in comparing these datasets but also identifying the most recurrent genes identified as related to Qp and growth (Tamo saitis and Smales, 2018). Here, we find that a number of the coding genes in this study were regulated in agreement with this meta-analysis study (Data S3). The gene ranking first in the list of the meta-analysis was Cd36, a multifunctional glycoprotein acting as receptor for a broad range of ligands of proteinaceous or lipidic nature , which was consistently downregulated in producers compared with the null pool in the DAVI dataset at both time points (Table 3). In addition, Cd36 was strongly downregulated in 3068 and 3080, the cell lines with the highest Qp in the JCE dataset, compared with 3936 and 4384 at both Day 4 and Day 7. Heat shock protein family A member 8 (Hspa8), a molecular   chaperone implicated in the protein quality control system and protection of the proteome from stress (Stricher et al., 2013), was downregulated at Day 12 against Day 4 in the DAVI dataset among every producer cell line and when compared with the null 3478 cell line at Day 12 (Table 3). An identical pattern of downregulation was shown by both Serpinh1, a collagen-specific molecular chaperone localized to the ER (Ito and Nagata, 2017), and vimentin (Vim), a type III intermediate filament protein responsible for maintaining cell shape and stabilizing cytoskeletal interactions (Musaelyan et al., 2018). Overall, Cd36, Hspa8, Serpinh1, and Vim were consistently downregulated in cell lines with higher Qp in our datasets in agreement with previous transcriptomics studies summarized in the meta-analysis (Tamo saitis and Smales, 2018).

Conclusion
In this work we present and make available to the community two RNA-Seq-derived transcriptomic datasets that comprehensively detail coding and non-coding transcript expression analysis of five IgG1-producing CHO cell lines and one null pool at different time points cultivated under fed-batch conditions in an ambrâ15 system. In particular, we use RNA-Seq to confirm the expression of lncRNAs in CHO cells and identify those whose expression is differentially regulated throughout fed-batch culture and between cell lines with different characteristics. The different time points for sample collection throughout culture and Qp of cell lines were clearly reflected in the PCA clustering of the transcript expression analysis and in the numbers of differentially expressed genes for the DAVI dataset. In addition, feeding was shown to be a significant source of variability even at early stage of culture, as shown by the comparison of seed train flask data used to inoculate the fed-batch process. These data show that there is a significant change in gene expression after only four days of culture among the same clone when the cells are still rapidly growing and dividing and before the major phase of production of the protein of interest. Taking into account this early transcriptional variability induced by changes in culture conditions, leading to differential regulation of the genome maintenance pathways, could be beneficial for a more effective cell engineering and clone selection in industrial settings.
KEGG functional enrichment analysis confirmed a tendency of pathways in the replication and repair domain to be differentially regulated in response to feeding when the seed train is inoculated in the fed-batch process and toward Day 12 in particular. The only exception was the clone with the highest Qp at Day 12, 3080, which interestingly did not show any differentially regulated pathways in the replication and repair domain. These data suggest that those cell lines that can maintain genome integrity and its surveillance may be better suited to prolonged culture and recombinant protein productivity. Although our datasets contained six different cell lines, we wanted to identify coding genes correlated to an increase in Qp across literature to improve robustness across different systems, leading to the identification of Cd36, Hspa8, Serpinh1, and Vim as genes negatively correlated with Qp in both our datasets and the most recent transcriptomics meta-analysis in CHO cells (Tamo saitis and Smales, 2018). Although the functions of those genes are heterogeneous, the conserved pattern of expression among very different experimental settings and cell lines suggests conserved roles with detrimental effects on Qp. Knock-out (KO) or knock-down (KD) strategies on the aforementioned genes with CRISPR or RNA interference could be implemented to investigate and validate any effects (Lee et al., 2015;Wu, 2009).
In conclusion, although many studies in CHO cells have investigated coding genes, our work aimed to unveil the non-coding transcriptome variation, specifically lncRNAs. We identified Adapt15, GAS5, and PVT1 among many others as lncRNAs linked to Qp. Although these genes lack an annotation in Chinese hamster genomes, their effects in model organisms and human diseases are well established (Colombo et al., 2015;Hollander et al., 1996;Ma et al., 2016). Although Adapt15 was initially identified in hamster and later linked to ER stress, DNA damage, and cell death, its role in CHO cells has never been further investigated.
Adapt15 upregulation in producing cell lines toward Day 12 of culture suggests an increasing stress on the ER as the recombinant protein is produced. In addition, we found higher levels of Adapt15 to be correlated to complete dysregulation of the genome maintenance pathways toward the end of culture in the producing clones. Our results suggest Adapt15 as a target for knock-down with RNAi or knock-out with CRISPR to alleviate its detrimental effects on productivity and genome integrity. GAS5 transcript accumulation occurs as a result of mTOR signaling repression, suggesting its use as a marker of translation repression and cell growth arrest in specific cell lines. Our data also suggest a correlation between the expression of the lncRNA PVT1 and of Myc; the relationship between the two depends on the stage of fed-batch culture and is reflected in the differential regulation of the PI3K/AKT pathway. When compared with an existent public database of gene expression (Singh et al., 2018), these lncRNAs are constitutively expressed in different CHO-S, CHOK1, and DG44 cell lines, suggesting a conserved functional role that could be exploited for cell engineering. Although lncRNAs investigation presents specific challenges in CHO due to the almost complete lack of annotation in the genome, our workflow allowed the identification, among many, of three strong candidates for cell engineering based on literature, structural predictions, transcriptional expression, and association with coding genes and pathways. Although we were able to propose functional mechanisms of action using secondary structure and RNA-protein interaction predictions, focused experimental studies on single transcripts will be required to assess their effects in mammalian cell factories.

Limitations of the Study
Although we have confirmed the expression of lncRNAs in CHO cells in this study by RNA-seq and qPCR and shown that the expression of some of these correlates with either cell growth or recombinant protein productivity, the conclusions around the potential impact of changes in the expression of specific lncRNAs that relate to productivity or growth across the different CGO cell lines are, at this stage, based upon correlation analysis only. Such correlation analysis does not imply cause or effect but merrily identifies these lncRNAs as potentially of interest. In order to determine any definitive link between the expression of specific lncRNAs and cell growth or productivity, it would be necessary to show experimentally that manipulation of the identified lncRNAs with expression profiles correlating with these phenotypes impacted the expected phenotype (growth or productivity) by, for example, knockdown/out or overexpression studies.

METHODS
All methods can be found in the accompanying Transparent Methods supplemental file.

DATA AND CODE AVAILABILITY
The accession number for the data reported in this paper is Gene Expression Omnibus (GEO) accession number GEO: GSE140671. The direct link to the data is www.ncbi.nlm.nih.gov/projects/geo/query/acc. cgi?acc=GSE140671.
The complete RNA-Seq dataset is thus provided as a reference for the community to investigate RNAseq data from CHO cells (see GEO accession number GEO: GSE140671 and Data S1).  Figure S2. Enriched KEGG pathways based on differentially expressed genes for each comparison between different cell lines at the same time point in the DAVI dataset. Each dot represents a pathway, with color shade representing the pvalue, size proportional to the overlap size (differentially expressed genes in the pathway) and x-coordinate recall (overlap size divided by the total number of genes in the pathway), Related to Figure 4. Figure S3. Enriched KEGG pathways based on differentially expressed genes for each comparison among different cell lines at day 4 and day 7 in the JCE dataset. Each dot represents a pathway, with color shade representing the p-value, size proportional to the overlap size (differentially expressed genes in the pathway) and x-coordinate recall (overlap size divided by the total number of genes in the pathway), Related to Figure 4. Figure S4. Enriched KEGG pathways based on differentially expressed genes for each comparison among the same cell line at day 4 against seed train culture in the JCE dataset. Each dot represents a pathway, with color shade representing the pvalue, size proportional to the overlap size (differentially expressed genes in the pathway) and x-coordinate recall (overlap size divided by the total number of genes in the pathway), Related to Figure 4.