Molecular determinants and heterogeneity underlying host response to EV-A71 infection at single-cell resolution

ABSTRACT

The pathogenic human enterovirus EV-A71 has raised serious public health concerns. A hallmark of EV-A71 infection is the distortion of host transcriptomes in favour of viral replication. While high-throughput approaches have been exploited to dissect these gene dysregulations, they do not fully capture molecular perturbations at the single-cell level and in a physiologically relevant context. In this study, we applied a single-cell RNA sequencing approach on infected differentiated enterocyte cells (C2BBe1), which model the gastrointestinal epithelium targeted initially by EV-A71. Our single-cell analysis of EV-A71-infected culture provided several lines of illuminating observations: 1) This systems approach demonstrated extensive cell-to-cell variation in a single culture upon viral infection and delineated transcriptomic differences between the EV-A71-infected and bystander cells. 2) By analysing expression profiles of known EV-A71 receptors and entry facilitation factors, we found that ANXA2 was closely correlated in expression with the viral RNA in the infected population, supporting its role in EV-A71 entry in the enteric cells. 3) We further catalogued dysregulated lncRNAs elicited by EV-A71 infection and demonstrated the functional implication of lncRNA CYTOR in promoting EV-A71 replication. Viewed together, our single-cell transcriptomic analysis illustrated at the single-cell resolution the heterogeneity of host susceptibility to EV-A71 and revealed the involvement of lncRNAs in host antiviral response.

KEYWORDS:

• Enterovirus
• single-cell sequencing
• host-pathogen interactions
• lncRNAs
• ev-A71
• enterocytes

INTRODUCTION

Several types of human enteroviruses are pathogenic to humans and have emerged as serious public health threats. Among the viruses, enterovirus A71 (EV-A71) reportedly causes severe neurological complications in infants and young children, despite a faecal-oral transmission route [1,2]. EV-A71 intrusion into host cells depends on specific viral receptors on the cell membrane, such as SCARB2 and PSGL-1 [3,4], and could also be facilitated by other cellular proteins, such as ANXA2 and hWARS [5,6]. Upon binding to the viral receptors, the virus enters host cells via endocytosis and subsequently releases its positive-stranded RNA genome to the cytoplasm [7]. The viral genomic RNA serves as the mRNA template for translating a precursor polyprotein, which undergoes further proteolytic processing by virus-encoded proteases to generate functional viral proteins in infected cells.

A hallmark of EV-A71 infection is the distortion of host antiviral responses in favour of viral replication in the infected cells [8]. In this context, EV-A71 is known to dysregulate host gene expression through [9] interfering with host mRNA processing, such as splicing and polyadenylation, and protein translation [10–12]. These molecular perturbations disrupt host cell biology and ultimately disarm the antiviral response. To unravel the host–virus interaction underlying the EV-A71 infection cycle, several high-throughput approaches, such as cDNA microarray, proteomic analysis, and RNA sequencing, have been exploited to dissect the dysregulated gene expression in the infected cells [13–16]. More recently, single-cell transcriptomic technology has emerged as a powerful tool to investigate infection-associated signalling and gene alterations at the single-cell resolution. This systems approach has been used to demonstrate cell-to-cell variation in a single culture upon viral infection [17], and to outline transcriptomic differences between the influenza A virus-infected and bystander cells, further illuminating host innate immunity and paracrine signalling during the infection [18]. Hence, the single-cell transcriptomics represents a feasible approach to delineate host antiviral immune response elicited by EV-A71 infection.

Aside from protein-coding genes, the non-coding RNAs, such as miRNAs and long noncoding RNAs (lncRNAs), are known to constitute an additional layer of regulation in the cellular response to virus infection. Even without coding potential, lncRNAs have been further linked to gene expression regulation and innate immune response by acting via multiple mechanisms [19,20]. Interestingly, their potential involvement in EV-A71 infection is revealed by the differential expression of lncRNAs in response to EV-A71 infection [21,22]. However, functional outcomes of these dysregulated lncRNAs in association with viral infection remain presently unclear.

Since most of the aforementioned profiling studies were performed in cell lines, the profiles of host gene alterations in a physiologically relevant context, such as the initially infected anatomical regions within the intestinal tracts, require further investigation. Drummond et al. previously demonstrated that transcriptomes of differentiated enteric epithelium are distinct from those of the undifferentiated cell line [13], supporting the notion that differentiated intestinal epithelial cells and enteroids are better models for recapitulating viral replication and host response during EV-A71 infection [13,23,24]. For that purpose, we have established an infection model by using human enteric epithelial cells differentiated from an intestinal cell line C2BBe1. We then set out to characterize at the single cell level dysregulation of host gene expression upon EV-A71 infection.

