Transcriptome & viral growth analysis of SARS-CoV-2-infected Vero CCL-81 cells

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a positive-sense, single-stranded RNA virus, which has a genome of about 29.9 kb in length12. It belongs to the genus Betacoronavirus (subgenus: Sarbecovirus) of the family Coronaviridae3. The case-fatality ratio of SARS-CoV-2 is lower than that of SARS-CoV-1 and Middle East respiratory syndrome coronavirus (MERS-CoV)4. However, ever since its first report from China5, the SARS-CoV-2 has claimed far more deaths than both SARS-CoV-1 and MERS-CoV combined6. Full genome sequences of the first two SARS-CoV-2 viruses were reported from India in February 20207, and as on July 20, 2020, there were reports of 1,118,043 confirmed cases with 27,497 deaths from India8.

The SARS-CoV-2 genome encodes for non-structural proteins [open reading frames (ORF1a, and ORF1b)] and four structural [envelope (E), nucleocapsid (N), membrane (M) and spike (S)] proteins along with six accessory proteins from ORFs 3a, 6, 7a, 7b, 8 and 10129, the function of which is yet to be identified. The SARS-CoV employs viral RNA as the template for the replication and transcription of the virus. The expression analysis of SARS-CoV involves the synthesis of sub-genomic RNA by discontinuous transcription along with the genomic positive-sense RNA9. The use of next-generation sequencing (NGS) approaches has facilitated the analysis of transcription, as demonstrated in a study that revealed differential filovirus transcription patterns in different cell lines10.

A recent study by Kim et al11 looked upon the transcriptome and the epi-transcriptome of the SARS-CoV-2 through the sub-genomic RNAs, and the modification on the viral transcripts in SARS-CoV-2-infected Vero CCL-81 cells at 12 h post-infection (hpi). The purpose of the present study was to understand the growth kinetics of the SARS-CoV-2 and the variation in the transcription pattern of SARS-CoV-2 in a Vero CCL-81 cell line at different time points. This study also compared the quasispecies present in the clinical sample and the fourth passage (P) of the virus isolate.

Material & Methods

The study was conducted during April to May 2020 in the Maximum Containment Complex of the ICMR-National Institute of Virology (ICMR-NIV), Pune, India, with prior approval of the protocol.

Cell culture: SARS-CoV-2 used was isolated from the throat/nasal swab samples of Italian tourists (GISAID number: EPI_ISL_420545) in India during March 202012. The monolayers of Vero CCL-81 cells grown in 24-well tissue culture plates were infected with SARS-CoV-2 with a 50% tissue culture infective dose (TCID₅₀) of 10^6.5/ml, Vero CCL-81 P-3.

The study was performed with a multiplicity of infection (MOI) of 0.1 to infect the Vero CCL-81 cells in triplicate sets. After centrifugation, the supernatant and cell pellet were harvested at indicated time points (15, 30 and 45 min and 1, 2, 3, 4, 5, 6, 7, 8, 10 and 12 hpi) and stored at −80°C till further use. The supernatants treated with RNase A (10 μg/ml) for one hour at 37°C were used as control. Studies using infectious viruses were performed in the biosafety level 4 (BSL-4) facility of the ICMR-NIV. SARS-CoV-2 was inactivated using Trizol reagent (Thermo Fisher Scientific, USA) before removal to lower containment laboratories for further analysis. The TCID₅₀ was calculated for both the supernatant and the cell pellet. The cycle threshold (C_t) values for each sample were determined by E and RdRp-2 genes-specific SARS-CoV-2 real-time reverse transcription-polymerase chain reaction (RT-PCR) using the protocol described elsewhere713.

