Parvovirus B19 and Cellular Transcriptome Dynamics in UT7/EpoS1 Cells

1. Introduction

In the Parvoviridae family, Parvovirus B19 (B19V) is a widely diffuse human virus with an ample pathogenic potential. B19V is characterized by a selective tropism for progenitor cells in the erythroid lineage (EPC), in which it exerts a cytotoxic effect with consequent blockade of erythropoiesis, leading to the typical pathologies of the hematopoietic compartment, such as erythroid aplastic crises and pure erythrocyte aplasia. The tropism for erythroid progenitor cells in the bone marrow is in strict dependence on the differentiation stage and proliferation rate of this cell population. This restriction is the result of a combination of cell susceptibility, linked to the presence on the cell membrane of specific receptor moieties for B19V virions, and permissiveness, linked to erythroid-specific intracellular signaling pathways. B19V can infect other cell types, notably endothelial or connective tissue cells; in these, internalization may not depend on receptor interaction, the replicative cycle is not normally productive, and consequences are related to the induction of an inflammatory response [1,2].

Studies on B19V have been hampered by its demanding in vitro growth requirements. Apart from primary erythroid progenitor cells (EPC) cultures, differentiated from peripheral blood mononuclear cells [3,4], a few cell lines are permissive for B19V and can be useful as model system for studying virus-cell interactions. The UT7/Epo cell line, and in particular the subclone UT7/EpoS1 [5,6], is widely used, although allowing only a restricted pattern of infection and a limited support to viral replication [7]; understanding the cell-type specific determinants of restriction to viral replication and the specific impact of virus on the cell expression profile, are topics of relevance [8].

High Throughput Sequencing (HTS) techniques are powerful tools for thorough investigation of genomic and transcriptomic layers in biological system, such as virus-cell systems. In particular, mRNAseq techniques can be used in a dual targeting approach, to trace both the dynamics of the viral transcriptome (outlined in Figure 1 for B19V) and the modification of host-cell transcriptome as a response to virus-induced stress. By these means, we previously developed a dedicated experimental and bioinformatics pipeline to investigate the B19V - EPCs system [9]; we now applied such experimental workflow to investigate the virus expression profile and the induced effects on host cell expression profile in B19V infected UT7/EpoS1 cells. A focus on UT7/EpoS1 cells is justified by their relevance as a standard cellular model for investigation on B19V lifecycle. In comparison to EPCs, UT7/EpoS1 cells constitute a more homogenous cell population, offering relevant terms for appraisal of variation in viral and virus-induced cellular expression profiles.

2. Materials and Methods

2.1. Cells and virus

Cells. UT7/EpoS1 cells were cultured in IMDM supplemented with 10% FBS, 2 U/ml Epo. The cells were maintained at 37°C and 5% CO2, at a maximal density of 10⁶ cells/ml.

Virus. B19V was obtained from a cloned synthetic genome, as described [12]. For infection, UT7/EpoS1 cells were incubated at a density of 10⁷ cell/mL, in the presence of B19V to a multiplicity of infection (moi, expressed as geq/cell) of 10² geq/cell, for 2 h at 37 °C. After removal of inoculum virus, cells were incubated at 37 °C in 5% CO₂ in complete growth medium, at an initial density of 10⁶ cells/mL.

2.2. Quantitative molecular analysis

Sampling and nucleic acids purification. Sampling was carried out for uninfected UT7/EpoS1 cells, as controls, and at 2-, 16- and 48-hours post-infection (hpi), in triplicate samples. For each sample, equal amounts of cell cultures, corresponding to 1.5x10⁵ cells, were collected. Pelleted cells were then processed by the Maxwell Viral Total Nucleic Acid kit on a Maxwell MDx platform (Promega), to obtain a purified total nucleic acid fraction in elution volumes of 150 µL.