On the basis of this experimental paradigm, we then made the following critical discoveries: In search of potential viral receptors in the host cells, we found that ANXA2, which is known to bind EV-A71 capsid protein VP1, is coordinately expressed with EV-A71 viral RNA (vRNA). In addition, compartmentalized expression analysis of innate immunity genes in relation to vRNA levels revealed the bystander effect in host cells during EV-A71 infection. Finally, by this single-cell transcriptomics approach, we identified a novel EV-A71 infection-responsive lncRNA, CYTOR, which might act as a positive regulator of viral replication.

MATERIALS AND METHODS

Cells and enteric cell differentiation

Enteric epithelial cell line C2BBe1 (ATCC CRL-2102) was cultivated in Dulbecco’s Modified Eagle Medium (DMEM) containing 10% foetal bovine serum (FBS), 0.25 M L-glucose, 10 mg/L holo-Transferrin human (Sigma-Aldrich, Cat. No. T0665), 1 mM sodium pyruvate (Gibco), and 4 mM L-glutamate (Gibco) at 37°C and 5% CO₂. To differentiate the C2BBe1 cells as described previously [24], the cells were cultivated with a medium mixed with the aforementioned C2BBe1 culture medium and Enterocyte Differentiation Medium (Corning) for 24 h, and then maintained in the Enterocyte Differentiation Medium containing 1% Insulin-Transferrin-Selenium-A (Gibco) for 48 h. RD cells (ATCC CCL136) were cultivated in DMEM containing 10% FBS and 1% non-essential amino acids at 37°C and 5% CO₂.

Virus and viral titration and infection

EV-A71 isolate of TW-00073-2012 (C4 subgenotype; Genbank accession number MG756724) was amplified in RD cells and titrated by plaque formation assay on monolayers of RD cells. For infection, differentiated C2BBe1 cells were washed with PBS before being incubated with the EV-A71 isolate at an MOI of 20 for 1 h in DMEM with 2% FBS and subsequently maintained in fresh C2BBe1 culture medium with 2% FBS. At the indicated time points post-infection, the infected cells were subjected to RNA extraction, cell lysate collection, or single-cell RNA sequencing.

Single-cell RNA sequencing and data analysis

Infected or mock-infected differentiated C2BBe1 cells were detached from the cell culture plates using trypsin-0.05% EDTA (Gibco), washed and then passed through a cell strainer to obtain single-cell suspensions at 1,000 cells/µl. Single-cell encapsulation was performed with beads using the 10× Genomics Single Cell v2 kit, based on the manufacturer instructions (10× Genomics, Pleasanton, CA). Single-cell gel beads in the emulsion were generated, in which reverse transcription was performed. This was followed by cDNA amplification, library construction, and multiplexing index tagging. After library quality control assessment by Bioanalyzer, sequencing was performed on an Illumina HiSeq 2500 (Illumina, Inc., San Diego, CA).

For data preprocessing and analysis, paired-end reads from Illumina sequencing were first aligned to a combination of the human hg38 and EV-A71 polyprotein gene using the Cell Ranger 3.0 (10× Genomics, Pleasanton, CA) pipeline. The resulting unique molecular identifier (UMI) count matrices were automatically filtered to remove background and non-cellular barcodes, and low UMI counting cells were excluded in the advanced analysis. Three samples representing mock-infection and infection at 12 h and 24 h were analysed, from which a total number of 2,201 cells, with a median of 4,219 genes per cell, were obtained after filtering and subjected to subsequent evaluation. All processed sequencing data were then imported into the Cell Loupe browser v4.1.0 (10× Genomics, Pleasanton, CA) for further analyses, including t-SNE plot, visualization, cell clustering, and differential gene expression profiling.

Reverse transcription and quantitative PCR

Total RNA of the mock or EV-A71-infected cells was extracted by TRIzol reagent (Invitrogen) following the manufacturer’s protocol and subjected to reverse transcription (RT) by SuperScript III reverse transcriptase (Invitrogen) with oligo-dT primer. The RT products were then subjected to quantitative PCR (RT-qPCR) reactions using KAPA SYBR FAST qPCR kit (Roche) with specific primers. The primer sequences used in this research are listed in Supplemental Table S1. The relative gene expression was calculated by the ΔΔCq method.

Immunoblotting and antibodies

Cellular protein extracts were separated by SDS-PAGE and transferred to PVDF membranes. The membranes were blocked with 5% skim milk and reacted individually with anti-EV-A71 3D (GeneTex) and anti-β-actin primary (Sigma-Aldrich) antibodies. After incubating with horseradish peroxidase-conjugated secondary antibody, the specific proteins were detected by a chemiluminescent substrate.

Determination of EV-A71 replication in differentiated C2BBe1 cells with siRNA and plasmid transfection

