key: cord-0864831-2zon869b
authors: Konishi, Kazuhiro; Yamaji, Toshiyuki; Sakuma, Chisato; Kasai, Fumio; Endo, Toshinori; Kohara, Arihiro; Hanada, Kentaro; Osada, Naoki
title: Whole-Genome Sequencing of Vero E6 (VERO C1008) and Comparative Analysis of Four Vero Cell Sublines
date: 2022-03-22
journal: Front Genet
DOI: 10.3389/fgene.2022.801382
sha: dc415aa90ddcd568df6d902ba975d9fd2e880ad4
doc_id: 864831
cord_uid: 2zon869b

The Vero cell line is an immortalized cell line established from kidney epithelial cells of the African green monkey. A variety of Vero sublines have been developed and can be classified into four major cell lineages. In this study, we determined the whole-genome sequence of Vero E6 (VERO C1008), which is one of the most widely used cell lines for the proliferation and isolation of severe acute respiratory syndrome coronaviruses (SARS-CoVs), and performed comparative analysis among Vero JCRB0111, Vero CCL-81, Vero 76, and Vero E6. Analysis of the copy number changes and loss of heterozygosity revealed that these four sublines share a large deletion and loss of heterozygosity on chromosome 12, which harbors type I interferon and CDKN2 gene clusters. We identified a substantial number of genetic differences among the sublines including single nucleotide variants, indels, and copy number variations. The spectrum of single nucleotide variants indicated a close genetic relationship between Vero JCRB0111 and Vero CCL-81, and between Vero 76 and Vero E6, and a considerable genetic gap between the former two and the latter two lines. In contrast, we confirmed the pattern of genomic integration sites of simian endogenous retroviral sequences, which was consistent among the sublines. We identified subline-specific/enriched loss of function and missense variants, which potentially contribute to the differences in response to viral infection among the Vero sublines. In particular, we identified four genes (IL1RAP, TRIM25, RB1CC1, and ATG2A) that contained missense variants specific or enriched in Vero E6. In addition, we found that V739I variants of ACE2, which functions as the receptor for SARS-CoVs, were heterozygous in Vero JCRB0111, Vero CCL-81, and Vero 76; however, Vero E6 harbored only the allele with isoleucine, resulting from the loss of one of the X chromosomes.

Cell lines established from mammalian tissues are often used for virus isolation and culture as well as for vaccine production. One of the cell lines frequently used for these purposes is the Vero cell line, which is an immortalized cell line established from kidney epithelial cells of an African green monkey (Chlorocebus sabaeus, AGM) by Yoshihiro Yasumura at Chiba University in 1962. After the cell line was established, it was brought to the National Institute of Allergy and Infectious Diseases in the US, the American Type Culture Collection (ATCC), and the Japanese Collection of Research Bioresources (JCRB) cell bank (Figure 1 ). Various sublines have been established through passaging (Earley and Johnson, 1988; Mizusawa et al., 1988; Terasima et al., 1988) , which have slightly different properties from one another. For example, Vero CCL-81 (ATCC CCL-81), rather than Vero E6 (a.k.a. VERO C1008), is capable of more propagating Japanese encephalitis virus under prolonged culture conditions . Likewise, the subline Vero E6 likely propagates severe acute respiratory syndrome coronavirus (SARS-CoV) 2 more efficiently, compared with other Vero sublines (Harcourt et al., 2020; Matsuyama et al., 2020) and therefore Vero E6 and its derivatives have been widely used as the sublines for host cells in SARS-CoV-2 research (Basile et al., 2020; Ogando et al., 2020) . However, the genetic factors that contribute to these phenotypic differences are largely unknown.