qPCR and qRT-PCR. A quantitative evaluation of target nucleic acids was carried out by qPCR assays in a Rotor-Q system (Qiagen). For the analysis of B19V DNA, aliquots of the eluted nucleic acids (corresponding to ~500 cells) were directly amplified in a qPCR assay (Maxima SYBR Green qPCR Master Mix, Thermo Scientific). For the analysis of B19V RNA, parallel aliquots were first treated with the Turbo DNAfree reagent (Ambion) before amplification in a qRT-PCR assay (Express One-step SYBR GreenER Kit, Invitrogen). Standard cycling programs were used, followed by a melting curve analysis to define the Tm of amplified products. Primer pairs (AppendixTable A1) were selected to allow quantitation of viral DNA, total viral RNA, and selected subsets of viral RNA. Quantitative evaluation of target was obtained by absolute quantitation on external calibration curves as described [13,14].

2.3. mRNAseq analysis

Sample preparation. Collected samples were processed by the Maxwell 16 SimplyRNA Cells Kit on a Maxwell MDx platform (Promega), to obtain a RNA fraction in elution volumes of 50 µL. The amount of purified RNA was determined with Qubit 4 Fluorometer (Invitrogen, Carlsbad, CA) using a Qubit RNA BR (Broad Range) Assay Kit. The RNA integrity was assessed by agarose gel electrophoresis, and further by the Bioanalyzer RNA assay (Agilent technologies, Santa Clara, CA) before HTS mRNA sequencing.

Sequencing. mRNAseq was carried by IGA Technology Services (Udine, Italy). Libraries were prepared using a CORALL Total RNA-Seq Library Prep Kit and a RiboCop rRNA Depletion Kit to remove ribosomal RNA from the samples. Final libraries were checked with both Qubit 2.0 Fluorometer and Agilent Bioanalyzer DNA assay. Sequencing was performed on paired-end 150 bp mode on NovaSeq6000 (Illumina, San Diego, CA). Raw data containing the base callings were demultiplexed and converted to FASTQ files by using Illumina’s Bcl2Fastq 2.20 software, and adapter sequences were masked with Cutadapt v1.11.

2.4. Data Analysis

Reads Trimming. Trim Galore! (v0.6.10) was used in paired end mode to remove the adapters and trim off the low-quality bases found at the ends of each read. FastQC (v0.12.1) was utilized to check the quality scores of the samples before and after trimming.

Reads Mapping. Trimmed reads were aligned to the synthetic B19V EC Genotype I consensus sequence (GenBank KY940273.1) with HISAT2 (v2.2.1) [15]. Aligned reads were then converted to BAM files using SAMtools (v1.22) [16] before sorting and indexing.

Features Mapping. Reads mapping to splice and cleavage-polyadenylation sites were extracted by string search using seqkit (v2.10.0) [17] (AppendixTable A2). Captured expression signals were collected into FASTQ files and later aligned to the B19V EC Genotype I consensus sequence (GenBank KY940273.1) with HISAT2. Sorted and indexed BAM files were used to generate coverage graphs with a custom Python script leveraging the pysam library (v0.23.3) together with matplotlib (v3.10.5) to generate the graphs.

Transcript Count. Salmon (v1.10.3) [18] was used to quantify the transcript abundance from the trimmed FASTQ files. GENCODE (release 48) [19] dataset was utilized. The quantified counts were then imported into the R statistical environment (v4.5.1) using the tximport R package (v1.36.0) [20] and the metadata file of the transcriptome containing the translation from transcript identifier to gene symbol. Around 22% of all transcript identifiers lacked an associated gene symbol; the affected transcripts were thus excluded from the downstream analysis.

Statistical Design. A multi-factorial linear model was designed, incorporating the variables of hours post-infection (hpi) and infection status, along with their interaction terms.

Within edgeR’s framework [21], the expected read count ${\mathsf{\mu }}_{gi}$ for gene g in sample i is modeled as a log-linear model of the form

$$\mathit{log}{\mathsf{\mu }}_{gi}=x_i^T{\mathsf{\beta }}_g+\mathit{log}{L}_i$$

(1)

where $x_i$ is a covariate vector encoding the experimental design for sample i, ${\mathsf{\beta }}_g$ is a vector of coefficients representing the experimental effects and log-fold-changes, and $L_i$ is the effective library size for sample i.

