Use of whole-genome sequencing in the molecular investigation of care-associated HCoV-OC43 infections in a hematopoietic stem cell transplant unit.

BACKGROUND: While respiratory viral infections are recognized as a frequent cause of illness in hematopoietic stem cell transplantation (HSCT) recipients, HCoV-OC43 infections have rarely been investigated as healthcare-associated infections in this population.
OBJECTIVES: In this report, HCoV-OC43 isolates collected from HSCT patients were retrospectively characterized to identify potential clusters of infection that may stand for a hospital transmission. STUDY
DESIGN: Whole-genome and S gene sequences were obtained from nasal swabs using next-generation sequencing and phylogenetic trees were constructed. Similar identity matrix and determination of the most common ancestor were used to compare clusters of patient's sequences. Amino acids substitutions were analysed.
RESULTS: Genotypes B, E, F and G were identified. Two clusters of patients were defined from chronological data and phylogenetic trees. Analyses of amino acids substitutions of the S protein sequences identified substitutions specific for genotype F strains circulating among European people.
CONCLUSIONS: HCoV-OC43 may be implicated in healthcare-associated infections.

Background

Although highly pathogenic human coronaviruses (HCoV), like SARS-CoV and MERS-CoV, have emerged recently, the majority of HCoV infections are caused by HCoV-229E, HCoV-NL63, HCoV−OC43 and HCoV-HKU1. These coronaviruses are mainly recognized as causative agents of community-acquired infections. With the exception of SARS and MERS-CoV, there are not many published studies about health-care associated HCoV infections [[1], [2], [3], [4], [5]]. HCoV-229E, HCoV-NL63, HCoV−OC43 and HCoV-HKU1 are usually associated with mild diseases in immunocompetent patients, but can cause severe respiratory tract infections in fragile populations. Indeed, development of molecular detection tests for diagnosis has shown that HCoV are clearly responsible for severe infections in immunocompromised patients, including hematopoietic stem cell transplant recipients (HSCT) [[6], [7], [8]].

HCoV−OC43, HCoV-229E, HCoV-NL63 and HCoV-HKU1 account for 6.7–15.4% of the viruses detected during respiratory infections in HSCT recipients [9]. Currently, HCoV−OC43 is detected around the world, circulating all year long with a slight predominance during winter in temperate countries. The efficient transmission via small droplets and a prolonged shedding by HSCT recipient’s contribute to virus dissemination and highlight the need for a better understanding of hospital-acquired HCoV infections [[10], [11], [12], [13]]. Although an outbreak of HCoV-NL63 respiratory infections in a long-term care facility has recently been reported, a limited amount data are available on the overall genomic characteristics of HCoV−OC43 healthcare-associated infections [14].

The main objective of this work was to identify potential clusters of infection that may stand for a hospital transmission. Next-Generation Sequencing (NGS) was used to generate whole-genome sequences from nasal swabs and phylogenetic analyses were then applied to characterize sequences. The secondary objective was to identify molecular signature of circulating strains.

Methods

Patients

Eight HSCT recipients with HCoV−OC43 detected in nasal swabs by multiplex respiratory viral PCR using Anyplex II RVS 16 detection kit at the University hospital of Lille between 2013 and 2015 were included in the study. Demographic and clinical data were collected. All patients were informed that they were included in a research cohort and agreed that their biological samples could be used for research purposes. Nasal swabs studied here were included into the biological collection of the University Hospital of Lille declared to the French Ministry of High Education and Research (reference number DC-2008-642). According to the French laws, this declaration implies acceptance by an ethic committee. Patients are named as MDSX for “Maladies Du Sang”X.

Sequencing

Viral genome sequencing was performed as described by Maurier et al. [15]. Briefly, extracted viral RNA was reverse transcribed, PCR amplified using multiple primer pairs and pooled before preparing sequencing libraries using NxSeq AmpFREE LowDNA (Lucigen) library kit and barcoded with Illumina-compatible adaptors. Libraries were paired-end sequenced in 2 × 300 cycles on Illumina MiSeq. Sequencing analysis and genome assembly was performed as described by Maurier et al. [15]

Genomes annotation