Whole-genome sequencing analysis of one of the sublines, Vero JCRB0111 was performed by Osada et al. (Osada et al., 2014) and an approximate 9-Mbp homozygous deletion on chromosome 12 was identified. The deletion was found to contain a cluster of type 1 interferon (IFN-1) genes that act as viral suppressors, as well as CDKN2A and CDKN2B, which are involved in the cell cycle control. These results indicate potential factors that resulted in the immortalization of Vero cells and may partly explain why viruses can readily multiply in these cells. Furthermore, whole-genome resequencing analysis of additional sublines, Vero CCL-81 and Vero 76 identified numerous genetic variations among the sublines as well as integrations of Simian endogenous retroviral sequences (SERVs) (Sakuma et al., 2018) ; however, the effects of the variants at the nucleotide level were not evaluated. More recently, Sene et al. reported the haplotype-resolved genome assembly of the Vero CCL-81 subline (Sène et al., 2021) . Hundreds of protein-coding genes, including ACE2, were identified as losing function. AGM is susceptible to SARS-CoV-2 under laboratory conditions (Woolsey et al., 2021) . Although the intrinsic function of ACE2 is an angiotensinconverting, it is also known as the host cell receptor for SARS-CoV and SARS-CoV-2. Overall, these studies have provided information for quality control and have produced novel engineered sublines, which will accelerate the development of vaccine manufacturing platforms.

In the present study, we performed genome sequencing of Vero E6, which was established from a clone isolated from Vero 76 cells in earlier stage and conducted a comparative genome analysis among four major Vero cell sublines, currently available from public cell banks: Vero JCRB0111, Vero CCL-81, Vero 76, and Vero E6 (Figure 1 ). The identification of genetic variants specific or enriched in particular sublines will contribute to the elucidation of factors responsible for the phenotypic differences among the sublines.

Vero E6 cells were obtained from ATCC (ATCC-CRL-1586). The short-read sequences of Vero JCRB0111, Vero CCL-81, and Vero 76 were downloaded from a public database (DDBJ: PRJDB2865). Paired-end sequences of Vero E6 were determined using an Illumina HiSeq 2500. The reads were mapped to the AGM reference genome (NCBI: GCF_000409795.2) using the BWAmem algorithm (Li and Durbin, 2009 ). The single nucleotide variants (SNVs) and short indels (<50 bp) were called using VarScan 2 software (Koboldt et al., 2012) . The thresholds for SNV detection were: ≥13 coverage, ≥2 mutation count, ≥15 average base quality. The minimum coverage was selected to cover at least 95% of the genome for each sample. The strand filter was applied to exclude the sites where ≥90% of the reads mapped to one strand. Only bi-allelic SNVs were considered in this study. Manta was used to identify large-scale structural variations (Chen et al., 2015) .

The GTF format file downloaded from the Ensembl database was used for gene annotation (Chlorocebus_sabaeus.ChlSab1.1.86.gtf) and snpEFF software was used for annotating the effect of each genetic variant (Cingolani et al., 2012) . The assignment for the gene function categories was performed using DAVID (Jiao et al., 2012) . For evaluating the impact of missense variants, we used PROVEAN, SIFT, and PANTHER-PSEP (Choi et al., 2012; Choi and Chan, 2015; Tang and Thomas, 2016) . The RNA-seq experiment data of Vero E6 TMPRSS2+ (SRR13091741-SRR13091746) and COS-7 cells (SRR1919325-SRR1919327) were downloaded from the public database (Zhang et al., 2021) . The RNA-seq reads were mapped to the reference genome using HISAT2 using the default parameters (Kim et al., 2019) .

We defined subline-specific/enriched SNVs as variants with a significantly higher frequency compared with the other sublines. For each SNV, we counted the number of reads supporting the reference and alternative alleles and performed Tukey's test (Tukey, 1949) . SNVs with p < 0.05 were considered sublinespecific/enriched SNVs after controlling for multiple-testing using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995) .

We used Control-FREEC software to identify copy number variations (CNVs) and loss of heterozygosity (LOH) (Boeva et al., 2012) . A window size of 50 kbp and a step size of 10 kbp were selected. To quantitate the genetic differentiation between sublines, we computed f 2 statistics between all pairwise sublines. F 2 was calculated as the sum of the squares of allele frequency differences between two sublines, divided by the total number of SNVs. The f 2 values were used as a genetic distance between two sublines and a neighbor-joining tree was constructed using MEGA X (Knyaz et al., 2018).