To specifically knockdown target gene (CYTOR), we used synthetic siRNAs with the following sequences: siCYTOR#1: 5ʹ-CACACUUGAUCGAAUAUGA-3ʹ and 5ʹ-UCAUAUUCGAUCAAGUGUG-3ʹ; siCYTOR#2: 5ʹ – GUCGUGUGUACAUCAUUGG and 5ʹ – CCAAUGAUGUACACACGAC-3ʹ; control siRNA: 5ʹ-GAUCAUACGUGCGAUCAGA-3ʹ and 5ʹ-UCUGAUCGCACGUAUGAUC-3ʹ. For constructing the CYTOR expression vector, the exonic region of CYTOR was amplified from the cDNA sample. The resulting PCR products were ligated into the cloning vector using the HE Swift Cloning Kit (TOOLS Life Science, Taiwan), and subsequently sub-cloned into the expression vector pcDNA3.1 (−). For siRNA and plasmid delivery, siRNAs or CYTOR overexpression construct were first mixed with Lipofectamine 3000 (Invitrogen) in Opti-MEM (Gibco) and incubated at room temperature for 25 min. The mixture was then added to C2BBe1 cell suspension in a medium mixed with Opti-MEM, C2BBe1 culture medium and Enterocyte Differentiation Medium. The cells were seeded to culture plates and incubated at 37°C and 5% CO₂ for 24 h. The cells were then maintained in Enterocyte Differentiation Medium with 1% Insulin-Transferrin-Selenium-A for 48 h. The siRNA-treated or CYTOR-overexpressing differentiated C2BBe1 cells were infected with EV-A71 TW-00073-2012 isolate at MOI of 0.1 or 0.001 in C2BBe1 culture medium with 2% FBS. At the indicated time points post-infection, total cellular RNA, conditioned media and lysates of the infected cells were collected for RT-qPCR and virus titration.

RESULTS

Establishment of differentiated enteric epithelial cells as an EV-A71 susceptible model

To investigate at the single cell level how the host responds to EV-A71 infection in the digestive tract, which is the initial site of enterovirus infection, we first set out to model the natural infection course in cell culture. To this end, we selected human enteric cell line C2BBe1, which can be differentiated into enteric epithelial cells as previously described [24]. C2BBe1 differentiation was characterized by the induction of marker genes, such as Villin (VIL1), E-cadherin (CDH1), and CDX2, after culturing in differentiating medium for 3 days (Fig. 1A and Supplemental Fig. S1A). The differentiated enteric epithelial cells were then infected with EV-A71 at a multiplicity of infection (MOI) of 20 for 12 and 24 h and analysed for the viral mRNA and EV-A71 3D protein expression, respectively, by RT-qPCR and Western blot methods (Fig. 1B and C). Further, the synthesis of EV-A71 protein was monitored by immunofluorescence and flow cytometry methods as a measure of the infection efficiency (Supplemental Fig. S1B and S1C). The robust detection of viral gene products demonstrated that the differentiated C2BBe1 cells could be readily infected by EV-A71. These results established the differentiated C2BBe1 cells as a feasible enteric cell model for subsequent single-cell transcriptomic analysis of host–virus interaction.

Figure 1. Infection of differentiated C2BBe1 cells by EV-A71. (A) C2BBe1 cells were cultivated in the differentiation medium for 1, 2, or 3 days. The mRNA expression of differentiation markers of enteric epithelial cells, Villin (VIL-1), E-cadherin (CDH1) and CDX2, were detected in differentiated C2BBe1 cells by RT-qPCR. (B) The differentiated C2BBe1 cells were infected with EV-A71 at MOI of 20 for 12h and 24 h. Lysates of the infected cells were collected and analysed by immunoblotting with anti-EV-A71 3D and anti-β-actin antibodies. (C) Total RNA of the differentiated C2BBe1 cells infected with EV-A71 were subjected to RT-qPCR for viral RNA detection. p.i., postinfection. *: p < 0.05; **: p < 0.01

Dynamic transcriptomic changes in EV-A71-infected culture at single-cell resolution

To obtain a global profile of EV-A71 infection and EV-A71-host interaction at the single-cell resolution, we performed single-cell RNA expression profiling (10× Genomics) of differentiated C2BBe1 cells infected with EV-A71 at MOI of 20. Three types of culture were collected for this analysis – mock, infection for 12 h, and infection for 24 h – from which we generated transcriptome data of 585, 646, and 970 cells, respectively. For a direct and comprehensive cross-treatment comparison, data from these three sequencing data were combined for subsequent analyses. t-Distributed stochastic neighbour embedding (tSNE) visualization of the single-cell data, followed by a clustering analysis (see Materials and Methods), was used to examine the global expression profile of the single cells. As shown by the tSNE plots, there was a clear distinction in the distribution of mock-infected vs. infected cells (12 and 24 h combined) in terms of transcriptome profiles (Fig. 2A upper panel and Supplemental Fig. S2A). Among the cells in the two infection time points (n = 1,616), a total of 1,371 cells were positive for the EV-A71 vRNA, indicating a considerable extent of EV-A71 infectivity in this differentiated enteric cell line (Fig. 2B upper panel and Supplemental Fig. S2B). However, infected cells, regardless of infection time (Fig. 2A; Supplemental Fig. S2A) or expression levels of the vRNA (Fig. 2B; Supplemental Fig. S2B), seemed to distribute in a broad cluster, indicating that EV-A71 infection may lead to generally concerted transcriptome perturbations.