Next-generation sequencing and data analysis: Total RNA was extracted from 0.5 ml of the infected tissue culture lysate, containing 105 cells, and 1 ml of the cell supernatant using the TriPure reagent (Sigma-Aldrich, USA) followed by column purification with a Qiagen viral RNA extraction kit (Qiagen, Germany) as per manufacturer's instructions. Poly-A-containing mRNAs were purified from total RNA using oligo-dT beads provided in the Illumina TruSeq Stranded mRNA LT Sample Preparation Kit (Illumina, USA). The cDNA libraries were prepared from the purified mRNAs using the same library preparation kit. The libraries were sequenced using the Illumina MiniSeq platform (Illumina, USA) using the mid-output cartridge (300 cycles, 150×2 paired-end sequencing). The uninfected Vero CCL-81 cells (used as control) and the RNase A-treated control supernatants were processed using the above-described steps. RNase A-treated control supernatants were used after depleting ribosomal RNA (rRNA) from the sample with the NEBNext® rRNA Depletion Kit (New England Biolabs, USA).

Complementary DNA (cDNA) libraries were made from extracted RNAs using the same protocol as above. Total reads generated from the libraries made were mapped with the reference SARS-CoV-2 isolate Wuhan HU-1 genome (accession number: NC_045512.1). Transcription of each gene was quantified using the reads per kilobase million (RPKM) using the reads generated from each libraries. The RPKM values normalize the total reads for sequencing depth as well as gene length, which, in turn, gives the number of reads mapped to each gene. The read coverage and RPKM values were determined using the QIAGEN CLC Genomics Workbench 20.014 (QIAGEN, Aarhus, Denmark). RPKM values were calculated for each time point at which the infected cells were harvested to determine gene expression over time. Insertion and deletion variation analysis was performed with a maximum cut-off set to one per cent to identify variants present in the genes using the QIAGEN CLC Genomics Workbench 20.0. The methods used to determine the read coverage of the genome, RPKM values of the expressed genes, and the variants detected are limited by the algorithms implemented in the software used for the analysis.

Results

Copy number and growth kinetics: The viral RNA copy number in the cells and the supernatants harvested at defined intervals post-infection were measured using E and RdRp-2 genes of SARS-CoV-2-specific real-time RT-PCR. The average of the three replicates performed to determine the E and RdRp-2 genes copy numbers in the cell pellet, and the supernatant is depicted in Figure 1A and B. It was observed that the viral load was approximately constant for the first six hours in the cell pellet. At seven hours post-infection, E and RdRp-2 copy numbers increased by 10 log-fold further increasing exponentially till the 12 hpi. The copy numbers for E and RdRp-2 genes at the 12 hpi were found to be 4.1×10⁹ and 1.8×10⁸ RNA copies per 10⁵ cells, respectively. The RdRp-2 copies started appearing within the cell supernatants by 7 hpi, which increased 10-fold by the 8 hpi. Towards the 12 hpi, the viral RNA load further increased, and the final E and RdRp-2 copy numbers were 1.9×10⁷ and 6.6×10⁶ per 100 μl, respectively.

Close

Fig. 1

Copy number of E and RdRp-2 genes in the (A) cell pellet cell pellet and (B) supernatant of SARS-CoV-2 infected Vero CCL-81 cells: Cell pellet and supernatant of the SARS-CoV-2-infected Vero CCL-81 cells were collected at defined time points, and the viral load was determined using real-time RT-PCR. The average copy number±SD of three replicates is plotted.

Export to PPT

The number of infectious viral particles in the cell supernatant was measured by calculating the TCID₅₀ at different times post-infection, depicted in Figure 2. It was observed that the supernatant did not have any detectable level of virus particles until 6 hpi. Infectious viral progeny was observed in the supernatant with a titre of 1×10^1.5 TCID₅₀/100 μl starting from 7 hpi. An increased titre of approximately 1.5-fold (1×10^3.5 TCID₅₀/100 μl) was observed after 8 hpi, which increased until the 12 hpi (1×10^5.5 TCID₅₀/100 μl). NGS analysis of the cell supernatant showed a parallel increase of viral reads through the course of infection (data not shown).

Close

Fig. 2

