Amplicon and Metagenomic Analysis of Middle East Respiratory Syndrome (MERS) Coronavirus and the Microbiome in Patients with Severe MERS

ABSTRACT Middle East respiratory syndrome coronavirus (MERS-CoV) is a zoonotic infection that emerged in the Middle East in 2012. Symptoms range from mild to severe and include both respiratory and gastrointestinal illnesses. The virus is mainly present in camel populations with occasional zoonotic spill over into humans. The severity of infection in humans is influenced by numerous factors, and similar to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), underlying health complications can play a major role. Currently, MERS-CoV and SARS-CoV-2 are coincident in the Middle East and thus a rapid way of sequencing MERS-CoV to derive genotype information for molecular epidemiology is needed. Additionally, complicating factors in MERS-CoV infections are coinfections that require clinical management. The ability to rapidly characterize these infections would be advantageous. To rapidly sequence MERS-CoV, an amplicon-based approach was developed and coupled to Oxford Nanopore long read length sequencing. This and a metagenomic approach were evaluated with clinical samples from patients with MERS. The data illustrated that whole-genome or near-whole-genome information on MERS-CoV could be rapidly obtained. This approach provided data on both consensus genomes and the presence of minor variants, including deletion mutants. The metagenomic analysis provided information of the background microbiome. The advantage of this approach is that insertions and deletions can be identified, which are the major drivers of genotype change in coronaviruses. IMPORTANCE Middle East respiratory syndrome coronavirus (MERS-CoV) emerged in late 2012 in Saudi Arabia. The virus is a serious threat to people not only in the Middle East but also in the world and has been detected in over 27 countries. MERS-CoV is spreading in the Middle East and neighboring countries, and approximately 35% of reported patients with this virus have died. This is the most severe coronavirus infection so far described. Saudi Arabia is a destination for many millions of people in the world who visit for religious purposes (Umrah and Hajj), and so it is a very vulnerable area, which imposes unique challenges for effective control of this epidemic. The significance of our study is that clinical samples from patients with MERS were used for rapid in-depth sequencing and metagenomic analysis using long read length sequencing.