Figure 2. Single-cell transcriptome profiling of EV71-infected C2BBe1 cells. (A) & (B) Based on gene expression profiles, cells were clustered and analysed by t-SNE method. Cells were classified into distinct subgroups, as indicated by colours, according to infection condition/time (A) or EV71 RNA expression (B). Numbers in the parentheses correspond to number of cells in each category. Heatmap representation of the genes differentially expressed between the indicated comparisons. Extent of the differential expression was based on the colour scale on the bottom

Next, we identified the genes that showed the greatest changes in mean expression (normalized to a unique molecular identifier, or UMI, counts) between different selected sub-population of the culture, as shown by the heatmaps in the lower panel of Fig. 2A and 2B, respectively. Distinct sets of the highly differentially expressed genes between the comparisons provided substantial evidence for enterovirus-mediated transcriptomic regulation. As independent validation, we also confirmed the infection-responsive differential expression of selected genes by RT-qPCR assay (Supplemental Fig. S3). Most notably, biological process analysis of differentially expressed genes in the mock vs. EV-A71-infected cells (including vRNA-negative and positive cells) showed the enrichment in genes related to metabolic, cell apoptotic and proliferation processes (Table 1). The pathway analysis via KEGG annotation further revealed a preponderance of genes involved in cancer pathways and metabolic processes, including small cell lung cancer signalling and glutathione metabolism (Supplemental Table S2). These unique gene annotations indicated that viral infection might be functionally correlated with metabolic and proliferation dysregulation in the host cells.

Table 1. Biological process analysis for the dysregulated genes among mock, EV-A71-vRNA-negative, and vRNA-positive cells

GO ID	Term	Count	P-value	Genes
42448	Progesterone metabolic process	4	9.23 x 10^−6	AKR1C3, AFP, AKR1C2, AKR1C1
43066	Negative regulation of apoptotic process	11	7.23 x 10^−5	IER3, ID1, YBX3, TGM2, RGN, BCL3, BCL2L1, BIRC3, MYC, CITED2, SOD2
48661	Positive regulation of smooth muscle cell proliferation	5	2.09 x 10^−4	ID2, SKP2, TGM2, ITGA2, MYC
8283	Cell proliferation	9	4.13 x 10^−4	DKC1, SKP2, ITGA2, BCL2L1, PDZK1, MYC, IGFBP4, CITED2, GNL3
6749	Glutathione metabolic process	4	2.58 x 10^−3	GSTA1, GSTA2, GSTA4, SOD2
30855	Epithelial cell differentiation	4	4.86 x 10^−3	GSTA1, GSTA2, AKR1C2, AKR1C1
8284	Positive regulation of cell proliferation	8	7.70 x 10^−3	AKR1C3, ODC1, AKR1C2, CEACAM6, LAMC2, PBX1, BCL2L1, MYC
122	Negative regulation of transcription from RNA polymerase II promoter	10	8.36 x 10^−3	CEBPA, CCND1, DKK1, ID2, ID1, YBX3, NR2F2, MYC, S100A1, CITED2
8152	Metabolic process	5	9.23 x 10^−3	GSTA1, GSTA2, ACSS1, GSTA4, ACAA1
32355	Response to estradiol	4	1.00 x 10^−2	CCND1, MMP15, NR2F2, MYC
8584	Male gonad development	4	1.09 x 10^−2	AKR1C3, YBX3, BCL2L1, CITED2
42493	Response to drug	6	1.65 x 10^−2	CCND1, ITGA2, GAL, MYC, DUSP6, SOD2
6364	rRNA processing	5	2.07 x 10^−2	CCDC86, NOLC1, DKC1, DDX21, WDR43
55114	Oxidation-reduction process	8	2.55 x 10^−2	TM7SF2, AKR1C3, GPX2, AKR1C2, HSD17B14, GLUD1, AKR1C1, SOD2
45892	Negative regulation of transcription, DNA-templated	7	3.51 x 10^−2	CEBPA, ID2, ID1, CEBPD, BCL3, NR2F2, CITED2
10628	Positive regulation of gene expression	5	3.94 x 10^−2	ID2, LDLR, ID1, MYC, CITED2

