RNA-Seq of Nasopharyngeal Swabs from ICU Patients Infected with SARS-CoV-2 Omicron Sublineages in Adelaide, Australia, 2022–2023

1. Introduction

SARS-CoV-2, the virus that causes COVID-19, has evolved substantially since its emergence in late 2019 [1], with the omicron lineage emerging in 2022 and markedly diversifying into a series of sublineages [2,3,4,5]. Ongoing mutations in the spike protein are central to the classification of sublineages and are associated with increased transmissibility, immune evasion and receptor binding [6,7]. Other SARS-CoV-2 genes [8] are also undergoing evolutionary changes [9] that impact a range of virus-host interactions [10,11,12]. Omicron lineage viruses were characterized by increased replication in the upper (versus the lower) respiratory track [12] and a great capacity to infect ciliated epithelial cells [13,14] when compared with earlier lineage viruses. Targeting these cells facilitates viral transit across the mucus layer, as well as disrupting mucocilary clearance [15,16].

SARS-CoV-2 infection is generally diagnosed by RT-PCR or protein-based diagnostic tests performed on material collected from patients via nasopharyngeal swabs [17,18]. Such sampling collects not only viral RNA and protein, but also host cells and cellular debris from the nasopharyngeal track. Such swabs thus provide a readily available source of both viral RNA and human cellular mRNA from a major site of infection in COVID-19 patients [19]. We thus sought to determine what information might be gleaned by RNA-Seq analysis of nasopharyngeal swabs, correlating the data with patient information, with the aim of increasing the understanding of nasopharyngeal infections by omicron sublineages in humans.

2. Materials and Methods

2.1. Ethics and Regulatory Statement

Informed consent was obtained from all patients involved in the study. Nasopharyngeal swabs were collected from 49 COVID-19 patients admitted to the ICU of the Royal Adelaide Hospital from October 2022 to August 2023. Swabs were placed into TRIzol® reagent (Invitrogen) and stored at -80°C prior to transport to QIMR Berghofer. Approval was obtained from the Central Adelaide Local Health Network (CALHN) Human Research Ethics Committee (CALHN13050). Approval was also obtained from the QIMR Berghofer Human Research Ethics Committee (Ref: P3600), with all samples and patient data deidentified prior to being sent to QIMR Berghofer. Research at QIMR Berghofer was approved by the Institutional Safety Committee.

2.2. RNA-Seq and Bioinformatics

RNA was extracted from TRIzol according to manufacturer’s instructions. DNAse I treatment (New England Biolabs) and RNA clean-up was undertaken using RNeasy MinElute Cleanup Kit (QIAGEN), according to the manufacturer's instructions. The Agilent 4200 Tape Station was used to provide RIN scores. Library preparations and ribosomal RNA depletions were undertaken using TrueSeq Stranded Total RNA Sample Preparation Kit, with Ribo-Zero (Illumina). Tape Station D1000 Screen Tape was used to assess the fragment size of the completed libraries; samples where cDNA amplicons could not be detected did not progress to sequencing. ThermoFisher Qubit 4 Fluorometer was used to measure the concentration of the completed libraries. Sequencing was performed using NextSeq 2000 using the flow cell P3 Reagent (200 cycles), providing 76 bp paired-end reads. All samples were run in a single sequencing run; Q30 was 89.97%.