ABSTRACT Middle East respiratory syndrome coronavirus (MERS-CoV) is a zoonotic infection that emerged in the Middle East in 2012. Symptoms range from mild to severe and include both respiratory and gastrointestinal illnesses. The virus is mainly present in camel populations with occasional zoonotic spill over into humans. The severity of infection in humans is influenced by numerous factors, and similar to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), underlying health complications can play a major role. Currently, MERS-CoV and SARS-CoV-2 are coincident in the Middle East and thus a rapid way of sequencing MERS-CoV to derive genotype information for molecular epidemiology is needed. Additionally, complicating factors in MERS-CoV infections are coinfections that require clinical management. The ability to rapidly characterize these infections would be advantageous. To rapidly sequence MERS-CoV, an amplicon-based approach was developed and coupled to Oxford Nanopore long read length sequencing. This and a metagenomic approach were evaluated with clinical samples from patients with MERS. The data illustrated that whole-genome or near-whole-genome information on MERS-CoV could be rapidly obtained. This approach provided data on both consensus genomes and the presence of minor variants, including deletion mutants. The metagenomic analysis provided information of the background microbiome. The advantage of this approach is that insertions and deletions can be identified, which are the major drivers of genotype change in coronaviruses. IMPORTANCE Middle East respiratory syndrome coronavirus (MERS-CoV) emerged in late 2012 in Saudi Arabia. The virus is a serious threat to people not only in the Middle East but also in the world and has been detected in over 27 countries. MERS-CoV is spreading in the Middle East and neighboring countries, and approximately 35% of reported patients with this virus have died. This is the most severe coronavirus infection so far described. Saudi Arabia is a destination for many millions of people in the world who visit for religious purposes (Umrah and Hajj), and so it is a very vulnerable area, which imposes unique challenges for effective control of this epidemic. The significance of our study is that clinical samples from patients with MERS were used for rapid in-depth sequencing and metagenomic analysis using long read length sequencing. KEYWORDS MERS-CoV, MinION, metagenomics, sequencing C oronaviruses were once described as the backwater of virology, as they did not cause extensive disease in humans. However, with the emergence of severe acute respiratory syndrome coronavirus (SARS-CoV) in China in 2003, Middle East respiratory syndrome coronavirus (MERS-CoV) in Saudi Arabia in 2012, and now SARS-CoV-2 originating in 2019 in China, this is clearly not the case. These viruses cause respiratory and gastrointestinal illnesses and share similar genome architectures and disease profiles. Patients with existing health care problems, such as underlying cardiovascular disease and genetic factors, may exhibit more severe symptoms and potentially have a fatal outcome (1)(2)(3)(4). As such, the presence of other respiratory pathogens may also be of critical importance in the etiology of these infections, and they may also be part of the normal healthy microbiome (5). Severe infection in humans is typified by cytokine storms (6,7), pneumonia, and kidney failure. MERS-CoV has been identified by the WHO as a prioritized disease and is on their list of pathogens for research and development in emergency contexts. With the advent of SARS-CoV-2, there is an urgent and unmet health care need to develop generic medical countermeasures to combat these infections and mitigate future outbreaks, particularly in contact tracing and understanding transmission dynamics. As a consequence, these viruses cripple national infrastructure, trigger city-wide transport curfews, and disrupt international travel and commerce. Many countries are at risk for both SARS-CoV-2 and MERS-CoV.
There are striking parallels between the emergence of all three of these viruses. For example, one of the major concerns with MERS-CoV is the potential spread of the virus and other pathogens during Hajj (8). This is an annual Islamic pilgrimage to Mecca in Saudi Arabia, involving some 2 million of the world's population and approximately 0.5 million Saudi residents. There are constant spillover events from camels into humans, and the geographical distribution of MERS-CoV in dromedaries is increasing either through spread and/or increased surveillance. Currently, antiviral therapies and vaccines have compassionate or emergency use for the prevention or treatment of SARS-CoV-2 infection in humans (9)(10)(11)(12); however, none exist for MERS-CoV. Several ongoing studies are evaluating medical countermeasures for MERS-coronavirus infection. They have included compounds already in use in the clinic, such as a combination of interferon-a2b and ribavirin (13). A phase 1 DNA vaccine, based on the viral spike glycoprotein, funded by US Department of the Army, has been recently evaluated in a human trial (14) and also in dromedary camels (15). Additionally, a phase 1 clinical trial has started in humans in Saudi Arabia assessing the safety and tolerability of a vaccine based on the ChAdOx1 MERS vaccine (16,17). The sporadic nature of the outbreaks and the lack of suitable animal models have hindered research.
Given the lack of medical countermeasures, shutting down transmission changes is essential to bringing such outbreaks under control (18), and sequencing the genomes of viruses can aid in this endeavor through molecular epidemiology, as was demonstrated with SARS-CoV in Singapore (19,20). Sequencing provides information on the rate of evolution of the virus (21), the origin of specific clusters (22), and whether an outbreak is caused by continuous human-to-human transmission (21,23); or it can be used to investigate zoonotic spillover (24). The use of long read length sequencing centered on the Oxford Nanopore platform has rapidly decreased sequencing turnaround time and can be applied to coronaviruses, for example with SARS-CoV-2 (25,26).
MERS-CoV is highly transmissible (27), and sporadic outbreaks are ongoing in the Kingdom of Saudi Arabia and surrounding regions. In order to aid in the diagnosis of MERS-CoV by rapidly generating a viral genome sequence and to map the microbiome in the nasopharyngeal sample, an amplicon-based and metagenomic MinION sequencing approach was used. During replication of coronaviruses, with an approximately 30kb positive-sense RNA genome, a nested set of subgenomic mRNAs are synthesized.
Subgenomic mRNAs located toward the 39 end of the genome are generally more abundant than transcripts nearer the 59 end (28,29), and these products can be identified by amplicon sequencing. This study used a mixture of amplicon sequencing and metagenomic approaches to characterize clinical samples from patients with MERS. The amplicon approach was based on generating 30, 15, or 8 overlapping fragments that spanned the genome but depended on the starting material. The rationale for using larger amplicons was to establish a system that could be used to identify deletions or recombination sites in the genome, which are hallmarks of coronavirus RNA synthesis. The data demonstrated that whole-genome or near-whole-genome information on MERS-CoV was rapidly obtained from samples taken from infected patients and sequenced using the amplicon approach. Furthermore, the ampliconbased sequencing approach provided data on both consensus genomes and the presence of minor variants, including deletion mutants. Metagenomic analysis provided information of the background microbiome and how this may be associated with outcome. Real-time analysis of the microbiome in patients and the identification of antibiotic resistance markers may provide better chemotherapeutic approaches to manage coinfections.