Annotation of the 8 newly sequenced genomes MDS2, MDS4, MDS6, MDS11, MDS12, MDS14, MDS15 and MDS16, was performed by comparison with the HCoV−OC43 isolate from Mexico (Genbank accession KX344031) [16] using BlastN [17]. This isolate served as a reference and allowed us to compute the coordinates of the following functional elements: ORF1a, ORF1b, ns2A, HE, S, M and N genes.

Phylogenetic analysis

For the whole genome evolutionary analysis of MDS2, MDS4, MDS6, MDS11, MDS12, MDS14, MDS15 and MDS16, we selected 40 published full-length human coronavirus genomes with diverse genotypes and geographical origins, making a total of 48 sequences. We also considered more specifically the phylogeny of S gene and included 12 additional sequences from partially sequenced genomes from Caen hospital (France) and one additional published sequence from China (Table supplementary data), making a total of 61 sequences. For each of the two datasets, three distantly related bovine coronaviruses served as outgroup and multiple-sequence alignment was built using MUSCLE from SeaView 4.6.4 followed by a manual correction step taking into account protein coding sequences for coding regions [18,19]. Based on the multiple sequence alignments, similarity between each pair of nucleotide sequences was computed using the Sequence Manipulation Suite [20]. Phylogenetic trees were then inferred using the maximum likelihood method implemented in PhyML 3.1 with a general time reversible (GTR) nucleotide substitution model, a proportion of invariant sites (+I) and a gamma rate heterogeneity (+Γ) with 1000 bootstrap replicates [21]. The substitution model was chosen by optimizing the AIC score using smart model selection [22]. A total of 30,815 nucleotides positions were included for the whole genome analysis and 4137 nucleotides positions for S gene. The software FigTree 1.4.4 was used to produce figures [23].

Substitution analysis

Amino acids comparative analyses for ORF1a, ORF1b, ns2A, HE, S, M and N proteins were done on 58 sequences strains (8 MDS + 40 previously included and 10 others) for which complete genomes were available and on 76 sequences for S gene only (8 MDS + 68 others) (Table Supplementary data). For this analysis, we added new sequences compared to the phylogenetic analysis. Indeed, we need to examine each position of the sequences separately, which requires a larger support. By doing this, we reduced the chance to select erroneous mutations. The goal is to ensure that the informative mutations found are characteristic of a particular genotype: they are found only in this genotype and are shared by all sequences known for this genotype. In both cases, the signature amino acid substitutions were determined by selecting a subset of informative sites using DIVEIN [24]. Colorization was performed manually to highlight conserved sites.

Spatio-temporal and evolutionary dynamics

We followed the methodology used in [25]. Notably, we added Asiatic sequences that were also used in [25] and [26]. We estimated divergence times using BEAST 1.8.4 [12] on both whole genome and S gene [27]. Exponential population size, relaxed clock with coalescent tree (BSP distribution for S and uncorrelated exponential distribution for whole genome) was tested to estimate the time of most recent common ancestor (tMRCA).

All information on the sequences used (identifiers, strains, bibliographical references) is summed up in the supplementary Table S1.

Complete genome sequences were deposited in Genbank under accession numbers MK303619 to MK303625, MK327281

Results

Patients and viral characteristics

HCoV−OC43 isolates were detected by RT-PCR from nasal swabs collected from 8 HSCT recipients. Table 1 shows the characteristics of these patients. The median time to HCoV−OC43 infection after HSCT was 239.5 days (range 13–1052 days). All patients had graft-versus-host diseases for which they received a steroid daily dose of higher than 1 mg/kg, except patient MDS4. All patients died except one, MDS4; however, in all the cases, death was not directly related to HCoV−OC43 infection.

Table 1

Clinical features of patients with HCoV−OC43 infections.

	Age	Sex	Hematology feature	Transplant numbers	Donor type	Days between transplant and HCoV-infection	date of HCoV detection	Cq value	Respiratory or infectious symptoms	Platelet count	Leukocyte count	Neutrophil count	Lymphocyte count	Monocyte count	steroids	GVH	Co-detection	Outcome
	y						m/d/y			x10exp9/L	x10exp9 /L	x10exp9 /L	x10exp9 /L	x10exp9/L