PCR was performed as described previously (Sakuma et al., 2018) . Genomic DNA (30 ng) was used as the template. The primers specified below were synthesized (Eurofins Genomics Inc., Tokyo, Japan): Fw1: 5′-GGAACACCTGAAGATCTATGTGTC TA-3′, Rv1: 5′-ATCAAATTCCTCTCTTCACATCTTCT-3′, Fw2: 5′-GGACATATTGTTATAAAAGTTCATGG-3′, Rv2: 5′-GAAACTATACCTATGATTTTGCCATAG-3′

We obtained short reads from Vero E6 and mapped them to the AGM reference genome. To compare the genomic features of the Vero sublines, we reanalyzed previously published genomic data (Vero JCRB0111, Vero CCL-81, and Vero 76). In total, we obtained 74-, 32-, 43-, and 78-fold coverage data of the Vero JCRB0111, Vero CCL-81, Vero 76, and Vero E6 genomes, respectively. A 20% frequency cutoff yielded 9,141,259 SNVs, 453,797 short insertions, and 547,598 short deletions. As expected, 96% of these variants were shared among all sublines. The number of SNVs shared among sublines is summarized in Figure 2 .

The patterns of CNV and LOH are shown in Figure 3 , and the estimated copy numbers for each chromosome interval are provided in Supplementary Table S1 . Overall, the patterns appeared similar among the four sublines, including the 9-Mbp deletion on chromosome 12. However, there were several marked differences; for example, a large part of chromosome 15 had three copies in Vero JCRB0111 and Vero CCL-81 but the region was normal with respect to copy number in Vero 76 and Vero E6. Interestingly, chromosome 21 is highly rearranged in Vero JCRB0111, Vero CCL-81, and Vero 76, but the copy number of chromosome 21 in Vero E6 is normal, except for the distal end. Vero E6 also showed monosomy for the X chromosome. This feature is partially observed in Vero 76, showing a mosaic copy number for the X chromosome. This indicates that the Vero 76 cell line is composed of a heterogeneous cell population, which may be distinguished by one or two copies of the X chromosome.

Next, we evaluated the genetic relationship among the four sublines. We used f 2 statistics to measure the genetic distance between the sublines. F 2 statistics are statistics that measure the difference of the frequency of variants between two cell populations (Patterson et al., 2012) . The neighbor-joining tree reconstructed using f 2 statistics is shown in Figure 4 . The tree indicates the sister relationship between Vero JCRB0111 and Vero CCL-81, and between Vero 76 and Vero E6 with 100% bootstrap support. The reconstructed tree was incongruent with the inferred history of the Vero cell lineages (Figure 1 ). This may have resulted from rapid turnover of cell lineages within the sublines. The frequency of clonal cell lineages may change over time because of the difference in the survival/proliferation rate or other unknown mechanisms, and the shift would lead the differentiation of cell phenotypes (Heng and Heng, 2021) . The results suggest that different lineages co-existed before the split of Vero JCRB0111 and Vero CCL-81, which increased in frequency independently in Vero JCRB0111-CCL81 and Vero 76-E6. Further studies of single-cell genomic sequencing will reveal the dynamics of the cell populations during subline divergence.

The study of Sakuma et al. (Sakuma et al., 2018) revealed many SERV sequences that are present in the genomes of Vero JCRB0111, Vero CCL-81, and Vero 76, however, the pattern of insertions was perfectly consistent among the three sublines except for one insertion, which was referred to as SVL27b. SVL27b is present in the genome of Vero 76 but absent in the genomes of Vero JCRB0111 and Vero CCL-81 (Sakuma et al., 2018) . Considering that the AGM individual sequenced for draft genome harbored the insertion as a heterozygous state, Sakuma et al. (2018) concluded that the insertion was lost in Vero JCRB0111 and CCL-81 but retained in Vero 76. In this study, we investigated the pattern of SERV insertions in Vero E6, following the method of Sakuma et al. (2018) and confirmed that the integration pattern was the same as Vero 76, except for the one on the X chromosome. The SERV insertion around chrX:47012301-47012900 is heterozygous in Vero JCRB0111, Vero CCL-81, and Vero 76, but absent in Vero E6. The result is consistent with the CNV analysis, which showed that Vero E6 lost one of the two X chromosomes in almost every cell.