RESULTS AND DISCUSSION
Single nucleotide polymorphisms, recombination, and resulting deletions (and potentially insertions) may account for the wide genome diversity observed in some strains of coronaviruses (30)(31)(32). These recombination events can led to new strains of coronaviruses (33) or potentially affect vaccine strategies (34,35). In order to identify these genome changes with sequencing for MERS-CoV, rapid long read length approaches offered by Oxford Nanopore were developed. They used longer amplicons selected through appropriate primer pairs located along the MERS-CoV genome (Fig. 1). Note that this approach would also use the nested set of subgenomic mRNAs to capture sequence information.
Validation of primers and generation of amplicons using total RNA purified from MERS-CoV-infected cells. To evaluate the utility of the selected primers for the amplification of viral RNA under controlled conditions, total RNA was purified from MRC-5 cells that had been infected with the EMC strain of MERS-CoV at a multiplicity of infection (MOI) of 5. This RNA was used as a template to prime cDNA synthesis using random hexamers. Amplification conditions are provided in Table 1; the MERS-CoV genome was amplified using either 30 amplicons ( Fig. 2A), 15 amplicons (Fig. 2B), or 8 amplicons (Fig. 2C), using the same set of conditions for each set of amplicons (see Table S1 in the supplemental material). The rationale behind using the same amplification conditions across all primer pairs is that amplification would be more efficient if a large-scale sample analysis was required. These data indicated that for the all of the approaches, PCR products were observed that spanned the MERS-CoV genome. However, the 15-and 8-amplicon approach products varied in efficiency.
Generation of amplicons from patients infected with MERS-CoV and derivation of consensus genome sequence. The 30-, 15-, and 8-amplicon approaches were evaluated on RNA extracted from nasal aspirates taken from patients with MERS. The data indicated that the 30-and 15-amplicon approaches could be used to generate fragments from clinical samples ( Fig. 3A and B, respectively). The 8-amplicon approach was not sufficient for obtaining coverage across the MERS-CoV from a clinical sample (data not shown), probably due to the lower quality of RNA than that obtained from cell culture. The PCR products generated in the 30-and 15-amplicon approach (from separate patients) were combined for each patient, barcoded, and sequenced on separate flow cells for each patient. Sequencing reads generated by the MinION were aligned to a reference sequence. The analysis showed that a complete genome sequence could be obtained from the 30-amplicon (Fig. 4A) and 15-amplicon approach (Fig. 4B). A consensus sequence was generated using the ARTIC bioinformatic pipeline. To identify the viral genome sequence that dominated the population, a custom perl script was used to count the number of each nucleotide against the reference sequence (GenBank accession number NC_019843.3). This sequence was used to generate a bespoke consensus genome for an individual patient.
Analysis of the minor variant population within patients. The minor variant population in infections has been shown to influence the kinetics of virus replication and be associated with patient outcome (36). Therefore, methodologies were developed that could be used to assess the minor variant frequency within a sample from a patient infected with MERS-CoV. The custom perl script used to call the consensus also revealed the nucleotide depth and the counts of each nucleotide at each position (Fig. 5). The depth was used to normalize the mutation frequency into a proportion instead of a raw count, allowing a comparison of samples of different read depths. Nucleotides that had a count less than 20 were removed from analysis. As proof of principle, this approach was applied to the sequencing data obtained from patients 10 and 115, as these were of higher quality. Patient 10 appeared to have more base changes than patient 115 (Fig. 6). Transitions (A . G, G . A, C . U, and U . C) were more frequently observed, and C . U seemed more prominent than other mutations. This finding is consistent with RNA editing by APOBEC. APOBEC3 family members have been shown to be involved in restricting the growth of human coronavirus NL63 (HCoV-NL63) (37) and identified as potential drivers of C-to-U transitions in SARS-CoV-2 (38). We note that assessing minor variant populations using data from MinION sequencing may be problematic due to the higher error rate in sequencing than for example Illumina-based approaches (39). However, with sufficient read depth in a sample, an overview of the minor variant frequency/population can still be obtained. Adjusted gene size was calculated by dividing each MERS gene by the length of the total viral genome (30,108 bp). Global indel polymorphisms were grouped by category and mapped to their corresponding location within the reference sequence (Fig. 7). The mean proportion of each observation (against depth) was multiplied by the adjusted gene size to produce a standardized indel proportion for each gene (Fig. 8).
Identification and analysis of deletions in the viral genome in samples from patients. Deletions within the genome of MERS-CoV were identified in the sequencing data from patients 10 and 115 (Table 2). Patient 10 had five deletions in Orf1ab and a deletion spanning N. Patient 115 was sequenced with the 15-amplicon approach and therefore generated amplicons over 2 kb in length. A deletion of 77 bases that spanned the orf4b and orf5 genes was identified in this patient. Deletions in this region may have implications on virus pathogenesis, ORF4A is able to inhibit early antiviral responses (interferon a/b [IFN a/b]) in the host (40), and likewise ORF5 may reduce inflammatory responses (41). This finding was consistent in both bioinformatic pipelines that were used. The presence of these deletions does not rule out that at the minor variant level full intact gene were present. However, the data do illustrate the plasticity of the MERS-CoV genome.
Identification of bacterial and viral sequences in samples from patients and antibiotic resistance markers. The presence of other infections in patients with MERS may also require clinical management and influence the outcome. A metagenomic approach was used to identify other microorganisms present in the clinical samples from patients with MERS. Species such as Acinetobacter baumannii, Pseudomonas aeruginosa, and Streptococcus pneumoniae were identified, of which  all can be associated with bacterial pneumonia (Fig. 9). Although the detection of RNA does not always infer active infection, identification is based on the active bacterial transcriptome. Patients with Acinetobacter transcripts were associated with a fatal outcome of the disease (Fig. 9). Acinetobacter infections are often exclusive to health care settings, especially in patients who have received ventilation support. Acinetobacter species are considered a serious multidrug-resistant pathogen and are encompassed in the ESKAPE acronym, referring to Enterococcus faecium, Staphylococcus aureus, Klebsiella pneumoniae, Acinetobacter baumannii, Pseudomonas aeruginosa and Enterobacter species (42). Evidence for antibiotic-resistant genes was identified in patient 6 within Klebsiella pneumoniae (TEM-4) and Acinetobacter baumannii species (LpxC, adeI, adeJ, mexT, adeN, adeK, and ADC-2) ( Table 3). This patient died from a severe MERS infection. Several groups have suggested that the widespread use of antibiotics as part of the management of patients with COVID-19 may effect antimicrobial resistance (AMR) (43). Metagenomic analysis of samples from patients with severe coronavirus infection, as illustrated in this study, can be used to track these markers. Human alphaherpesvirus 1 transcripts were identified in patient 5 through this approach. Ongoing research throughout the SARS-CoV-2 pandemic suggests that coronaviruses may be able to reactivate herpes simplex virus and cytomegalovirus during severe disease (44).
The abundance of bacteria was compared between fatal cases and nonfatal cases using DESeq2 (Fig. 10). In general, transcripts from bacteria in the Proteobacteria phyla were higher in abundance in the fatal cases than those in nonfatal cases, including species from the Acinetobacter and Klebsiella genera. Previous studies have suggested interactions between viral infections and bacterial communities. An increase in the abundance of Proteobacteria was observed following an experimental rhinovirus infection in chronic obstructive pulmonary disease (COPD) patients (45). Patients with 2009 influenza A H1N1 infections and pneumonia were found to have an expansion of Proteobacteria compared with nonviral pneumonias (46)(47)(48).
The data indicated that an amplicon-based approach could be used to rapidly sequence MERS-CoV from clinical samples and provide information on genetic diversity and insertions/deletions. This approach was complemented with a metagenomic approach that was able to resolve the microbiome present in the clinical samples. This analysis identified bacteria associated with ventilation and also antibiotic resistance markers. Overall, the research demonstrates the utility of rapid long read length sequencing for characterizing MERS-CoV infection in samples from humans with MERS.