MDS-2	19	F	acute myeloid leukemia	1	unrelated	2014	3-11-14	35	yes	196	6.45	3.6	1.3	1.4	y	y	no	death
MDS-4	30	M	Multiple myeloma	2	auto/unrelated	13	2-15-13	34	yes	71	6.35	4.5	1	0.7	n	y	no	alive
MDS-6	64	F	myeloproliferative syndrome	1	unrelated	359	1-22-13	37	yes	14	1.59	1	0.4	0.1	y	y	no	death
MDS-11	45	F	atypical myeloproliferative syndrome	1	unrelated	144	1-30-14	20	yes	21	1.49	0.3	1.1	0.1	y	y	no	death
MDS-12	65	F	relaspe of multiple myeloma	1	unrelated	430	2-20-14	27	yes	73	1.26	0.7	0.4	0.1	y	y	Influenza A	death
MDS-14	33	F	relapse of lymphoblastic acute leukemia	2	unrelated	335	1-8-13	25	yes	102	6.92	4.8	1.3	0.6	y	y	Influenza A + Bocavirus	death
MDS-16	33	M	relaspse of mycosis fungoides	1	related	1052	1-16-15	27	yes	256	14.05	11.3	2.1	0.7	y	y	no	death
MDS 15	40	F	acute myeloid leukemia	1	unrelated	103	2-11-14	24	yes	125	2.04	0.5	1.1	0.5	y	n	no
Means						494.4				95.3	4.5	3.0	1.0	0.5

Fig. 1 presents the timeline of events for patients involved in HCoV−OC43 clusters of infections, that why patients MDS15 and MDS16 are not mentioned. Patients were frequently admitted in hospital for consultations in the outpatient clinic or inpatient hospitalizations. They shared the same hospital units and the same healthcare worker teams. At the time of nasal swab collection, acute respiratory symptoms were present for all the patients. Other respiratory pathogens were detected in nasal swabs for only 2 patients (Table 1). No concomitant co-infection with bacteria or fungi was detected. All episodes occurred during winter.

Fig. 1

Timeline of events. NS+: Nasal swab positive for HCoV−OC43.

Characteristics of whole genome sequences

To characterize the overall diversity of HCoV−OC43 circulating among the 8 HSCT patients, 7 full-length genome sequences were obtained from nasal swabs of HSCT patients using NGS method. The sequence obtained from MDS15 patient was incomplete since a part between positions 7473 and 7689 in orf1ab was missing. Considering the small fragment of sequence missing, the partial full-length sequence obtained from this patient was also included in all phylogenetic analyses. To genotype the HCoV−OC43 isolates, phylogenetic trees were reconstructed by the maximum-likelihood method using sequences of full-length genomes and S genes obtained and compared to those retrieved from GenBank (Fig. 2 A/B). The tree obtained from full-length genomes showed that MDS sequences are divided into 4 known genotypes B, F and G (Fig. 2A) [25,26,28]. MDS4, MDS6, MDS14 and MDS12 sequences, belonging to the genotype F, are grouped into the same cluster separated from other sequences of genotype F with high bootstrap value (79.6 %). These sequences were obtained from samples collected in 2013 except for the MDS12 sequence that was obtained from sample collected in 2014. This cluster was named “cluster 2013” (Fig. 2A). Among the cluster 2013, MDS6 and MDS14 sequences are closely related (69.2 % bootstrap value), MDS12 sequence is separated from the MDS 6 and MDS14 sequences (100 % bootstrap value) and MDS4 sequence is separated from the 3 other sequences (98.5 % bootstrap value). On the tree, MDS2 and MDS11 sequences are divided within the sequences of genotype G with too low bootstrap values to separate firmly the sequences of this cluster (Fig. 2A). MDS2 and MDS11 sequences were collected during the year 2014 and will be named cluster 2014. Overall, the tree obtained from S gene sequences shows the same distribution of MDS sequences. Notably, the cluster 2013 appears closely related with previous published sequences collected from French patients in 2012 and 2013. In sum, the HCoV−OC43 sequences obtained from clinical samples in this study belong to genotypes B, E, F and G, and 6 of them appear clustered into two distinct clusters.