The observed read count $y_{gi}$ is assumed to follow a mixture distribution across biological replicates whose variance relates to the mean quadratically

$$var\left(y_{gi}\right)={\mathsf{\sigma }}_g^2{\mathsf{\mu }}_{gi}+{\mathsf{\psi }}_g{\mathsf{\mu }}_{gi}^2$$

(2)

where ${\mathsf{\sigma }}_g^2$ represents the technical variability, or quasi-dispersion, and ${\mathsf{\psi }}_g$ captures biological variability, commonly referred in the field as the Biological Coefficient of Variation (BCV) [22]. Synchronic or diachronic contrasts, which are linear combinations of the statistical model coefficients, were defined in order to perform comprehensive tests for gene dysregulation. Design matrix and contrast creation was performed using the limma R package (v3.64.1) [23].

Differential Expression Testing and Modelling. Differential Gene Expression Analysis (DGEA) was carried out in the R environment using the edgeR (v4.6.2) package [24]. The gene-level count matrix was first imported into a DGEList object which encapsulates the raw count data as well as the sample grouping information and the experimental design matrix. Genes with a low expression level across all samples were removed to reduce the background noise in the data. Library sizes normalization was carried out using Trimmed Mean of M-values (TMM).

Exploratory Data Analysis. Principal Component Analysis (PCA) was performed on the log-transformed Counts per Million (logCPM) gene counts and visualized using the factoextra R package (v1.0.7). Multidimensional Scaling (MDS) was applied and visualized with the plotMDS function of edgeR in order to assess the expression profiles of the different samples and replicates.

Model Fitting. The BCV of each gene was computed and compared to its average log-transformed CPM value with the plotBCV function of the edgeR package. The robust quasi-likelihood negative binomial generalized log-linear model from edgeR was trained using the glmQLFit function in order to determine profiles of differential gene expression. Subsequently, the quasi-likelihood dispersion of the model was plotted using the plotQLDisp function to diagnose whether the model was able to describe the data correctly or not.

Contrast Testing. Differential gene expression in different conditions was tested by evaluating individual contrasts under the conditions of (i) no log-fold change (logFC) threshold, employing glmQLFTest, and (ii) with a logFC threshold of 1.6 using glmTreat. The results from both methods were processed with decideTests, a function within edgeR which classifies genes as upregulated, downregulated or not significantly dysregulated; a p-value < 0.05 was used as the threshold for significance. The default Benjamin & Hochberg (FDR) procedure was employed to correct for multiple testing. Visualizations were created in the form of volcano plots using the ggplot2 R package (v3.5.2) and in the form of heatmaps and UpSet plot visualizations using the ComplexHeatmap R package (v2.24.0) [25].

Gene Set Enrichment Analysis. Competitive gene set testing was performed using the Camera method as implemented in the edgeR package [26], to evaluate the dysregulation of biological pathways found in the MSigDB hallmark collection [27] whilst accounting for inter-gene correlation. Lollipop plots were created with ggplot2.

3. Results

3.1. Course of infection of B19V in UT7/EpoS1 cells

The course of B19V infection in UT7/EpoS1 cells was monitored at 2, 16 and 48 hpi by qPCR analysis, to quantitate variation in the amount of viral DNA, and by qRT-PCR, to quantitate variation in the amount of total (mRNA1-5) and of relevant subsets of viral mRNAs, in particular the proximally-cleaved, unspliced mRNAs (mRNA1) or the distally cleaved, spliced mRNAs (mRNA3-5). The amount of viral DNA detected in the cell population remained stable from to 2 to 16 hpi, and showed a net increase of 1.2 Log from 16 to 48 hpi, confirming the productive although restricted infection pattern in UT7/EpoS1 cells. Total viral mRNA, barely detectable at 2 hpi, showed increasing amounts at 16 and 48 hpi, also confirming the productive pattern of infection (Table 1). By comparing the abundance of mRNAs subsets between these two last time points, it was confirmed a relatively higher abundance of proximally cleaved unspliced transcripts, coding for NS1 protein, at 16 hpi, and a relative higher amount of distally cleaved spliced transcripts, coding for VP1/2 and 11 kDa protein, at 48 hpi. Therefore, the biphasic early/late expression pattern of B19V expression was confirmed in the presented experimental setting, constituting a framework for interpretation of mRNAseq data.