MATERIALS AND METHODS
Ethics statement. Ethical approval was obtained from the Institutional Review Board no. 18-102, King Fahad Medical City, Riyadh, Saudi Arabia. Samples for diagnostic purposes from nasopharyngeal aspirates (NPAs), oropharyngeal swabs, tracheal aspirates, and throat swabs were used for this study.
Sample collection and processing. NPAs, oropharyngeal swabs, tracheal aspirates, and throat swabs were collected from MERS-CoV-positive patients (45 to 74 years old) admitted to different hospitals within Saudi Arabia. MERS-CoV diagnosis was confirmed by reverse transcriptase PCR (RT-PCR) (bioMérieux Diagnostics). For this study, the NPAs with a confirmed MERS-CoV diagnosis from the Ministry of Health (MOH) Saudi Arabia, had no identifying information. The sampling NPAs, oropharyngeal swabs, tracheal aspirates, and throat swabs was carried out as per the MOH's guidelines. Samples were stored at 280°C until used. RNA from the NPAs was extracted using an EZ1 virus minikit v2 (955134; Qiagen). The RNA concentration was measured by the Qubit RNA broad-range (BR) assay (32852; Qiagen). Information on samples from patients used in this study is included in Table 4. This information includes sex, age, hospital, specimen type, threshold cycle (C T ) value of the E gene and ORF1AB, comorbidities, outcome, and whether the patient was in an intensive care unit.
Virus stock generation and infection. To prepare RNA from virus-infected cells as a control for amplification, the MRC-5 cell line was infected with MERS-CoV (Erasmus strain) at an MOI of 5. Total RNA was purified using the TRIzol method.
DNase Treatment of RNA. RNA was extracted from respiratory samples taken from 15 patients with severe MERS, including combined nasal and oropharyngeal swabs (n = 9), a throat swab (n = 1), nasal swabs (n = 4), and tracheal aspirate (n = 1). Metagenomic approaches have been used to identify the microbiome in samples from patients and correlate them with patient outcome (49). Sequence-independent, single-primer amplification (SISPA) was used to identify bacterial and viral transcripts present within the clinical samples. These transcripts were assessed using the real-time analysis pipelines EPI2ME with AMR or Kraken and Phyloseq packages. RNA samples were treated with Turbo DNase (Invitrogen, Vilnius, Lithuania) by adding 1 ml of Turbo DNase with 5-ml 10Â buffer in a 56-ml reaction. The reaction mixture was incubated at 37°C for 30 minutes. A total of 5 ml of inactivating agent was added and the mixture incubated at room temperature for 2 minutes. The reaction was centrifuged at 10,000 Â g for 90 s, and the RNA supernatant was transferred into a new tube.
Primer design and synthesis. Primers for the generation of overlapping amplicons were designed using the Primer3Plus platform and rechecked again by Primer BLAST (NCBI) to avoid primers with a high self-complementary score. Primers were synthesized by Eurofins Genomics (the sequence of the primers is shown Table S1). The stock concentration was 100 mM, and primers were diluted in DNase/ RNase-free H 2 O to make a 10 mM working concentration. Twenty MERS-CoV genome sequences were aligned with a reference sequence (NC_019843.3; the Erasmus Medical Centre [EMC] sequence). These sequences represented viruses collected in regions of Saudi Arabia and countries that reported cases, including South Korea. Primer binding sites were chosen from conserved regions after alignment so that  Table S1) that could be used to walk across the MERS-CoV genome with overlap between each generated amplicon. Alternatively, primer pairs could be selected that allowed the generation of larger amplicons that spanned more of the genome in contiguous reads (Fig. 1). To evaluate the utility of the selected primers for the amplification of viral RNA under controlled conditions, RNA was purified from MRC-5 cells that had been infected with the EMC strain of MERS-CoV at an MOI of 5. Infection was carried out under CL31 conditions and total RNA purified from these infected MRC-5 cells at 16 h postinfection. This RNA was used as a template to prime cDNA synthesis using random hexamers. Primer pairs (Table S1) were tested by gradient PCR to determine the optimum annealing temperature (data not shown).
cDNA synthesis and PCR. Superscript IV reverse transcriptase (18090010; ThermoFisher) and random hexamers (2.5 mM) were used to generate cDNA templates from RNA. cDNA was amplified using a Q5 high-fidelity DNA polymerase (M0491; New England BioLabs [NEB]).
Amplicon analysis and sequencing. PCR products were run on a 1% agarose gel in Tris-acetate-EDTA (TAE) buffer to confirm the presence of amplicons. The amplicons generated from an individual patient were pooled and cleaned using AMPure XP beads (A63882; Beckman Coulter) following the manufacturer's instructions before preparing the sequencing library with the ligation sequencing kit (SQK-LSK109; Oxford Nanopore Technologies). The sequencing library was added to a flow cell connected to a MinIT device, and sequencing was initiated through MinKNOW. An initial rapid base-calling setting was used.
SISPA. RNA was reverse transcribed with SuperScript IV reverse transcriptase (Invitrogen) using Sol-PrimerA (59-GTTTCCCACTGGAGGATA-N9-39), followed by second-strand DNA synthesis with Klenow  Takes into consideration the mapping quality scores, where a value greater than 10 has higher confidence.