Despite the broad clustering of the infected cells, information relevant to this study could still be extracted by collapsing these clusters according to the expression of viral genes as well as the infection time. For this purpose, we re-assigned cell groups based on the infection time-points and viral RNA detection (upper panel of Fig. 3A and Supplementary Fig. S2C), and clusters-based differentially expressed genes were then analysed and represented by a heatmap with correspondence to infection states (lower panel of Fig. 3A). A subset of genes differentially expressed between clusters was also confirmed by independent RT-qPCR assay (Supplemental Fig. S4). Interestingly, we first noticed that a considerable proportion of cells infected for 24 h expressed detectable vRNA as compared to those infected for 12 h (Fig. 3B), in line with the notion that prolonged infection is correlated with an increased level of viral products. In addition, while the majority of the differentially expressed genes exhibited consistent expression alterations in response to a viral infection – either downregulated or upregulated at both time points of infection, a particular set of genes was specifically activated at 24 h post-infection (e.g., LCN2 and LTB). Interestingly, the distinct gene expression patterns between EV-A71-expressing vs. EV-A71-negative cells demonstrated transcriptomic dysregulation elicited specifically by the enteroviral infection, and further raised the possibility of the bystander effect of host viral response. These single-cell analyses provided better resolution for interpreting gene regulation in host cells conferred by enteroviral infection.

Figure 3. Sub-group expression profiles of EV-A71 vRNA. (A) Further sub-grouping of the cells was done based on both the presence of EV71 RNA expression (vRNA) as well as the infection time points (12 vs. 24 hrs), and labelled accordingly in the t-SNE plot. Heatmap representation depicts the genes differentially expressed between the indicated comparisons. Colours correspond the extent of the differential expression, based on the colour scale on the bottom. (B) The numbers of cells positively detected with EV-A71 vRNA in the infection cultures

Expression correlation between EV-A71 RNA and viral receptors in differentiated enterocytes

We then examined whether the expression of known EV-A71 receptor and facilitating factor genes, including SCARB2, PSGL-1 (SELPLG), ANXA2, and WARS, was coordinated with EV-A71 RNA expression. As shown in the left panel of Fig. 4A, two of the genes, SCARB2 and ANXA2, were generally expressed in the differentiated C2BBe1 cells, whereas the expression of WARS and SELPLG was low or hardly detected. Incidentally, ANXA2 expression was significantly higher than the other receptor genes (Fig. 4A, right panel). Further expression analysis revealed that ANXA2 and SCARB2 were downregulated by enteroviral infection (Fig. 4B), while WARS expression remained unchanged during the course of infection. Although ANXA2 expression levels were similar between vRNA-positive cells and those devoid of vRNA (Fig. 4C), the cells highly expressing ANXA2, but not SCARB2, were found to harbour abundant EV-A71 vRNA (Fig. 4D; Supplemental Fig. S5). These correlations in the expression are a strong indication that ANXA2 is a critical host determinant underlying EV-A71 infection to differentiated enteric cells.

Figure 4. Expression correlation between EV-A71 RNA and viral receptors in differentiated enterocytes. (A) Numbers of cells positive for SCARB2, SELPLG, ANXA2, and WARS expression as detected by the single-cell sequencing are shown in the bar graph on the left, and the corresponding read counts of respective genes are depicted in right panel. (B) Average sequencing read counts of SCARB2, SELPLG, and WARS genes in the Mock and infected clusters were determined and presented, and statistical significance was assessed by student’s t-test. (C) Cells were clustered based on detectable EV-A71 RNA expression, and average read counts of SCARB2, SELPLG, and WARS genes in these clusters are shown. (D) Pearson’s correlation analysis for coordinated expression between ANXA2 and EV-A71 RNA in the infected culture was performed, with the regression line (red dotted line) and statistical significance (p-value) of the analysis also depicted

EV-A71-elicited host responses of differentiated enterocytes profiled at single-cell level

To assess how host immune signalling changes in response to viral infection in each cell, we next assessed expression alterations of antiviral response and cytokine/chemokine genes at the single-cell level. Towards this end, our gene expression profiling was based on the following comparisons: infection (combined two time points) vs. mock, vRNA-negative vs. vRNA-positive cells in the infection groups, and vRNA-negative infected cells vs. control mock cells. As shown by the heatmap in Fig. 5A, which summarizes the fold changes in the relative abundance of selected immune response effectors, there was a general upregulation of these genes in cells subjected to infection, regardless of the presence of vRNA. This observation lends strong support to the bystander effect of viral infection, through which the neighbouring non-infected cells also induce immune response signalling. In addition, when we compared the vRNA+ vs. vRNA- cells at either 12 or 24 h postinfection, we found that several genes involved in innate immunity, such as IFIT1, IFITM1, DDX58 (RIG-I) and ISG15, were expressed to a lesser extent by the vRNA+ cells. As additional validation, we confirmed the expression of selected genes differentially altered between clusters by independent RT-qPCR assay (Supplemental Fig. S6). However, expression of these genes did not decline further in the cells with high levels of vRNA (Fig. 5B), suggesting that the progressive replication of viral RNA might not be associated with repressed immunity.