The outcome of infection in the cell population was further confirmed by FISH detection of viral nucleic acids at 48 hpi [28]. The fraction of positive, productively infected cells was evaluated at about 2%, a low value in accordance with the known restrictive pattern of infection in UT7/EpoS1 cells (Figure 2). Therefore, subsequent interpretation of mRNAseq data, at least for viral mRNAs, should be considered as relative to such subset of cells.

Figure 2. Flow-FISH for viral nucleic acids in control and infected UT7/EpoS1 cells at 2 and 48 hpi; FITC labelling, X-Axis. FSC, Y-Axis. FISH positive cells (inlet) are in UR quadrant.

Figure 2. Flow-FISH for viral nucleic acids in control and infected UT7/EpoS1 cells at 2 and 48 hpi; FITC labelling, X-Axis. FSC, Y-Axis. FISH positive cells (inlet) are in UR quadrant.

3.2. mRNAseq analysis

For mRNAseq experiments, analyzed samples included uninfected UT7/EpoS1 cells as a baseline control, and B19V infected cells collected at 2, 16, 48 hpi. Each sample was processed in triplicate for RNA purification, and mRNAseq analysis was conducted on about 1 µg RNA per sample. The rough output yielded 2.9-3.5 x10⁷ reads per sample.

Read mapping was first conducted on the synthetic B19V reference genome (GenBank KY940273.1). No reads from the uninfected control samples mapped to the B19V genome, as expected. Less than 10² reads per sample were mapped to B19V at 2 hpi; 3.1-3.5 x10⁴ reads mapped to B19 genome at 16 hpi (0.01% of total); 1.3-1.4 x10⁶ reads mapped to B19V at 48 hpi (0.53% of total). For each sample, mapping of residual reads to the reference human genome (build GRCh38.p14) returned between 69-80% correctly mapped reads.

3.3. Viral transcriptome

Mapping of reads to B19V genome was visualized using a custom Python script leveraging the pysam library together with matplotlib (Figure 3).

At 2 hpi, the scattered mapping of very few reads testified the early onset of viral mRNA synthesis, but did not allow any further characterization of the viral transcriptome. At 16 and 48 hpi, the increasing number of mapped reads allowed for such characterization. Reads were distributed on the genomic template, in different abundances depending on the genomic region and reported inclusion in exons or introns, in different patterns for the different time points. The first increase in read number at 16 hpi correlated to a prevalent mapping to the left-side genomic region, within the proximal cleavage-polyadenylation sites. The substantial increase in reads at 48 hpi correlated to a prevalent splicing of the first intron, coupled to an increased, although irregular, representation of the right-side genomic region, up to the distal cleavage-polyadenylation site (Table 2).

A closer inspection was conducted to assess the relevance and pattern of usage of known mRNA processing sites: 1) start of transcription at nt. 530; 2) donor/acceptor splice sites at nt. 586/2089-2209 (D1/A1.1-A1.2); donor/acceptor splice sites at nt. 2363/3224-4883 (D2/A2.1-A2.2); the pAp1 cleavage-polyadenylation site at nt. 2842, the pAp2 site at nt. 3142; the pAd distal cleavage polyadenylation site at nt. 5189. To this purpose, reads spanning reported splice junctions and cleavage-polyadenylation sites were specifically selected and mapped, and their relative abundance determined (Figure 4).

Considering the splicing process, the first intron showed a prevalent frequency of splicing events, increasing from early to late time points (13% vs. 5% unspliced vs. 87% to 95% spliced transcripts), and usage of proximal to distal acceptor sites. The second intron showed a more balanced and constant frequency of splicing (~40% unspliced transcripts) and usage of proximal and distal acceptor sites. Considering the cleavage-polyA process, usage of the pAp1 site resulted prevalent at both early and late time-points, the pAp2 site was not represented over background, and the pAd site showed increasing frequency at late time point (Table 3).