Amplicon and Metagenomic Analysis of MERS-CoV
July/August 2021 Volume 6 Issue 4 e00219-21 msphere.asm.org 13 (NEB). Reaction conditions for round A were as follows: 1 ml of Sol-PrimerA (40 pmol/ml) was added to 4 ml of sample RNA, which was heated at 65°C for 5 min and then cooled on ice for 5 min. Then, 7 ml of SuperScript master mix (4 ml 5Â SuperScript IV buffer, 1 ml 100 mM dithiothreitol [DTT], and 1 ml RNaseOUT [Invitrogen]) and 1 ml SuperScript IV reverse transcriptase (200 U/ml) were added and incubated at 23°C for 10 minutes and then at 50°C for 10 minutes. The reaction was inactivated by incubating at 80°C for 10 minutes. For second-strand synthesis, Klenow fragment (NEB) was used. The cDNA from each reaction mixture was purified using a 1:1 ratio of AMPure XP beads (Beckman Coulter, USA). A total of 5 ml cDNA was amplified with Q5 high-fidelity DNA polymerase (NEB) according to manufacturer's instructions with Primer Sol-B. Cycling conditions were as follows: 98°C for 30 s; followed by 30 cycles of 98°C for 10 s, 55°C for 30 s, and 72°C for 1 min; with a final extension at 72°C for 10 min. All amplified cDNA was purified using a 1:1 ratio of AMPure XP beads (Beckman Coulter, USA) and resuspended in 25 ml. Purified cDNA was quantified using the Qubit double-stranded DNA (dsDNA) high-sensitivity (HS) assay (Q32851; Invitrogen) according to the manufacturer's instructions.
Bioinformatics. Fast5 files were base called again using Guppy (v.3.5.2) to mitigate for the lower accuracy associated with the rapid base calling setting initially used on the MinIT device. To investigate deletions in the MERS-CoV genome, library reads were filtered by the expected amplicon size and then aligned to the NCBI MERS reference NC_019843.3 using minimap2 (v.2.17-r941) as per the artic pipeline, and SVIM (v.1.4.2) was used to identify deletions (26).
In addition, SAMtools (v.1.10) was used to sort and index the alignment files, picard MarkDuplicates (v.2.23.4) was used to remove amplification duplicates, and then a custom perl script was used count the nucleotides at each position on the reference genome, when the mapping quality was above 10. This step ensures that every nucleotide on the derived consensus genome has been sequenced at least 10 times. Data were visualized using R studio using ggplot2. To consider the minor variation, the nucleotide depth at a single position with less than 20 counts was not taken forward into analysis to mitigate random sequencing errors.
For rapid assessment of the metatranscriptome approach, Fastq files were uploaded to the Oxford Nanopore Technology (ONT) cloud-based pipeline EPI2ME (Fastq antimicrobial resistance; WIMP [rev. 3.4.0], and ARMA CARD (rev. 1.1.6)) workflow to retrieve taxonomy classification and antimicrobial resistance information from MERS-CoV clinical samples.
In addition, metatranscriptomic reads were assessed with Kraken2 (v.2.2.1). Host removal was carried out to remove human reads. The human genome assembly GRCh38.p13 was used as a reference. The human genome assembly was indexed with minimap2 using the parameter "-ax map-ont" (50). The quality-controlled Oxford Nanopore fastq reads were mapped to the indexed human genome assembly with minimap2 using default parameters. The resulting SAM file was sorted and converted to bam with the command "samtools sort" (51). A fastq file with unmapped reads was produced with the "samtools fastq" command with the parameter "-f 4" (51).
The quality-controlled, host-removed reads were classified using Kraken2 (52). The standard Kraken2 output and sample report output files were produced with the options "-output" and "-report,"    (54). DESeq2 was used to compare the difference in abundance of species between fatal (n = 7) and nonfatal (n = 7) MERS-CoV infections (55).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. TABLE S1, DOCX file, 0.02 MB.

ACKNOWLEDGMENTS
This work was supported by the intramural research fund from Research Center, King Fahad Medical City, Saudi Arabia, under grant number 019-003 named "Elucidating the viral biology of MERS-CoV and the host response using high resolution sequencing" awarded to W.A. This work was also supported by the US Food and Drug Administration contract number 5F40120C00085 named "Characterization of severe coronavirus infection in humans and model systems for medical countermeasure development and evaluation" awarded to J.A.H.
We declare that we have no conflict of interest.