Figure 5. EV-A71-elicited host responses of differentiated enterocytes profiled at single-cell level. (A) Clustering of single cells was done on the basis of the infection time-points and EV-A71 RNA expression, and the expression profiles of immune-associated genes among the indicated clusters are shown by the heatmap representation. (B) Relative gene expression (read counts per million) of the immune-associated genes in the indicted cluster comparison is quantitatively shown by the bar graphs, and statistical significance was measured by student’s t-test

Since previous studies hypothesized that infection by EV-A71 could dysregulate host innate immune response to facilitate viral replication, we examined the expression of type I and III IFNs, which are the frontlines against virus infection. However, in accordance with the previous findings, we did not find significant induction in mRNA expression of IFN-β, λ1, and λ2 upon EV-A71 infection in the differentiated C2BBe1 cells, either in infected cell population as a whole (Supplemental Fig. S7A) or at the single-cell level (Supplemental Fig. S7B). Collectively, this down-regulation of host defence might be attributed to a virus-mediated antagonism of the cellular gene expression machinery.

Subsequently, we set out to investigate the implications of the vRNA expression-dependent gene expression alterations elicited in infected cultures. For that purpose, we compartmentalized single cells into subgroups based on the abundance of viral RNA and assigned six different vRNA-expressing clusters of approximately equal cell number, from vRNA (+1) to vRNA (+6) (Fig. 6A and Supplemental Fig. S8). Next, we analysed cellular gene networks that change coordinately with the extent of viral infection. The heatmap representation shown in Fig. 6B summarized the genes and corresponding changes in expression that are correlated with viral infectivity (i.e. EV-A71 vRNA expression). In particular, expression alteration of certain genes (e.g. AFP and FAM3B) was correlated with the abundance of vRNA, implying viral-mediated regulation of host susceptible factors for infection. Further biological process analysis of these highly variable genes revealed that genes involved in cell adhesion, apoptosis, and inflammatory response were correlated to vRNA abundance (Table 2). Interestingly, pathway (KEGG) analysis of the identified genes showed that cancer-related genes, such as those in the Hippo, NF-kB and TNF signalling pathways, might be the prominent cellular changes underlying host response to viral infection.

Figure 6. Host gene expression perturbations in association with the extent of viral RNA abundance. (A) Cells were sub-grouped into seven clusters according to the levels of EV71 RNA expression, and depicted according in the t-SNE representation. Split view for the clusters is depicted in the lower panel. (B) Heatmap representation for the top differentially expressed genes among the sub-groups. Colours correspond to folds of difference

Table 2. Biological process and KEGG pathway analyses of EV-A71 RNA clusters in Figure 6

GO ID	Term	Count	P-value	Genes
32355	Response to estradiol	4	1.52 x 10^−3	CCND1, MMP15, NR2F2, MYC
97421	Liver regeneration	3	2.38 x 10^−3	CCND1, RGN, MYC
45785	Positive regulation of cell adhesion	3	5.17 x 10^−3	SAA1, TGM2, ITGA2
43066	Negative regulation of apoptotic process	6	5.32 x 10^−3	IER3, ID1, TGM2, RGN, BIRC3, MYC
48661	Positive regulation of smooth muscle cell proliferation	3	9.86 x 10^−3	TGM2, ITGA2, MYC
30593	Neutrophil chemotaxis	3	1.18 x 10^−2	CCL20, SAA1, CXCL3
50729	Positive regulation of inflammatory response	3	1.43 x 10^−2	LDLR, TGM2, ITGA2
187	Activation of MAPK activity	3	2.94 x 10^−2	DUSP5, SAA1, DUSP6
42493	Response to drug	4	4.01 x 10^−2	CCND1, ITGA2, MYC, DUSP6
KEGG ID	Pathway	Count	P-value	Genes
hsa05222	Small cell lung cancer	4	2.7 x 10^−3	CCND1, ITGA2, BIRC3, MYC
hsa04390	Hippo signaling pathway	4	1.33 x 10^−2	AFP, CCND1, ID1, MYC
hsa04064	NF-kappa B signaling pathway	3	3.37 x 10^−2	RELB, BIRC3, LTB
hsa05200	Pathways in cancer	5	3.92 x 10^−2	CCND1, SLC2A1, ITGA2, BIRC3, MYC
hsa04668	TNF signaling pathway	3	4.92 x 10^−2	CCL20, CXCL3, BIRC3

Functional characterization of lncRNA dysregulation in association with EV-A71 infection