3.4. Cellular transcriptome

Analysis of cellular transcriptome compared triplicate samples of not-infected to infected cells collected at 2, 16 and 48 hpi. Principal Component Analysis (PCA) (Figure 5A) performed on the log-transformed CPM values showed significant gene expression differences between sample groups, coupled to low variability within groups, a scenario also confirmed by Multidimensional Scaling (MDS) (Figure 5B).

Figure 5. A. PCA plot; the two main dimensions account for 61% of total variation. B. MDS plot; the two main dimensions account for 52% of total variation. In both PCA and MDS plots, the samples cluster according to the hpi variable, while the tight clustering of the replicates with one another certifies the similarity between them.

Figure 5. A. PCA plot; the two main dimensions account for 61% of total variation. B. MDS plot; the two main dimensions account for 52% of total variation. In both PCA and MDS plots, the samples cluster according to the hpi variable, while the tight clustering of the replicates with one another certifies the similarity between them.

Gene count dispersion assessed with a Biological Coefficient of Variation (BCV) plot (Figure 6A), comparing the BCV to the average log-transformed CPM value of each gene, confirmed that experimental data provided enough statistical power to distinguish differential expressed genes (DEGs) with high reliability. In the Quarter-Root Mean Deviance plot (Figure 6B), the tight clustering around a relatively smooth trend line indicates that the model was able to estimate correctly the quasi-likelihood dispersion values and that the variance structure was well estimated by the applied statistical model.

Differential Gene Expression Analysis (DGEA) identified dysregulated genes across the tested synchronic (each time point vs. control) or diachronic (time point vs. time point) contrasts, visualized as volcano plots, comparing the log-transformed fold change of each gene to its negative log-transformed p-value (Figure 7). These plots highlight the effects of the progress of infection. Synchronic contrasts (A-C) show a progressive increase in gene dysregulation, from the 2 hpi sample showing minimal dysregulation to the 48 hpi exhibiting the most differential gene expression. Diachronic contrasts (D-E) show a substantial difference in the set of differentially expressed genes starting from 2 hpi compared to the 16 and 48 hpi samples, whilst comparison of the 16 hpi to the 48 hpi conditions reveal a subtle difference in gene expression, suggesting that most of virus-induced variation in gene expression is set at early times post-infection.

This scenario is also confirmed by UpSet plots (Figure 8), that visualize the size of the intersections among the sets of differentially expressed genes identified across synchronic or diachronic contrasts. The intersection profiles observed in the synchronic contrast plot (A) reveal that 343 genes (~12% of all dysregulated genes) are shared between the 16 hpi and 48 hpi condition, 186 upregulated and 157 downregulated. The next greatest intersection is between the 2 hpi and 48 hpi condition in both up- and down- regulation. Interestingly, 35 genes (~3% of all downregulated genes) are consistently downregulated across all conditions. Intersection in the diachronic contrasts (B) show that most genes (880 genes, ~30% of all dysregulated genes) are common between the 2-16 and 2-48 hpi, whilst additional 203 (~7% of all dysregulated genes) add between 16 and 48 hpi. The number of genes dysregulated in contrasting sense is constantly low.

To provide an insight on the gene dysregulation during the course of infection, the fifty most dysregulated genes were visualized in a heatmap representation (Figure 9) in which each row represents a gene and each column represents a tested contrast. Results for both synchronic and diachronic contrast indicate a clear separation between two main clusters, either down- or up-regulated.

Gene Set Enrichment Analysis (GSEA) was carried out on the list of recognized genes using CAMERA (Competitive Gene Set Tests for Digital Gene Expression Data) on the hallmark collection of the Molecular Signature Database (MSigDB) and visualized as a lollipop graph (Figure 10). The top sets were selected based on the enrichment of the set among the different contrasts and their p-value score. For synchronic contrasts, GSEA mainly returned down-regulation of gene sets, whilst only a few gene sets showed up-regulation. Of these, at 2 hpi, unfolded protein response; at 16 hpi, mitotic spindle and glycolysis; at 48 hpi, hypoxia and heme metabolism. On the opposite, as emerging from diachronic contrasts, more gene sets were upregulated from 2 to 48 hpi, including hypoxia, heme metabolism, glycolysis, estrogen response and bile and fatty acid metabolism.