Whole-genome sequence analysis also showed that SVL27b is present in the Vero E6 genome, which is consistent with the idea that Vero E6 is a clonal derivative from Vero 76 cells (Earley and Johnson, 1988) . This was verified by genomic PCR experiments as shown in Figure 5 . When a PCR primer set matching a region~2-kbp away from the SVL27b inserted site, an approximate 3.7 kb DNA fragment corresponding to AGM chromosome 27, without the SERV insertion at this position, was amplified from all four Vero cell lines as well as the AGM control ( Figure 5 ). This indicates that the Vero cell lines have an allele(s) that do not contain SVL27b. Notably, a DNA fragment with an~7 kb SVL27b insertion was amplified from Vero E6 (Figure 5) . When a PCR primer was designed proximal to the SVL27b inserted site, a DNA fragment with the SVL27b insert was amplified from Vero 76 and Vero E6 cells, but not from the others ( Figure 5) . The DNA band amplified from Vero E6 was much denser than the band from Vero 76 ( Figure 5) . It should also be pointed out that, under the latter PCR conditions, amplified DNA without the SVL27b insertion was too short to be visible using agarose gel electrophoresis. Together with the results of the previous study (Sakuma et al., 2018) , these results strongly suggest that Vero 76 cells are a mixture of cell types with and without the SVL27b insertion and support that Vero E6 is derived from the Vero 76 lineage, which stably has SVL27b. Furthermore, these results indicate that SVL27b is a good genomic marker to distinguish the Vero 76-Vero E6 lineage from the Vero CCL-81 and JCRB0111 lineages (Figure 1 ).

We next looked into the genetic variants specifically observed or exhibiting a high frequency in each subline. We defined a sublinespecific/enriched variant as a variant with higher frequency Table S4) .

We further surveyed genes related to autophagy, apoptosis, antiviral activity, or innate immune response and harbored subline-specific/enriched missense or LOH variants in a single subline. For this purpose, we excluded the structural variants identified by Manta, because candidate variants could contain a nonnegligible number of false-positive variants and required additional experimental validations (Kawamoto et al., 2020) . In total, we identified eight genes with subline-specific/ enriched missense SNVs (Table 3) . On the other hand, none of the genes with subline-specific/enriched LOH variants have functions related to the above functional categories.

We focused on genetic variants specific/enriched in Vero E6, which efficiently propagate SARS-CoV-2. Missense variants in IL1RAP (L387F) and BR1CC1 (D696H) were almost exclusively found in Vero E6 in a heterozygous state, whereas those in TRIM25 (S418P) and ATG2A (E757Q) had a statistically high frequency in Vero E6.

IL1RAP is an auxiliary receptor of IL1R1, a receptor for interleukin (IL)-1α and IL-1β. Virus-induced cell death releases IL-1α in the cytoplasm (Malik and Kanneganti, 2018) , and the released IL-1α binds to IL1R1 on the plasma membrane. The intracellular Toll IL-1 receptor (TIR) domain of IL1R1 binds to IL1RAP, leading to the activation of key transcription factors and kinases associated with the inflammatory and immune response, such as NF-κB, AP1, JNK, MAPK, and ERK (Mantovani et al., 2019) . The leucine residue at 387 of the IL1RAP protein is highly conserved among mammals and prediction programs reported that the variant, leucine to phenylalanine, would be deleterious. Therefore, the L387F amino acid change in IL1RAP may disrupt the downstream cascade, suppressing the inflammatory and immune response and facilitating the increase in SARS-CoV-2 infection.

RB1CC1 constitutes a part of the ULK1 complex, which is required for the initiation of autophagy (Bello-Perez et al., 2020) . The ULK1 complex phosphorylates PI3KC3, a subunit of the phosphatidylinositol 3-kinase complex, which triggers the formation of a phagophore, and known as the sequestration membrane (Dikic and Elazar, 2018) . Phagophores envelop viruses and virus-derived antigens to form autophagosomes, and the fusion of lysosomes and autophagosomes results in the formation of autolysosomes and degradation of their contents (Ahmad et al., 2018) . Autophagosomes fuse with endosomes to form amphisomes, and amphisomes fuse with lysosomes to form autolysosomes (Zhao and Zhang, 2018) . The endosomes may contain SARS-CoV-2 that has entered via endocytosis (Bian and Li, 2021) . Previous studies have shown that ORF3a of SARS-CoV-2 inhibits two pathways that form autolysosomes, suggesting that it prevents itself from being degraded by inhibiting autophagy (Miao et al., 2021) . The prediction programs suggested that the D696H variant would be deleterious and deteriorate its function. Therefore, the occurrence of non-synonymous variants in RB1CC1 may have inhibited the initiation of autophagy and suppressed the degradation of SARS-CoV-2 by autophagy.