Fig. 2

Maximum-likelihood trees of HCoV−OC43 strains. A, B, C, D, E, F, G genotypes/clades were represented. MDS sequences are in bold. A. Tree obtained from whole-genome sequences. B. Tree obtained from S gene sequences.

The sequences of the cluster 2013 are very conserved at the nucleotide level of the full-length genomes (up to 99 %). As observed on the tree, MDS 6 and 14 sequences are the closest (with 99.93 % of identity), MDS12 sequence is less similar to the MDS6 and MDS14 sequences (99.11 % of identity) and the MDS4 sequence is the most divergent from the other of the cluster 2013 (99 % of identity). However, MDS4 sequence presents a high level of homology with MDS6 and MD14 sequences (99.81 %). Likewise, MDS2 and MDS11 sequences are very similar with 99.86 % of identities. MDS15 presents a range of homologies comprise between 95.4 % and 96.2 % and MDS16 between 98.4 and 99.2 % with the sequences of the clusters 2013 and 2014. These data reinforce the results of the phylogenetic trees and suggest that sequences of each cluster share the same origin.

To reinforce the hypothesis of common origin of clinical isolates, the divergence times of MDS sequences were estimated using BEAST method. The estimated mean evolutionary rate was 4.0 × 10−4 (3.4 × 10−4 – 4.5 × 10−4) substitutions/site/year for the full-length genome and, for the S gene, 5.5 × 10−4 (4.5 × 10-4 – 6.5 × 10−4) substitutions/site/year. Based on the full-length genome and the S gene data, the time of emergence for genotype A was estimated in the 1960s and in the late 1990s to early 2000s for genotypes B to D (Table 2 ). Similarly, emergence of genotypes D to G was estimated in the late 2000s to early 2010s. Based on the evolutionary estimates of the full-length genome and the S gene, the common ancestor of the sequences of the cluster 2013 was dated back to 2010.9 (2009.6–2012.2) and to 2010.9 (2009.8–2011.7), respectively. Furthermore, the common ancestor of the sequences of the cluster 2014 was dated back to 2011.5 (2010.7–2012.5) and to 2011.0 (2010.1–2011.7) for the full-length genome and the S gene, respectively. These estimates obtained from two different databases were consistent and show that the strains of the clusters 2013 emerged from a common ancestor in 2010.9 and those of 2014 emerged from another in 2011 +/-0.5.

Table 2

Time of the most recent common ancestors (tMRCA) with 95 % highest posterior density (95 % HPD) for HCoV−OC43 genotypes A to G and clusters based on the spike (S) gene and full-length genome.

	tMRCA (95 % HPD)
Genotype	S gene	Full-length genome
A	1963.1 (1956.6–1966.7)	1967.0 (1966.9–1967.1)
B	1999.3 (1995.1–2002.3)	2001.6 (2000.3–2002.7)
C	1993.2 (1988.7–1996.5)	2002.9 (2001.8–2003.7)
D	1998.1 (1994.8–2001.2)	2003.0 (2002.0–2003.7)
E	2006.5 (2005.7–2007.5)	2007.2 (2005.9–2008.4)
F	2007.7 (2007.0–2008.3)	2009.4 (2008.0–2010.5)
G	2010.8 (2010.0–2011.4)	2010.9 (2010.3–2011.5)
Cluster 2013	2010.9 (2009.8–2011.7)	2010.9 (2009.6–2012.2)
Cluster 2014	2011.0 (2010.1–2011.7)	2011.5 (2010.7–2012.5)
Mean evolutionary rate (nt subst/site/year)	5.5 × 10−4 (4.5 × 10−4 – 6.5 × 10−4)	4.0 × 10−4 (3.4 × 10−4 – 4.5 × 10−4)