Finally, within the gene sets previously identified, an exploration of interaction networks was conducted by using the STRING database and a clustered analysis with the Markov Cluster Algorithm (MCL) [29]. Results are shown in Appendix A, Table A3 for the most relevant clusters, according to samples (synchronic or diachronic contrasts). By reference to the Reactome database [30], the principal pathways involved were: at 2 hpi, immune system signaling by interleukins and cellular responses to stress; at 16 hpi, cytokine signaling, immune system signal transduction, and cell cycle regulation; at 48 hpi, immune system signal transduction and cell cycle regulation; in the 2-48 hpi time interval, mainly cell cycle regulation and cellular responses to stress.

3.5. Cellular Transcriptome, Comparison of UT7/EpoS1 to EPCs

Direct availability of a transcriptome analysis of B19V infected erythroid progenitor cells (EPCs) [9] allowed for a comparison of the virus-induced expression profile perturbations in these two cellular systems. A simple comparison of UT7/EpoS1 versus EPCs, at the same time points of 2, 16 and 48 hpi, only showed a substantial divergence of the two cellular systems. To elucidate the variations specifically induced by viral infection, the model had to determine the differential variations in the expression profile in the infected and in the basal, not infected, conditions, for both cellular systems (Figure 11). This allowed for individuation of convergent or divergent dysregulated gene sets, eliminating a background of gene sets not significantly affected by the virus (Figure 12).

Comparing the two systems, most genes showed different expression patterns (upset plot, row set size), a result highlighting the difference in the two cell population. A minor fraction of genes showed concurrent virus-induced dysregulation (upset plot, column set size), and most of these coherently across the different time points, either upregulated or downregulated in the two systems. GSEA was carried out on MSigDB, followed by exploration of interaction networks by using the STRING database and a clustered analysis with the Markov Cluster Algorithm (MCL). Results are shown in Appendix A, Table A4 for the most relevant clusters. With reference to the Reactome database, the principal pathways showing differential regulation between UT7/EpoS1 and EPCs were: at 2 hpi, mitotic G1 phase and G1/S transition, transcription by RNA polymerase II; at 16 hpi, transcription regulator activity, response to cytokine, Interferon signaling; at 48 hpi, cell cycle, extracellular matrix organization, cytokine signaling, Interferon alpha/beta signaling, Interleukin signaling, DNA Repair.

4. Discussion

Investigation of virus-cell interactions can take substantial advantage from the implementation of HTS techniques in addition to the commonly used quantitative molecular techniques. Potential advantage of HTS techniques are their unbiased targeting and output, open to genome- and transcriptome-wide investigation. Specifically, when investigating virus-cell systems, experimental output includes both viral and cellular mRNAs and is therefore especially suited to comparative analysis and correlation studies between viral transcription and induced modifications in the cell expression profile at population level.

In our present work, we focused on the system formed by B19V and UT7/EpoS1 cells, a commonly used cell line for in vitro studies on B19V lifecycle. UT7/EpoS1 cells are a subclone of the parental UT7/Epo cells, an erythropoietin-committed sublineage originally derived from the multipotent myeloblastoid UT7 cells [5,6]. Different subclones derived from the parental UT7/Epo show different degrees of susceptibility and permissiveness to B19V infection, correlating to the commitment to progression though cell cycle [8]. Thus, while population heterogeneity may not be apparent from cell phenotype, it is however be assumed as a basic property accounting for the restrictive properties towards B19V, since as only a small percentage of cell actually can support a productive infection. Thus, all information obtained either by quantitative molecular techniques or HTS techniques, in the absence of single cell analysis, is necessarily averaged over the cell population.