TRIM25 is a ubiquitin ligase that regulates RIG-1 (DDX58), which detects viral RNA and triggers the innate immune responses (Zeng et al., 2010) . TRIM25 causes polyubiquitination in a region of RIG-1 called CARD, which activates RIG-1 (Gack et al., 2007) . Once activated, RIG-1 recognizes viral RNA, it triggers the activation of NF-κB, IRF3, and IRF7, which induces the expression of the antiviral protein Viperin (Schneider et al., 2014) . Viperin catalyzes the production of ddhCTP, which interferes with RNA polymerase (RdRp), to promote the degradation of viral proteins, and interferes with the transport of viral proteins (Rivera-Serrano et al., 2020) . However, the amino acid site 418 in TRIM25 is relatively variable among other vertebrate orthologs and PROVEN software predicted that the variant of S418P would be functionally neutral.

ATG2A encodes a protein involved in lipid transfer between membranes (Osawa et al., 2019) and required for phagophore membrane expansion (Kotani et al., 2018) . Previous studies have shown that silencing ATG2A results in the accumulation of unclosed autophagosomes (Velikkakath et al., 2012) . The E757Q change in ATG2A was predicted as functionally tolerated by PROVEN and SFIT but possibly damaging by PANTHER score.

We investigated the genetic variants in the ACE2 gene of Vero E6. Sene et al. (2021) reported that the ACE2 gene in Vero CCL-81 potentially harbors structural variants causing LOF. However, using the publicly available RNA-seq data (Zhang et al., 2021) , we confirmed that ACE2 mRNA is properly expressed in the Vero E6 TMPRSS2+ cells both in the presence and absence of SARS-CoV-2. This indicates that ACE2 is expressed in Vero cells, although the angiotensin-converting enzymatic activity of the gene product is lost as shown by Sene et al. (2021) . We additionally investigated whether ACE2 is expressed in the COS-7 cell line, which is a kidney-derived permanent cell line of AGM and was established independently of the Vero cells (Jensen et al., 1964; Gluzman, 1981) . Unlike the Vero E6 cell line, a public RNA-seq data of COS-7 did not show any detectable gene expression of ACE2. We identified one missense SNV that showed marked differences between Vero E6 and the other sublines. The missense variants of V739I in ACE2 were heterozygous in Vero JCRB0111, Vero CCL-81, and Vero 76, but Vero E6 harbors only the isoleucine allele. The difference results from the fact that Vero E6 exhibits monosomy X.

In this study, we determined the whole-genome sequence of Vero E6, which has been widely used for the study of SARS-CoV-2 and performed comparative genomics on four different sublines that make up the major Vero cell lineages. Genomic resources for the Vero cell lines will benefit quality control of vaccine-producing cell substrates. In addition, finding candidates genes contributing to the different phenotypes of the cell lines will facilitate the identification of mechanisms of viral proliferation and the development of effective and safe substrates for vaccine production. The primary goal of this study was to present a whole-genome sequence of Vero E6 as research resources and catalog a list of candidate variants that potentially affect the phenotypic differences among the Vero sublines. The validation of each effect using additional sequencing and experiments will be necessary, although it is beyond the scope of the present study. Despite of these limitations, we provide a list of genetic differences among the four sublines, as well as variants specific or enriched in particular sublines, which represent a valuable resource for quality control of cell lines and understanding the mechanisms of viral proliferation.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ddbj.nig.ac. jp/, DRX311507.

Ethical review and approval was not required for the animal study because we used a cell line derived from AGM. The cell line is public and the ethical review is not applicable. 

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.801382/ full#supplementary-material 