Using the prototype ATCC VR759 as the reference strain for amino acid positions, sequence MDS16 shares 19 amino acid substitutions with other sequences of genotype B and presents 18 specific amino acid substitutions mapped across the whole genome relative to the other sequences of the genotype B (Fig. 3 A). Likewise, sequences of the cluster 2013 show the specific substitution of the genotype F (Y177 H) and 9 additional amino acids substitutions relative to the others of the genotype F [25]. Notably, MDS6, MDS12 and MDS14 sequences show 12 additional substitutions relative to the genotype F, which are absent on the MDS4 sequence. MDS15 sequence, for which only the S region could be analyzed, shares the same amino acid substitutions with the other sequences of the genotype E and shows only 1 specific substitution relative to the other sequences of the same genotype. MDS2 and MDS11 sequences appear similar to the others of genotype G.

Fig. 3

Signature amino acid substitution differences. A. Amino acid substitution differences across the full-length genome. B. Amino acid substitution differences among the S gene sequences.

Characteristics of HCoV−OC43 strains

Since S protein sequences of other French HCoV−OC43 isolates have been previously published, phylogenetic analyses show that MDS sequences of the cluster 2013 are clustered to other French sequences on the phylogenetic tree, whereas sequences of the cluster 2014 appear more related to Malaysian, Chinese and American sequences (Fig. 2B) [29]. When signature amino acid substitutions were analyzed, for the genotype B, 3 amino acid substitutions appear on the sequence MDS16: absence of substitution D471 N, substitutions L482 F that is common to other sequences of genotype B, C, F and G and L505S that is also observed on sequences of genotype C. For the genotype E, the amino acid substitution N767I is observed on the sequence MDS15 and had never been described previously. For the genotype F, the amino acid substitution P38S is observed on the MDS4, MDS6, MDS12 and MDS14 sequences which are shared with genotype B and G. The substitutions D1264H and C1319S are absent from the MDS4, MDS6, MDS12 and MDS14 sequences as from the other French sequences K963240, K963243 and K963244. Finally, 2 amino acid substitutions (A1054S and D1252 G) are observed on MDS6, MDS12 and MDS14 sequences and also on previously published sequences K963243 and K963244 (Fig. 3B).

These data show that any amino acid substitution appear specific to sequences of genotype B and G circulating in France. But, sequences circulating among French people belonging to the genotype F had a peculiar molecular signature (absence of substitution D1264H and C1319S, substitution A1054S and D1252 G).

Discussion

In this work, we describe retrospectively 8 HCoV−OC43 infections having occurred during the winters 2013 and 2014 at the HSCT Unit of the Lille University hospital. Six patients were divided into 2 clusters, respectively named 2013 and 2014 on the basis of the year of occurrence and genomic characteristics.

HCoV−OC43 cases described in this report could be classified as healthcare-associated infections defined as a group of infections emerging among patients who come from a community with a history of previous exposure to healthcare, but do not fit the nosocomial infection criteria [30]. The 6 patients divided into the two clusters of cases regularly attended the hospital for in- and outpatient cares, as consultations, day and inpatient hospitalizations. So, they were frequently in contact with health care professionals and hospital visitors that may be considered as multiple sources of HCoV infections, during periods of usual circulation of HCoV−OC43 in France [31]. Outside the hospital, they resided in geographically non-related places and no common activity could be identified. So, the hospital admission seems the only link between the patients, even though we were unable to identify any direct contact between the patients. HCoV−OC43 can be transmitted either by an asymptomatic patient with prolonged shedding, by healthcare worker asymptomatic or with an upper respiratory tract infection [10,32]. Phylogenetic analyses showed that whole-genome sequences included in a same cluster were very closely related, had very high levels of similarities and shared same common ancestors. All these elements provide additional evidence that these infections could be classified as healthcare-associated infections. Similarly, to what is observed for influenza viruses, parainfluenza virus type 3 and respiratory syncytial virus, our data show that HCoV−OC43 is involved in healthcare-associated infections [[33], [34], [35]].

Because we were unable to identify the index case of HCoV−OC43 infections of the cluster 2013, the MDS12 patient infection during the winter 2014 with a strain highly related to those circulating the previous year is surprising. However, previous studies have shown that genomes of other high and low pathogenic coronaviruses presented low variation rates during hospital-acquired outbreaks [36,37]. Our data show that sequences of HCoV−OC43 appear also very stable.