STAR aligner was used to align the RNA-Seq reads to a combined SARS-CoV-2 BA.5 (GenBank: OP604184.1) and human (GRCh38, version 38) reference genome [20]. Samtools v1.16 [21] was used to generate viral read counts expressed as viral reads per million reads (RPM). RSEM v1.3.1 was used to generate read counts for host genes; differentially expressed genes (DEGs) were identified using DESeq2 v1.40.2 [22]. DEGs were analyzed using Ingenuity Pathway Analysis (IPA, v84978992) (QIAGEN). Relative abundance of specific cell types was estimated via cellular deconvolution using the SpatialDecon package in R [23,24], with cell-type expression matrices obtained from the human lung cell atlas [25] (https://data.humancellatlas.org/explore/projects/c4077b3c-5c98-4d26-a614-246d12c2e5d7).

2.3. RT-qPCR of Nasopharyngeal Samples

cDNA synthesis (from the purified RNA from the nasopharyngeal swab samples) was undertaken using ProtoScript II First Strand cDNA Synthesis Kit (New England Biolabs). RT-qPCR was performed using iTaq Universal Probes Supermix (Bio-Rad) and SARS-CoV-2 (Sarbeco-E) primers; Forward: 5′-ACAGGTACGTTAATAGTTAATAGCGT-3′, Reverse: 5′-ATATTGCAGCAGTACG CACACA -3’ [26]. The florescent probe was 5’-FAM-ACACTAGCCATCCTTACTGCGCTTCG- ZEN™/Iowa Black® FQ-3’ (Integrated DNA Technologies Australia). BioRad CFX96 was used to perform RT-qPCR with the following cycling conditions: incubation at 50°C for 10 min, followed by 3 min at 95°C, then 40 cycles of 15 secs at 95°C and 30 secs at 60°C.

2.4. Phylogenetic Tree Generation

The Integrative Genomics Viewer (IGV; v.2.9.4) [27] was used to align the viral sequences retrieved from RNA-Seq to the original SARS-CoV-2 strain. Consensus spike sequences were extracted from IGV, followed by translation via Expasy (online tool: https://web.expasy.org/translate/). MegaX (v.11) was used to generate the phylogenetic tree. Sublineages were classified further based on their spike mutations using GISAID.

2.5. Statistics

Correlations were undertaken using SPSS Statistics (v23) (IBM). Fischer’s exact tests were performed using JNP Pro (SAS Institute, Cary, North Carolina, USA).

3. Results

3.1. Omicron Sublineages and Viral Read Counts from COVID-19 ICU Patients

Human nasopharyngeal swab samples were obtained from consented COVID-19 patients admitted to the Intensive Care Unit (ICU) at the Royal Adelaide Hospital between October 2022 and August 2023. The virus sublineages from patients were classified by the diagnostic laboratory at the hospital using RT-PCR of nasopharyngeal swab samples, with viral infections grouped into four omicron sublineages (Figure 1a, x axis).

Nasopharyngeal swab samples from 49 patients were collected into TRIzol® reagent and were transported to QIMR Berghofer. RNA was isolated, with 38 patient samples yielding sufficient RNA (>100 ng) for library preparation and RNA-Seq (Table S1). Viral read counts, expressed as viral reads per million (RPM), for this patient cohort ranged across >5 logs, with levels >5000 viral RPM only seen in patients infected with BA.5 and BA.2-like sublineage viruses (Figure 1a, yellow box). These higher viral loads were associated with age > 64 years (see below).

RNA-Seq of 19/38 patient samples yielded sufficient spike sequence data to confirm the prior RT-PCR-based sublineage identification (Figure 1a, black circle outlines). The variant or sublineage for four samples could not be established by either RNA-Seq (insufficient reads) or by the hospital diagnostic laboratory (Figure 1a, Unknown).

3.2. Quantitation by RT-qPCR Correlated Well with Viral read Counts

The same purified RNA samples that were used for RNA-Seq were also analyzed by quantitative reverse transcription polymerase chain reaction (RT-qPCR) using the standard SARS-CoV-2 Sarbeco-E primers [26]. Viral RPM obtained from RNA-Seq correlated well with Cq values from RT-qPCR; viral RNA quantitation by these two methods thus provided highly comparable results.

3.3. Phylogenetic Tree of Omicron Sublineages

Of the 38 samples analyzed by RNA-Seq, 19 provided sufficient spike sequence data to allow construction of a phylogenetic tree (Figure 1c). Clear clustering of the sublineages was evident, reflecting the groupings applied in Figure 1a. As might be expected, XBF, a recombinant of BA.5.2 and BA.2.75 [28], grouped with the BA.2-like viruses.

3.4. Progression of Omicron Sublineages from Oct 2022 to Aug 2023

The number of patients infected with an identified sublineage was tracked month by month from Oct 2022 to Aug 2023, and is illustrated as a percentage of each sublineage per month (Figure 1d). A clear progression was evident from BA.5 to BA.2-like to XBB-like and finally XBC (Figure 1d), consistent with the global evolution of omicron sublineages [2,3,4,5,29].

3.5. Correlations of Viral Read Counts with Age and Other Patient Data

A range of patient data were collected and potential correlates sought. High viral RPM were significantly associated with patients who were > 64 years of age (Figure 2a), in agreement with previous reports [30,31]. Corticosteroid and baricitinib (JAK inhibitor) treatment showed no significant effect on viral RPM. Corticosteroid treatment is well established for reducing mortality [32], although a delay in viral clearance has been reported in some [33], but not other studies [34]. Baricitinib treatment also reduces mortality [35] and has been reported not to affect viral clearance [36]. All patients were treated with Remdesivir (SARS-CoV-2 antiviral).

A series of other correlations failed to show significance. Whether the patients had received a COVID-19 vaccine did not significantly influence the mean viral RPM; however, only 3 patients did not receive a vaccine (Figure S1a). High viral RPM were associated with patients that had one or more comorbidities; however this did not reach significance, with only 5 patients free of comorbidities (Figure S1b). Previous studies have shown that patients with high viral loads often have comorbidities [37]. The length of stay in ICU or in hospital also did not correlate with viral RPM (Figure S1c), consistent with previous studies [38]. The viral RPM were not significantly different between males and females (Figure S1d). Viral loads are generally reported to be higher in females [39].

3.5. RNA-Seq Analysis of Human Gene Expression

There was sufficient RNA (>100 ng) in 38/49 samples for RNA-Seq. RNA integrity number (RIN) scores ranged from 2.3 to 7.8 (mean 5.4 + SD 1.4), illustrating that RNA in the nasopharynx, as might be expected, often suffers from a level of degradation. This is a recognized issue in many clinical samples, and can be ameliorated by depletion of ribosomal RNA, rather than enriching by poly(A) capture, prior to library generation [40,41]. Approximately 30 million reads were obtained for each sample with a mean of 75% aligning to the human reference genome (Table S1), illustrating that human gene expression data can be obtained from most nasopharyngeal swabs.

DeSeq2 was used to identify differentially expressed genes (DEGs) using pairwise comparisons between samples from patients infected with the different sublineages. This approach provided only a small number of DEGs for each comparison, with BA.5 vs. XBC providing the largest number of DEGs (n=43 DEGs with gene name annotations) (Table S2). No cogent pathways were identified by inter alia Ingenuity Pathway analysis (IPA) [23,42,43]. However, the top upregulated gene was Growth and differentiation factor 15 (GDF15) (Table S2) (a member of the TGF-β superfamily), which in COVID-19, increases with tissue damage [44]. As the mean viral RPM was ≈17 times higher in BA.5 compared with XBC (Figure 1a), more tissue damage might be envisioned in the BA.5-infected cohort.

To mitigate against the influence of large differences in viral loads, patient samples with comparable viral RPM were compared for human gene expression (Figure 3a, yellow shaded box). The exclusion of samples with high and very low viral RPM, left four XBC and four BA.5 samples (Figure 3a). These 8 samples had slightly higher RIN scores (range 3.8 - 7.8, mean 5.6 + SD 1.2). DeSeq2 again only identified a small number of DEGs; however 9/10 of the mRNA species that were down-regulated in XBC samples were associated with cilia (Figure 3b). This analysis argues that more ciliated epithelial cells were infected and disrupted in the nasopharynx of XBC-infected patients than in BA.5-infected patients. The genes were the cilia-associated protein CFAP95 [45] and CFAP65 (expressed in lung [46]), ROPN1L [47], MAPK8IP1 [48], RAB36 (paralog of RAB34) [49], CES1 (highly expressed in ciliated epithelium [50]), LDLRAD1 [51], UBXN10 [52] and UNC119 [53].

Using the data from the same samples (Figure 3a, yellow shading), but using the entire expression matrices (rather than just DEGs), cellular deconvolution analysis was undertaken to estimate the relative abundance of different cell types. Cell-type expression matrices for cells from the human nasopharynx (upper respiratory track) are not available, so cell-type expression matrices for human lung (lower respiratory track) were used. The abundance of ciliated cells emerged as significantly lower in XBC vs. BA.5 samples (Figure 3c), consistent with the analysis of individual DEGs in Figure 3b. Of all the cell types identified by cellular deconvolution, only the abundance of ciliated cells was significantly different between XBC and BA.5 samples (Figure 3d).

4. Discussion

Herein we present viral and human gene expression data derived from nasopharyngeal swabs using RNA-Seq. Ribosomal depletion rather than poly(A) capture was used to ameliorate the problem of lower RNA quality (lower RIN scores, indicating some RNA degradation) [35,36]. The viral RPM correlated well with standard RT-qPCR viral load quantitation, providing a level of cross-validation of RT-qPCR and RNA-Seq data. High viral RPM were significantly associated with age >64, and patients with high viral RPM also often had comorbidities, data consistent with previous studies [30,31,37]. Immunosuppression with corticosteroids and/or baricitinib did not significantly influence viral RPM, results again consistent with previous studies [34,36]. Taken together, these data illustrate that viral RPM obtained via RNA-Seq provides reliable information on viral loads in the upper respiratory track. However, in this context RT-qPCR is currently clearly cheaper and more accessible and unlikely to be replaced by RNA-Seq.

The advantage of RNA-Seq is that it also provides gene expression data for the human host, with the data presented herein illustrating that a very respectable number of human mRNA reads can be obtained from nasopharyngeal swabs (Table S1). The very large range in viral loads may be due to a number of factors including sampling at different times during the course of the infections, differing effectiveness of anti-viral treatment or anti-viral immunity, and/or differing levels of the initial infectious viral inoculum. Although perhaps to be expected in any cohort of virus-infected patients, such large variability complicates the ability to generate salient comparisons of human gene expression between samples from patients infected with different sublineages. Large differences in the level of virus infection will have a much greater effect on gene expression than the more subtle differences that might emerge between infections with different subvariants of the same virus. Nevertheless, we illustrate that choosing patient samples with comparable viral RPM, allows generation of human differential gene expression data that reveals information about the biology of omicron sublineages. Specifically, that XBC infection in the nasopharynx resulted in significantly lower ciliated epithelial cell mRNA signatures when compared with BA.5 infection. Conceivably, the reduction in ciliated cell signatures might be associated, at least in part, with reprograming of gene expression in the infected cells [13]. However, XBC may simply infect and ultimately kill ciliated epithelial cells (via virus-induced cytopathic effects [54,55]) more efficiently than BA.5, leaving fewer of these cells available to be collected by the nasopharyngeal swabs.

5. Conclusions

RNA-Seq of human nasopharyngeal swabs provided evidence that the omicron sublineage XBC targets ciliated epithelial cells in the human upper respiratory track more effectively than BA.5. A key difference between the omicron lineage viruses and previous SARS-CoV-2 lineages was increased infection of ciliated epithelial cells in the upper respiratory track [13,14]. The data presented herein suggest that this evolutionary trend of increased targeting of ciliated cells continued as the omicron sublineage evolved from BA.5 to XBC.