Despite the lack of protein coding potential, lncRNAs have emerged as critical determinants in the regulation of cellular response against virus infection. Given the comprehensive transcriptomic profiling provided by our single-cell analysis, we next sought to explore potential virus infection-associated lncRNAs. Towards this end, we have focussed on the annotated lncRNAs and analysed at the single-cell level their expression changes upon EV-A71 infection (Fig. 7A). Our interrogation first uncovered several lncRNAs with detectable abundance and infection-responsive alteration. Intriguingly, as shown by the hierarchical clustering analysis, several lncRNAs such as MALAT1 and NORAD were found to express inversely to the abundance of vRNA in the infected cells, suggesting that there might be a unique temporal regulation of lncRNA expression during viral infection. One of the differentially expressed lncRNAs, CYTOR, was identified as an upregulated gene in the infected cells in our single-cell analysis, of which enhanced expression in the infected cells was further confirmed by RT-qPCR experiment (Supplemental Fig. S9). Owing to this unique expression pattern and its hitherto uncharacterized link to enterovirus infection, we sought to investigate the role of lncRNA CYTOR during EV-A71 infection. To this end, we knocked down lncRNA CYTOR expression in differentiated C2BBe1 cells by synthetic duplex siRNA, targeting efficiency of which was confirmed by a RT-qPCR experiment (Fig. 7B, left panel). We subsequently examined virus replication and found that reducing the level of CYTOR led to a decreased production of EV-A71 vRNA at indicated postinfection (MOI of 0.1 and 0.001; Fig. 7B, right panel). In parallel, viral titration was also reduced by CYTOR knockdown at 48 h postinfection (Fig. 7C). We further performed ectopic overexpression of CYTOR in differentiated C2BBe1 cells and subsequently confirmed the upregulation of EV-A71 RNA production as well as the viral titre in CYTOR-overexpressing cells (Fig. 7D). These results collectively demonstrated that lncRNA CYTOR acts as a positive regulator of EV-A71 RNA replication in differentiated enteric cells, and further supported the significance of noncoding RNA-mediated regulation in the viral life cycle progression.

Figure 7. Knockdown of CYTOR lncRNA affects replication of EV-A71 in differentiated C2BBe1 cells. (A) The heatmap shows the top differentially expressed lncRNAs among the sub-groups as described in Figure 6A. Colours correspond to folds of difference. (B) Control and CYTOR-specific siRNAs were transfected into differentiated C2BBe1 cells and subsequently infected at the indicated viral titration (MOI). Expression profiles of CYTOR (left) and viral RNA (middle) in the infected cells were analysed by RT-qPCR assay. (C) As another readout for infection extent, cells in (B) were harvested at 48 h p.i. and analysed for viral titration by the plaque assay. (D) CYTOR-overexpressing differentiated C2BBe1 cells were infected at the indicated viral titration (MOI), and subsequently subjected to the measurement of CYTOR and viral RNA expression (left and middle), and of viral titration by plaque assay (right)

DISCUSSION

Host transcriptome landscape is considerably dysregulated in the wake of EV-A71 intrusion, which either directly targets cellular transcription by viral proteins or trigger innate antiviral responses. Several platforms have been applied previously to explore the full extent of these molecular changes, with next-generation sequencing and mass spectrum-based methodology as the most notable approaches. While these tools report comprehensive profiles of the cell culture under study, changes of gene expression at the single-cell level are insufficiently characterized, thus underestimating critical aspects of viral infection, such as the heterogeneity in cellular responses. In this study, we set out to revise earlier experimental design by applying a single-cell RNA sequencing approach on a differentiated enterocyte model. Our single-cell analysis of EV-A71-infected culture provided a global and dynamic view of the host transcriptomic response to viral infection. This information at single-cell resolution illustrated host gene alterations in response to EV-A71 infection, as well as the molecular distinctions of both infected and non-infected cells at the population level. In addition, transcriptomic events occurring in the non-infected neighbouring cells and/or EV-A71-negative infected cells, which constitute the bystander effect, could be thoroughly dissected. Based on this approach, our current study revealed not only the heterogeneity of host cell susceptibility to EV-A71 enterovirus infection but also the cell-to-cell variation in gene expression profiles. Moreover, the identification of the host transcripts co-varied with EV-A71 vRNA may illuminate the molecular distinctions between pro-viral pathways and protective responses within the host cells.

Because the gastrointestinal epithelium represents the primary target for initial EV-A71 infection, the differentiated enteric epithelial cells from human enterocyte cell line C2BBe1 were chosen in this study as the experimental model. An alternative choice of platform is the human induced pluripotent stem cells (iPSCs), which can be differentiated into multiple somatic cell lineages that ultimately form three-dimension organoids [25,26]. Given that intestinal organoids derived from iPSCs mimic the environments of the gastrointestinal tract, they are presumably a viable in vitro model for studying the enteric pathogen infection and replication [27]. However, because of the heterogeneity in the cellular make-up of the intestinal organoids, which comprise variable proportions of enterocytes, goblet cells, and Paneth cells, the viral susceptibility and host response to infection could vary to a large extent among the cell types and even between experiments. The differentiated intestinal C2BBe1 cells used by the present study have invariably apical villin localization and are more homogenous in cell state. The results obtained from this model may thus be representative of the infection occurring in differentiated enterocytes.