Phylogenetic trees based on the whole genome and S gene HCoV−OC43 show that genotypes B, G, E and F currently circulate among French patients, in agreement with other report [29]. Even if HCoV−OC43 was previously detected among patients with haematology malignancy, the involvement of HCoV−OC43 genotypes in this type of infection was rarely investigated [[38], [39], [40]]. Indeed, the analysis of more HCoV−OC43 strains from other healthcare-associated infections will reveal the relative prevalence of each genotype in this type of infections in different localities.

The divergence time of sequences included in the phylogenetic analyses was calculated to accumulate to collect additional information showing the common origin of clinical isolates. Our results are similar to those previously described, validating the tMRCA method used [25,28]. Indeed, our results show that each cluster of sequences has diverged from a distinct ancestor in 2010.9 for the cluster 2013 and in 2011.5 for the cluster 2014.

Due to limited HCoV−OC43 full-length genome published in databases, phylogenetic analyses were conducted on S gene sequences. All results obtained from whole genome and S gene sequences were similar. However, phylogenetic analyses from the numerous published S gene sequences allowed comparison of HCoV−OC43 sequences isolated in Lille with those isolated in different parts of the world. Strains described here represent probably a part of strains currently circulating among European people. However, our data show that both kinds of HCoV−OC43 strains circulate among patients included in the study: those belonging to the genotype G closely related to Asian strains and other belonging to the genotype F more related to other European strains previously described [25,26,28,29,41]. Amino acid substitution analyses showed that some substitutions especially in the S protein could be specific of strains circulating in European population but the number of strains analyzed here and the number of sequences obtained from European isolates of HCoV−OC43 are too small and further studies will be needed to confirm our data.

We are aware that our study presents some limitations. As previously mentioned, we were unable to identify the index case of the HCoV−OC43 infections in the HSCT unit and determine precisely the contact history of the patients. Thus, we cannot exclude that other patients were involved in the HCoV−OC43 healthcare-associated infections described. In the same way, we were unable to detect infection or asymptomatic portage in healthcare workers of the HSCT unit. Some reports have shown that different HCoV−OC43 strains circulated during a same winter, suggesting that there may be multiple channels of HCoV introduction in hospital [25,42,43]. In fact, the HCoV strains of healthcare-associated infections were probably contracted from one of the strains circulating in the community at that time. For all these reasons, cases reported in this study did not correspond to the definition of nosocomial transmission of pathogens and we used the term of coronavirus healthcare-associated infections.

Another weakness of our work concerns the determination of pattern of amino acid substitutions specific of strains circulating among French people. Indeed, our data suggest that some substitutions might be specific to strains circulating in Europe but the analysis of a larger number of strains would be necessary to confirm our hypothesis.

To conclude, we retrospectively described the molecular investigation of HCoV−OC43 infections in a HSCT unit using whole-genome sequencing combined with advanced phylogenetic analyses. These tools allowed us to associate HCoV−OC43 infections with healthcare. Moreover, we suggest that two kinds of HCoV−OC43 strains circulate among the French population, one sharing common ancestors with Asian strains and the other closely related to European lineage.

Funding sources

This work was supported by a Lille Hospital University grant (EPI-CoV) and by and Lille University (PEPS 2015).

Transparency document

CRediT authorship contribution statement

Delphine Beury: Investigations, Formal analyses. Léa Fléchon: Investigations, Formal analyses. Florence Maurier: Investigations, Formal analyses, Data curation. Ségolène Caboche: Conceptualization, Investigations, Data curation, Writing - original draft. Jean-Stéphane Varré: Conceptualization, Data curation, Writing - original draft. Hélène Touzet: Conceptualization, Data curation, Writing - original draft. Karine Faure: Conceptualization, Resources, Writing - original draft. Jean Dubuisson: Conceptualization, Writing - review & editing, Funding acquisition. David Hot: Conceptualization, Data curation, Writing - original draft. Benoit Guery: Conceptualization, Writing - review & editing, Funding acquisition. Anne Goffard: Conceptualization, Resources, Writing - original draft, Funding acquisition, Project administration.

Declaration of Competing Interest

The manuscript had not been published or presented elsewhere, and the authors have no conflict of interest to declare.