50% tissue culture infective dose (TCID₅₀) of SARS-CoV-2-infected Vero CCL-81 cell line: Vero CCL-81 cells in the 24-well culture plates were infected with the SARS-CoV-2 strain (GISAID accession number EPI_ISL_420545). The supernatant was collected at defined times. The average of TCID₅₀/100 μl is obtained from the three replicates along with standard deviation as the error bars. The maroon line depicts the moving average of three periods.

Export to PPT

Next-generation sequencing and data analysis: The transcriptional data obtained for the different intervals were mapped against the SARS-CoV-2 isolate Wuhan HU-1. During the initial time points, only a few reads mapped to the N and ORF10 genes near the 3’ end of the genome. The mapped reads corresponding to the infected cells gradually increased towards the 5’ end. The expression of all viral genes was detected at 6 hpi in the infected cells. The supernatant of infected cells, which was used as a control in this experiment, did not show differential gene expression as that observed in the cell samples (Fig. 3).

Close

Fig. 3

Expression of SARS-CoV-2 virus genes in Vero CCL-81 cell line at a defined time point with respect to control. The CDS encoding different genes is depicted in the upper panel. The second panel determines the viral RNA present in the original supernatant used to infect the Vero CCL-81 cell line. The bottom two panels are the 1^st and the 6^th time point Vero CCL-81 expression data.

Export to PPT

Figure 4 depicts the RPKM values for structural, non-structural, and accessory genes at different times post-infection in the cell pellets. As expected, the cell pellets expressed a higher amount of N gene compared to any other genes. The presence of N gene transcripts at higher levels is suggestive of its importance in the viral replication and expression. The expressions of other genes are normalized with respect to the N gene so that their values lie between 0-1. The expression of N transcript was followed by the other two accessory proteins, ORF8 and ORF7a (Fig. 4B). The ORF1ab was the least expressed gene at all the periods post-infection (Fig. 4A). ORF7b transcripts were observed to be expressed towards the end of the replication cycle.

Close

Fig. 4

Reads per kilobase million (RPKM) values of cell pellet of (A) structural, (B) non-structural, and accessory transcripts of SARS-CoV-2 gene expressed in Vero CCL-81, collected at a defined time point. This value is normalized against the nucleocapsid (N) RPKM values at every point. S, Spike; E, envelope; M, membrane.

Export to PPT

Variant and quasispecies: The consensus nucleotide sequences from the clinical sample and isolate (P-0 and P-4) were compared, and a few common mutations were observed between them at 5’ UTR [genomic position (GP): 241]: ORF1ab (GP: 3037, 4809 and 14408) and S gene (GP: 23403). Nucleotide substitutions at GPs: 4809, 14408 and 23403 are non-synonymous changes causing mutation in the amino acid SàF, PàL, and DàG, respectively. The consensus nucleotide of the P-4 isolate had an additional mutation at the S gene (GP: 23607) compared to the P-1 and clinical sequence (P-0).

Quasispecies were observed in the sequencing data of the P-4 when compared with P-1 (GISAID number: EPI_ISL_420546) data. These nucleotide variations may have aroused among the virus population due to adaptation to host cells. Further, it was observed that the nucleotide reads at GP position 28254 (ORF8) of the P-4 isolate had a single-nucleotide deletion, leading to a frameshift (Fig. 5). The frameshift led to an increase in the length of the protein encoded by the ORF8 gene by four amino acids. Three of these quasispecies changes at GPs 8947, 27145 and 28902 are non-synonymous changes.

Close

Fig. 5

Per cent nucleotide variations observed in the SARS-CoV-2-infected Vero CCL-81 cell pellet retrieved after a 12 h time point of P-4. The variant was obtained using the basic variant tool present in the QIAGEN CLC Genomics Workbench 20.0. Dark blue colour shows the per cent of alleles present in the isolate sequences, whereas the light blue colour indicates the nucleotide observed in the reference sequence. The per cent is expressed as the variation in the length of the different nucleotides observed.