Within this framework and limits, our investigation yielded information on virus and cell transcriptome modulation in the course of B19V infection in the UT7/EpoS1 cells. mRNAseq analysis of the viral transcriptome returned information in agreement with what obtained by quantitative molecular methods, including information on differential transcript abundance and post-transcriptional processing [13,14]. The leader sequence from the start of transcription is normally detected at high abundance, and splice boundaries are neatly defined for both first and second intron, with respective alternative processing patterns. Concerning the left-hand genomic cassette, mRNA1 transcripts, encoding the NS1 protein, are expressed during all stages of infection, at higher relative abundance at early stages and lower relative abundance at late stages. mRNA2 transcripts, which may encode the putative 7.5 kDa and a 9 kDa protein, are highly abundant during all infectious phases as previously known. According to the string analysis, the first proximal polyadenylation site pAp1 is more easily detected than the second proximal site pAp2, which does not emerge over background. Concerning the right-hand genomic cassette encoding for VP and 11 kDa proteins, the increase at late compared to early times is less evident, and usage of the pAd cleavage-polyadenylation site less defined. More accurate mapping is likely prevented by the sequencing strategy employed, which did not yield a uniform coverage over this distal region of mature transcripts.

Information on cellular transcriptome compared variation induced by the virus in the course of a replicative cycle, either with respect to a basal, not-infected state for single time points (synchronic variations); or with respect to subsequent different time points (diachronic variations). Within a very complex interaction network, downregulation of genes involved in the immune/cytokine/interleukin response is prominent from earlier times through the whole time course of infection (2 hpi, cluster #1; 16 hpi, cluster #1; 48 hpi, (cluster #1 and #2). On the opposite, upregulation of genes involved in cell stress response is found at 2 hpi (cluster #2), and genes involved in cell cycle regulation are modulated at 16 hpi (cluster #2) and 48 hpi (cluster #3). Likely, the early response is distributed over the whole cell population and may contribute to definition of the restrictive characteristics of this cell population; response at later times may thus be restricted to the productively infected subset of cell, and shape the cellular environment to determine the outcome of infection. Given the inherent heterogeneity of UT7/Epos1 cells and the restrictive pattern of permissiveness to B19V, single-cell analysis will necessarily be required to better dissect the viral-induced modulation of cellular environment.

The present dual transcriptome analysis conducted on UT7/EpoS1 cells can be compared to the recent analysis carried out on differentiated EPCs, that constitute the cellular system more closely representative of the natural target cells in bone marrow [9]. Viral transcriptome showed quite comparable profiles through the time points, thus corroborating the validity of UT7/EpoS1 cells as a model system supporting B19V replication. On the other hand, concerning the cellular transcriptional landscape, the differences in the two systems are prevalent, such that some comparative information can be obtained only by including in the analytical model all variables – cell type, infection status, and time point post-infection. Given the interaction network, UT7/EpoS1 show a relative upregulation of genes involved in cell cycle regulation (e.g., 2 hpi, cluster #1, and 48 hpi, cluster #1); and a relative under-expression of genes involved in cytokine and interleukin signaling (e.g., 16 hpi, cluster #2, and 48 hpi, clusters #3, #5). Differences between the two systems thus involve key cellular processes, where the impact of virus may differ substantially; thus, caution should be exerted when extrapolating results obtained in the myeloblastoid UT7/EpoS1 cells to the natural target EPCs.

5. Conclusions

In our present work, a characterization of both viral and cellular expression profile in the course of B19V infection of UT7/EpoS1 was obtained, reconstructing the viral transcriptome, and highlighting viral-induced variations in the cellular transcriptome. A comprehensive analysis of the modulation of the expression profile in a model cell population is presented, also providing a comparison to the in vitro differentiated EPCs, which are most closely representative of the target erythroid cells in bone marrow. Similarities and differences in the two systems emerged, thus posing a note of cautions when extrapolating experimental data obtained in UT7/EpoS1 to EPCs. Further, compared to bulk techniques, a single-cell analysis approach, now more easily attainable, will better elucidate the complex virus-cell relationship and the impact of virus infection on target cells. In turn, this will allow a better comprehension of the pathogenic mechanisms of infection and a better definition of potential antiviral strategies.