In our initial exploration of the possible expression correlation between known EV-A71 receptors and entry facilitation factors, such as SCARB2, PSGL-1, ANXA2, and WARS, we made the surprising finding from the single-cell transcriptome analysis that ANXA2 was abundantly expressed in differentiated C2BBe1 cells and closely correlated with the viral RNA expression in the infected population. These results are in line with the scenario that ANXA2 is an important factor for EV-A71 infection in the enteric cells; however, its upregulation in response to the EV-A71 infection remains presently unresolved. In addition to viral receptors, previous studies have shown that the release of EV-A71 particles from enterocytes may be mediated by exosomes [24]. Interestingly, we discovered that CD63, an endosome-specific marker protein, was highly expressed in differentiated C2BBe1 cells (Supplemental Fig. S10), and significantly correlated in expression with EV-A71 RNA levels at single-cell resolution. These concerted expression profiles between host factors and viral RNA provided strong support the notion that transcriptomic events accompanying EV-A71 infection might facilitate virus particle entry and release, and consequently the overall infection propagation.

A number of studies have demonstrated that infection by EV-A71 could down-regulate cellular antiviral response. For example, activation of type I IFN expression is reportedly interfered by proteins encoded by EV-A71 [28–30]. A recent research conducted with primary human enteroids additionally showed that IFN-related antiviral signalling was not activated to a significant extent upon EV-A71 infection [13]. Consistently with these findings, our single-cell transcriptomic analysis of EV-A71-infected differentiated enteric C2BBe1 cells also indicated that IFN-β, λ1, and λ2 were not significantly induced upon infection (Supplemental Fig. S7). In addition, the analysis of EV-A71-infected group showed that expression of several ISGs, such as IFIT1, IFITM1, ISG15, and DDX58, was reduced in the vRNA-positive cells as compared to the cells without vRNA (at 12 h post-infection; Fig. 5B). These expression patterns together attested to the notion that EV-A71 may suppress antiviral response in the early stage of infection as a means to establish its replication in the primary infection site. However, while a subset of cells in the infected culture did not express vRNA, they were found to produce higher levels of ISGs in comparison to the mock-infected cells. It thus remains formally possible that antiviral response may be stimulated to a greater extent in the cells neighbouring the EV-A71-infected cells, even though the infected cells did not express sufficient amounts of type I and III IFNs. Viewed together, our current findings based on the single-cell RNA sequencing technology constitute a strong basis for illuminating the bystander effect of IFN response as part of the host–virus interaction.

The roles of lncRNAs as critical and multifaceted regulators of gene expression have been extensively documented. Despite the lack of protein-coding potential, the majority of these genes are transcribed by RNA Pol II and therefore expressed with poly(A) tails, an attribute that facilitates their detection by the existing single-cell RNA-seq platform. We further relied on these data to catalogue at the single-cell resolution the dysregulated lncRNAs in the context of EV-A71 infection. Towards this end, one of the key findings of our study is the identification of CYTOR as an infection-associated, differentially expressed lncRNA and also its potential implication in promoting EV-A71 replication. Previously characterized as a regulator in tumorigenesis, CYTOR was found upregulated in multiple types of cancers, including colon and gastric cancers and nasopharyngeal carcinoma [31–33]. Our results, therefore, broadened the functional significance of this lncRNA to virus replication, and strongly evidenced a critical host non-coding factor underlying enteroviral infection.

Acknowledgments

This research was supported by grants from the Ministry of Science and Technology, Taiwan (MOST105-2314-B-182-061-MY4 and MOST107-2320-B-182-042-MY3 to BC-MT) and the Chang Gung Memorial Hospital (CMRPD1H0023, CMRPD1H0263, and BMRP960 to BC-MT; BMRPC09 to RL-K). This work was also financially supported by the Research Center for Emerging Viral Infections from the Featured Areas Research Center Program within the framework of the Higher Education Sprout Project by the Ministry of Education (MOE) in Taiwan and the Ministry of Science and Technology (MOST), Taiwan (MOST108-3017-F-182-001 and MOST109-2634-F-182-001). The bioinformatics analyses were performed by the Bioinformatics Core Laboratory, Molecular Medicine Research Center, Chang Gung University, Taiwan (grant CLRPD1J0012).

Supplementary material

Supplemental data for this article can be accessed here.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Chang Gung Memorial Hospital, Linkou [BMRP960]; Chang Gung Memorial Hospital, Linkou [BMRPC09]; Chang Gung Memorial Hospital, Linkou [CMRPD1H0263]; Chang Gung Memorial Hospital, Linkou [CMRPD1H0023]; Chang Gung University [grant CLRPD1J0012]; Ministry of Science and Technology, Taiwan [MOST108-3017-F-182-001]; Ministry of Science and Technology, Taiwan [MOST109-2634-F-182-001]; Ministry of Science and Technology, Taiwan [MOST107-2320-B-182-042-MY3]; Ministry of Science and Technology, Taiwan [MOST105-2314-B-182-061-MY4].