Export to PPT

Discussion

This study was aimed to analyze the growth kinetics and the transcription pattern of SARS-CoV-2 in the Vero CCL-81 cells. The determination of viral titres suggested that one complete virus replication took 7 hpi, while the expression of the complete gene set was observed at 6 hpi. The final concentration of the viral particle in the cell and supernatant was approximately two-fold more than the first release at 7 hpi. The TCID₅₀ value for 1 MOI demonstrated a similar trend (data not shown).

The ratio of gene expression was determined at 96 hpi in the infected cells with an MOI of 1 (data not shown) in the Vero CCL-81 cells. It was observed that M:E:S was expressed in the ratio of 0.22:0.09:0.08 with respect to N gene, whereas non-structural genes ORF8:ORF7a:ORF6:ORF3a had a ratio of 0.35:0.22:0.17:0.12 at the 7 hpi in the supernatant. The data revealed comparatively higher expression of ORF8, ORF7a, and M gene transcripts. ORF8 gene of SARS-CoV-2 is reported to diverge from the ORF8b gene of SARS-CoV, which is involved in the activation of innate immune response by inducing the NLRP3 (NOD-, LRR- and pyrin domain-containing protein 3) inflammasome15. The presence of a nucleotide deletion in the ORF8 transcript of the P-4 isolate reads was intriguing. The deletion in the ORF8 needs further investigation on its role played in eliciting host cell innate response and pathogenicity of the virus to the host cells. The role of ORF7a of SARS-CoV-2, which is 17.9 per cent diverged from the ORF7a gene of SARS-CoV, also needs further investigation. ORF7a of the SARS-CoV is reported to be a potent inhibitor of the bone marrow stromal antigen 2, which is involved in restricting the virus release from the cell membrane16.

Transcription of coronavirus occurs in both continuous and discontinuous methods17. Transcription regulatory sequences (TRSs) present at the 5’ end of each gene play an important role in the discontinuous method of transcription of sRNAs (sg mRNAs). The TRSs act to slow down the replication complex by controlling the transcription of sg mRNAs18. Alignment of the SARS-CoV-1 (accession number: NC_004718), SARS-CoV-2 (accession number: NC_045512), and MERS-CoV (accession number: NC_019843) reference sequences revealed the conservation of TRS leader and TRS sequences between different genes of SARS-CoV-1 and SARS-CoV-2. This possibly indicates a comparable level of transcriptional control for both SARS-CoV-1 and SARS-CoV-2. The higher expression of N transcript can be related to the different functions it plays in SARS-CoV-1 as RNA chaperone19, cell cycle regulation, pathogenesis, signal transduction, assembly, and stress20. The higher expression of an accessory protein, such as ORF8 and ORF7a, next to N, needs to be evaluated. Su et al21 demonstrated reduced replication fitness of the SARS-CoV-2 lined to the 382-nt deletion in the region spanning the ORF8 and its TRS.

The emergence of virulent or non-virulent subpopulation strains from the existing strains occurs in a process where a virus is under the influence of selection to adapt to a particular host22. The current study indicated the need for analyzing the different quasispecies present or developed during the passages of the virus in the cell line along with the changes in the clinical samples of a different host. The identification of quasispecies is important in a step towards designing drugs and vaccines for restricting the virus23 and understanding the wider range of variants present that need to be checked24.

The time period of growth kinetics needs to be enlarged to understand the kinetics of the virus which is one of the study limitations. The time point expression analysis of the subgenomic RNA needs to be looked upon to identify early expressing RNAs.

In conclusion, our findings for the growth kinetics indicate that SARS-CoV-2 completes its replication cycle in seven hours. Considering this observation, it could be hypothesized that the replication cycle of the SARS-CoV-2 would approximately be the same in humans. The quasispecies present in the clinical samples needs to be assessed to understand the changes in the virus population that occur due to host adaptation.