Cell-based Culture Informs Infectivity and Safe De-isolation Assessments in Patients with Coronavirus Disease

Canonical and Noncanonical Autophagy as Potential Targets for COVID-19

Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing

Angiotensin-converting Enzyme 2 (ACE2): SARS-CoV-2 Receptor and RAS Modulator

Control-FREEC: A Tool for Assessing Copy Number and Allelic Content Using Next-Generation Sequencing Data

Manta: Rapid Detection of Structural Variants and Indels for Germline and Cancer Sequencing Applications

Provean Web Server: A Tool to Predict the Functional Effect of Amino Acid Substitutions and Indels

Predicting the Functional Effect of Amino Acid Substitutions and Indels

A Program for Annotating and Predicting the Effects of Single Nucleotide Polymorphisms

Mechanism and Medical Implications of Mammalian Autophagy

The Lineage of the Vero, Vero 76 and its Clone C1008 in the United States

TRIM25 RING-finger E3 Ubiquitin Ligase Is Essential for RIG-I-Mediated Antiviral Activity

SV40-transformed Simian Cells Support the Replication of Early SV40 Mutants

Severe Acute Respiratory Syndrome Coronavirus 2 from Patient with Coronavirus Disease, United States

Karyotype Coding: The Creation and Maintenance of System Information for Complexity and Biodiversity

Infection of Human and Simian Tissue Cultures with Rous Sarcoma Virus

DAVID-WS: A Stateful Web Service to Facilitate Gene/protein List Analysis

Identification of Characteristic Genomic Markers in Human Hepatoma HuH-7 and Huh7.5.1-8 Cell Lines

Graph-based Genome Alignment and Genotyping with HISAT2 and HISAT-Genotype

Varscan 2: Somatic Mutation and Copy Number Alteration Discovery in Cancer by Exome Sequencing

The Atg2-Atg18 Complex Tethers Pre-autophagosomal Membranes to the Endoplasmic Reticulum for Autophagosome Formation

MEGA X: Molecular Evolutionary Genetics Analysis across Computing Platforms

Fast and Accurate Short Read Alignment with Burrows-Wheeler Transform

Function and Regulation of IL-1α in Inflammatory Diseases and Cancer

Interleukin-1 and Related Cytokines in the Regulation of Inflammation and Immunity

Enhanced Isolation of SARS-CoV-2 by TMPRSS2-Expressing Cells

ORF3a of the COVID-19 Virus SARS-CoV-2 Blocks HOPS Complex-Mediated Assembly of the SNARE Complex Required for Autolysosome Formation

Cell Line Vero Deposited to Japanese Cancer Research Resources Bank

SARS-coronavirus-2 Replication in Vero E6 Cells: Replication Kinetics, Rapid Adaptation and Cytopathology

The Genome Landscape of the African green Monkey Kidney-Derived Vero Cell Line

Atg2 Mediates Direct Lipid Transfer Between Membranes for Autophagosome Formation

Viperin Reveals its True Function

Comparative Characterization of Flavivirus Production in Two Cell Lines: Human Hepatoma-Derived Huh7.5.1-8 and African green Monkey Kidney-Derived Vero

Novel Endogenous Simian Retroviral Integrations in Vero Cells: Implications for Quality Control of a Human Vaccine Cell Substrate

Interferon-stimulated Genes: A Complex Web of Host Defenses

Haplotype-resolved De Novo Assembly of the Vero Cell Line Genome

Panther-psep: Predicting Disease-Causing Genetic Variants Using Position-specific Evolutionary Preservation

History of Vero Cells in Japan

Comparing Individual Means in the Analysis of Variance

Mammalian Atg2 Proteins Are Essential for Autophagosome Formation and Important for Regulation of Size and Distribution of Lipid Droplets

Establishment of an African green Monkey Model for COVID-19 and protection against Re-infection

Reconstitution of the RIG-I Pathway Reveals a Signaling Role of Unanchored Polyubiquitin Chains in Innate Immunity

SARS-CoV-2 Hijacks Folate and One-Carbon Metabolism for Viral Replication

Autophagosome Maturation: An Epic Journey from the ER to Lysosomes

Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest