To address this task we used two datasets, GSE151513 and GSE158930. GSE151513 comprises gene-expression profiles from Calu-3 human lung adenocarcinoma cell-lines infected with SCOV2 virus at a multiplicity of infection (MOI) of 2. Poly-A RNAseq gene-expression profiles were acquired with the ′Illumina HiSeq 2500 (Homo sapiens)′ platform (15,761 of them are included in the provided original dataset as gene-filtering was applied). The cell-line samples were collected at six hours post-infection (hpi) time-points: 0, 1, 2, 3, 6 and 12 hours, in triplicates, (n=18, three replicates for each infected sample). No normalization was performed as the original data were already normalized². The results of the study are published in (57). We applied gene filtering on the available RNAseq data, removing 40% of gene-transcripts, ranked by their maximum expression level across all samples, leaving 9,457 gene-transcripts for further analysis. The results from the experiment conducted with the GSE151513 dataset are presented in section 3.1.1, and demonstrate a two-stage core molecular fingerprint that contrast between early and late SCOV2 infection stages. In order to further confirm our findings, we used the GSE158930 RNAseq dataset. GSE158930 comprises gene-expression profiles from ALI (air-liquid interface) cultures of primary human bronchial epithelial cells (HBECs) infected with SCOV2 virus at a MOI of 1. RNAseq gene-expression profiles were acquired with the ′Illumina NextSeq 500′ platform (57,905 Ensembl gene-transcripts). The cell-line samples were collected at five hpi time-points: 4, 24, 48, 72 and 96 hours (in du-/tri- or quadr-uplicates, n=19 infected samples). The results of the study are published in (58). We also applied gene filtering, removing gene-transcripts with CPM<5 across all samples, leaving 9,833 gene-transcripts for further analysis. In addition, sample filtering was also applied, keeping only the duplicate samples (i.e., replicates 1 and 2, as not all samples were profiled in tri- or quadru-plicates), leaving 10 infected samples to analyze. The filtered gene-expression counts were transformed using the variance stabilizing transformation (VST) method (59). The results from the experiment with the GSE158930 dataset are also presented in section 3.1.1.

To contrast the molecular profiles underlying SCOC2 and SCOV1 infections we used two datasets, GSE33267 and GSE148729. GSE33267 includes the gene-expression profiles of Calu-3 cell-line samples infected with the wild-type (icSARS) or the DORF6 (mutant not expression ORF6 protein) SCOV1 strains at a MOI of 5. The gene-expression profiles were acquired with the ′Agilent Whole Human Genome Microarray 4x44K G4112F (Feature Number version)′ microarray platform (33,631 gene-probes). Samples were collected in triplicates at eleven different hpi time-points, namely: 0, 3, 7, 12, 24, 30, 36, 48, 54, 60 and 72 hours. The results from this study are published in (60). No data normalization was performed as the original gene-expression profiles were already quantile normalized and log2 transformed. We retained just the icSARS/SCOV1 samples (n=33); the mock-treated and DORF6 samples were discarded. The feature numbered probes were mapped to human gene nomenclature committee/HGNC³ official gene symbols, relying on the respective iDEP’s mapping process. Gene filtering was also applied where, 40% of genes ranked by their maximum expression level across all samples were discarded, leaving 10,353 genes for further analysis. The results from the experiment conducted with the GSE33267 dataset are presented in section 3.1.2, and show a similar two-stage core molecular fingerprint underlying both SCOV2 and SCOV1 infections. Confirmation of the findings was done using the GSE148729 dataset. GSE148729 comprises RNAseq gene-expression profiles of Calu-3 human cell-line samples infected with SCOV1 (Frankfurt strain) and SCOV2 (patient isolate BetaCoV/Munich/BavPat1/2020| EPI_ISL_406862) viruses at a MOI of 0.33, using bulk and single-cell polyA-RNA, smallRNA, and totalRNA sequencing. In our experiments we focus on the total-RNAseq gene-expression data. Gene-expression profiles were acquired with the ′Illumina NextSeq 500/HiSeq 4000 (Homo sapiens)′ platform (40,648 Ensembl gene-transcripts) over three hpi time-points: 4, 12 and 24 hours, in duplicates (n=12, 6 samples for each virus infection). The results of this study are published in (61). Gene filtering was also applied where, gene-transcripts with CPM<5 across all samples were discarded, leaving 10,535 Ensembl gene-transcripts for further analysis, with their count values being VST transformed. The results from this experiment are also presented in section 3.1.2.

To address this task, we used the GSE47960 dataset that is widely used in various relevant studies. It contains the gene-expression profiles of human airway epithelial (HAE) cell cultures infected with the wild-type (icSARS) SCOV1 strain (n=34), and the H1N1 influenza strain (n=20), with samples collected at nine different hpi time-points: 0, 12, 24, 36, 48, 60, 72, 84 and 96 hours for SCOV1, and seven hpi time-points: 0, 6, 12, 18, 24, 36 and 48 hours for H1N1. The larger extend of SCOV1 hpi time-points is expected, as the respective immune/defense responses, in contrast to INFL/H1N1, are delayed in SCOV1 infection. This is mainly due to the larger extent of the respective incubation periods; a mean of 2 days for the INFL/H1N1 2009 pandemic, 2-7 days for SCOV1 and 4-12 days for SCOV2 (62). Gene-expressions were profiled with the ′Affymetrix Human Gene 1.0 ST Array′ microarray platform (32,067 gene-probes). The results of the study are published in (63). In our experiment we kept just the samples infected with the icSARS/SCOV1 strain; the samples infected with other strains (SARS-dORF6 and SARS-BatSRBD) were discarded. No normalization was performed, as the gene-expression profiles in the original dataset were already quantile normalized, and no gene-filtering was applied. The results of this experiment are presented in section 3.1.3 and demonstrate a stronger, compared to SCOV1, immune/defense response profile during the early H1N1 infection stage.

To address this task we used two datasets, GSE152075 and GSE156063. GSE152075 comprises RNAseq profiles (acquired with the Illumina ′NextSeq 500′ gene-level platform; 35,784 genes) of nasopharyngeal (NP) swabs from 430 SCOV2 infected individuals and 54 negative controls. In our experiments we retained and focused just on the adult samples (age >= 20). SCOV2 viral-load were assessed by N1Ct, the PCR cycle threshold (Ct) of the SCOV2 nucleocapsid N1 target. As suggested in the original study publication (64), N1Ct values were discretized into three intervals that reflect respective viral-load levels: LOW (25.1 ≤ N1Ct ≤ 30.5, average N1Ct = 26.4, n=59) and HIGH (12.3 ≤ N1Ct ≤ 18.0, average N1Ct = 16.2, n=72); the samples with medium viral-load levels (18 ≤ N1Ct ≤ 25) were disregarded as we are interested to contrast between extreme viral-load levels. Gene filtering was applied where, genes with CPM<2 in at-least 59 samples (the number of samples with low viral-load) were filtered-out, leaving 8,130 genes for further analysis. The original gene read-counts were VST transformed. The results from the experiment with the GSE152075 dataset are presented in section 4.1, and indicate that low viral-load cases exhibit a suppressed immune response profile similar to the one underlying the early SCOV2 infection stage. Confirmation of findings was performed using the GSE156063 dataset. GSE156063 comprises RNAseq profiles acquired with the ′Illumina NovaSeq 6000′ ensemble gene-transcript level (15,979 Ensembl transcripts are included in the provided original dataset as gene-filtering was applied) of nasopharyngeal/oropharyngeal (NP/OP) swabs from 234 patients with acute respiratory illnesses (ARIs), infected either by SCOV2 or by some other virus or bacteria. The results of the study are published in (65). The infection type for each sample was acquired early in the infection stage by metagenomic sequencing (mNGS). mNGS integrates infection transcriptional signatures, and has proved its significance as a diagnostic tool for the assessment of ARIs, including those caused by SCOV2 infection (66–69). Viral-load is quantified in reads-per-million (RPM) values. In our experiments, only the samples with log₁₀RPM>0 were retained for further analysis. The log₁₀RPM values were quantile discretized into four nominal values, namely: too-high (average log₁₀RPM=5.4, n=17), high-to-medium (4.3, n=16), low-to-medium (2.3, n=16), and too-low (1.15, n=17); samples with medium viral-loads (average log₁₀RPM=3.0) were disregarded. Discretization was performed using scikitlearn’s KBinsDiscretizer python implementation following a quantile strategy. Too-high and high-to-medium samples were assigned to the HIGH viral-load class (3.75 ≤ log₁₀RPM ≤ 5.89, average = 4.85, n=33), and the too-low and low-to-medium samples to the LOW class (0.14 ≤ log₁₀RPM ≤ 2.55, average = 1.69, n=33), respectively. Gene filtering was additionally applied where, genes with CPM<5 over all samples were filtered-out, leaving 3,551 gene-transcripts for further analysis. The original gene read-counts were VST transformed. In both experiments we focused on the SCOV2 infected samples in order to differentiate between HIGH and LOW viral-load levels. The results from the experiment with the GSE156063 dataset are also presented in section 4.1.

To address this task we utilised two datasets, GSE166190 and GSE161731. GSE166190 includes RNAseq profiles of a total of 20 SCOV2 infected individuals, 11 adult (≥20 years old) and 9 children (≤16 years old). RNAseq profiles were acquired with the Illumina ′HighSeq 4000 (Homo sapiens)′ platform; 58,825 Ensembl gene-transcripts. For each individual the days post onset of symptoms/DPOS, as reported by each individual, were recorded and divided into five intervals: interval 1 (0-5 DPOS; the early infection stage), interval 2 (6-14 days), interval 3 (15-22 days), interval 4 (23-35 days), and interval 5 (36-81 days). The original publication of this study (70) report results that contrast between the response profiles of children and adults. Gene filtering was applied where, gene-transcripts with CPM<3 in at-least half of the samples were discarded, with the read-counts being VST transformed. From them, gene-transcripts with official HGNC gene symbols were retained for further analysis (13,027 unique genes). The results from the experiment with the GSE166190 dataset are presented in section 4.2, and show that early responders exhibit a robust antiviral immune response in the early stages of infection which, inverses the standard two-stage SCOV2 infection profile. To further confirm our findings, we used the GSE161731 dataset. GSE161731 includes RNAseq profiles of a total of 77 whole-blood samples from COVID-19 patients (the RNAseq gene-expression profiles were acquired via the ′Illumina NovaSeq 6000′ platform; 60,675 Ensembl transcripts). Whole blood samples were collected between 1-35 days post infection and divided based on disease severity and time from symptom onset. In particular, as reported in the original study publication (71), the individuals were assigned to three classes according to their reported time from onset of symptoms, early (≤10 days, the EARLY responders), intermediate (>10-21 days, the intermediate/MID responders) and late (>21 days, the LATE responders), with some of them being hospitalized. Each individual was sampled at different post infection time-points. For each of the multi-sampled individuals we kept just one sample, the one that corresponds to the individual’s earliest sampling. That is, if an individual is assigned to the early class (i.e., reported in its first hospital visit, ≤10 days from the onset of symptoms), and was then re-sampled at later time-points, the sample of the first hospital visit is retained. In this way we assure that the gene expression profile of each retained sample corresponds to an early, intermediate or late response of each individual. The samples were assigned to the following classes: EARLY_NO (not hospitalized early responders, n=9), EARLY_YES (hospitalized early responders, n=5), MID_NO (not hospitalized intermediate responders, n=19), MID_YES (hospitalized early responders, n=6) and LATE_NO (not hospitalized late responders, n=6). With this process 45 unique SCOV2 infected individuals were retained for further analysis. A new class variable is introduced, ′onset_hosp′, that holds the respective combined (days from symptoms onset and hospitalization status) nominal values. Gene filtering was also applied where, the transcripts with CPM<5 in at-least 10% of the retained samples, and which are mapped to official HGNC gene symbols, were retained for further analysis (11,497 unique genes), with the read-counts being VST transformed. The results from the experiment with the GSE161731 dataset are also presented in section 4.2.

The union of seven lists of genes being included in the core molecular fingerprints induced from specific experiments aforementioned above compose a set of 52 genes that were utilized for the devise and assessment of SCOV2 and COVID-19 diagnostic and prognostic classifiers, respectively. The seven fingerprint gene used are: (a) the three gene sets that characterize the progression of SCOV2 infection and resulted from the experiments with datasets GSE151513, GSE158930 and GSE148729 (tasks 1 and 3), (b) the two gene sets that characterize the low viral-load profile of infected individuals and resulted from the experiments with datasets GSE152075 and GSE156063 (task 4), and (c) the two gene sets that characterize the profile of early responders to the infection and resulted from the experiments with datasets GSE166190 and GSE161731 (task 5) (please refer to section 5). The 52 genes were assessed for their power in the diagnosis of SCOV2 infection and the prognosis of COVID-19 when used as descriptors in respective classifier models. The diagnostic classifiers aim to differentiate between SCOV2 and other (viral or non-viral) acute respiratory illness (ARIs). Three studies and respective datasets were utilized for this task, namely: GSE156063 (65), GSE188678 (72) and GSE163151 (73). The performances of the induced classifiers were compared with the results reported in the respective publications. The prognostic classifiers aim to predict the different SCOV2/COVID-19 phenotypes with respect to different clinical outcomes (e.g., severe/critical vs. mild-moderate/non-critical) or symptomatology (symptomatic vs. asymptomatic) of real patient cases. Four studies and respective datasets were utilized for this task: GSE152418 (74), GSE178967 (75), GSE172114 (76) and GSE177477 (77). The performance of the induced prognostic models was also assessed and compared with the results reported in the respective publications.

For RNAseq datasets delivered in gene count formats a minimum CPM (counts per million) cutoff was set for a predetermined number of samples (e.g., CPM ≥ 3 in half of the samples). In all other cases, a filtering of gene-probes or genes was performed according to a preset minimum expression value. For these experiments, 40% of low expressed genes (based on their maximum value across all samples) are discarded. In addition, probes or Ensembl-transcripts not mapped to HGNC official symbols are also filtered-out. Mapping of gene-probes or Ensembl transcripts to HGNC official symbols was done either by the automatic iDEP’s mapping service or, with the help of the g:profiler server⁴. Not normalized RNAseq data were normalized with DESeq’s variance stabilizing transformation/VST method (59), as implemented in the iDEP server. The specific gene filtering and normalization processes followed in each experiment were outlined in previous section 2.1.1 where, in the presentation of each task the data-preprocessing details for the respectively used datasets are presented (refer also to Supplementary Table 1 at the end of the provided Supplement File (Supplement.pdf).

A central component of our analysis is the identification of SCOV2 infection stages and/or phenotypes based on the available gene-expression profiles. To do so, the genes of the input dataset were clustered using the k-means algorithm. The pre-defined number of clusters was determined with the Elbow method − an effective method used by various biomedical and gene-expression analysis studies (78, 79). In all experiments, k-means was applied on a set of 500-2000 most variable mean-centered genes. The genes of each cluster were then passed through an enrichment analysis process in order to identify and focus on the cluster being heavily enriched in GO-biological-processes (geneontology.org) that relate to host immune and defense responses. With a careful inspection of the PCA plot of the samples constrained on the genes included in this cluster, the clearly separable sample groups are identified. We take these groups as representatives, either of the different infection stages (e.g., early/high) or, the different viral-load phenotype classes (e.g., high/low), with the samples of each group to be assigned to the respective infection stage or viral-load class. The specifics of the designation of infection stages or phenotypic classes followed in each experiment are detailed in the respective sections and paragraphs. For a detailed example of the methodology followed for the designation of infection stages please refer to the paragraph ′Designation of SCOV2 infection stages′ in section 3.1.1 and its supporting (Figures 3B–D). Such an elaborate and thorough designation of infection stages and phenotypes, represents a novel rationally designed approach for the analysis of relevant gene-expression data. Specifically, in the case of SCOV2 infection, it highlights and enables the emergence of critical differential genes and molecular fingerprints that characterize the infection.

For each experiment, the already normalized or VST transformed gene-expression profiles of samples are assigned to respective infection stages or phenotype classes which, were then passed to a differential expression gene (DEG) analysis process. The well-known limma R package (80, 81), as implemented in the iDEP server, was utilized for this. Limma was preferred for DEG analysis as it is applied directly on transformed/normalized gene-expression profiles and not on read counts (also suggested by the developers of iDEP). DEGs are induced by setting FDR (false discovery rate) to 0.05 and varying fold-change (FC) thresholds in the different experiments. In most cases the minimum FC threshold was set equal to 2, but for datasets with a large number of induced DEGs, and in order to ease their biological interpretation, higher FC thresholds were chosen. The specifics of DEG analysis followed in each experiment are detailed in the respective sections and paragraphs. Consult also Supplementary Table 1 ate the end of the provided Supplement File (′Supplement.pdf′).

The identified up-/down-regulated DEGs were subjected to enrichment analysis for the following entries: GO-biological processes (82, 83), KEGG (www.genome.jp/kegg/pathway.html) and REACTOME (reactome.org) molecular pathways, as well as for specific gene hallmark signatures from the Human MSigDB collections (http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp) (84). The Benjamini–Hochberg adjusted hypergeometric test was used and applied for the identification of enriched entries, using the fgsea R package (85) as implemented in the iDEP server.

We heavily relied on the STRING server (string-db.org, version 11.5) (86), and in particular on the respective clustering services it offers, in order to form coherent networked gene clusters from the identified DEGs. We consider and interpret these clusters of networked genes as representative core molecular fingerprints for the designated SCOV2 infection stages or phenotypes. STRING collects, scores and integrates a spectrum of publicly available sources (e.g., ENSEMBL, GeneCards, KEGG, NextProt, RefSeq and UniProt) of physical and functional protein-protein interaction information, coupling them with computational predictions, to offer operations that ease the formation of comprehensive protein/gene networks (86). Downstream analysis and interpretation of results was supported by the MCL (Markov Clustering) clustering algorithm which was applied on the STRING formed network of the identified DEGs. MCL is a fast and scalable unsupervised graph clustering algorithm (87). It has proven its superiority in extracting clusters from interaction networks (88), and its effectiveness for protein association network analysis (89). MCL is applied on the weighted network of induced DEGs, with the weight of connecting gene edges acquired from relevant information sources (text-mining, experimental, databases, co-expression, neighborhood, gene-fusion and co-occurrence). By retaining the highly confident connections (the highest cutoff of 0.9 was used) the formed robust and coherent clusters of gene networks offer informative hints about the key molecular fingerprints underlying SCOV2 infection, and serve the biological interpretation of our findings. Exploiting the rich sources of molecular knowledge managed by STRING to generate gene association networks, and in particular, our focus on the coherent clusters of networked genes as key molecular fingerprints underlying infection, represent a novel and at the same time, rational and well-designed analytical approach for corresponding efforts in the field. The STRING visualization services were used in order to picture the constructed gene network clusters which, were then manually adjusted to be included in the reported figures.

The devise of diagnostic and prognostic classifiers is based on the union over the induced fingerprint genes that differentiate and characterize: the early SCOV2 infection stage, the low viral-load molecular profile of infected individuals, as well as the profile of early responders to the infection. This union consists of 52 unique genes that were used as descriptors for the devise of the respective SCOV2 diagnostic and COVID-19 prognostic classifier models. The Weka open-source machine-learning environment (www.cs.waikato.ac.nz/ml/weka/) was utilized for the devise and performance assessment of the classifier models (details are presented in see section 5).

The designation of infection stages and the assignment of samples into respective classes was performed by following a k-means/PCA-plot process on the most variable genes (for details refer to section 2.1.2/Analysis pipeline/′Designation of infection stages and phenotypes′). (Figure 3B) shows the Elbow-plot where, two clusters found to be the optimal choice for k. From the two induced clusters, one (consisting of 1,199 gene-transcripts) found to be largely dominated by enriched GO-biological processes that are directly related to cell-mediated host immune/defense responses, including: ‘Response to type I interferon’, ‘Defense response to virus’, ‘Cellular response to type I interferon’, ‘Negative regulation of viral process’, ‘Response to virus’, ‘Interferon-gamma-mediated signaling pathway’, ‘Negative regulation of viral genome replication’ and ‘Regulation of viral life cycle’. (Figure 3C) shows the PCA plot of samples on the expression-profiles of the 1,199 cluster genes. With a careful inspection of the plot, and by taking in consideration that both PC1 and PC2 components are strongly correlated with the hpi variable (p=1.24e-14 and p=7.71e-12 for PC1 and PC2, respectively), two separable groups of samples could be clearly identified − one with samples in hpi time-points 0, 1 and 2 (n=9), and one with samples in hpi 12 (n=3). These hpis designate the respective EARLY and LATE SCOV2 infection stages, with the samples in each group to be assigned to the corresponding stage class (a total of n=12 samples). As the scope of our study is to intensely contrast between the extreme stages of the infection course, the samples at intermedia hpis 3 and 6 are discarded. The clear separation between the two infection stage classes (when the intermedia hpis are removed) could be observed by inspecting the PCA plot of the retained 12 samples on all 9,457 genes (Figure 3D).

The gene-expression profiles of the retained SCOV2 EARLY/LATE samples were subjected to DEG analysis. Note that the analysis is performed on all, retained after filtering, 9,457 genes, and not just on the 1,199 genes of the enriched cluster used to designate the infection stages. By setting FDR≤0.05 and FC≥2, a total of 132 gene-transcripts (DEGs) were found to significantly differentiate between the LATE and EARLY stages; 79 up- and 53 down-regulated in the EARLY infection stage. All induced DEGs induced by the experiment with the GSE151513 dataset (and for all other experiments reported in the sequel) are included in the DEGs.xlsx file that is deposited in the Supplement Material accompanying the paper (consult sheet ′1_GSE151513_SCOV2_2-stage′). Here it is crucial to note that the differentiation polarity is inverted, that is, down-regulated DEGs in the EARLY stage are up-regulated in the LATE stage, and vice-versa. In an effort to gain insight into the functional roles and the interactions between the 53 down-regulated DEGs we followed a clustering approach using the MCL graph-based clustering algorithm (refer to section 2.1.2/Analysis pipeline/′Designation of infection stages and phenotypes′). Keeping only high-confidence gene interactions (i.e., a cutoff of 0.9 was used in the construction of the STRING gene network), 23 DEGs, all protein coding genes, were coherently grouped into three interconnected clusters (Figure 4A); the three clusters are indicated with different colors). We consider these genes as representatives for the two-stage SCOV2 core molecular fingerprint (for the GES151513 dataset). In addition, as shown in (Figure 4B), these genes exhibit low and high expression levels in the EARLY and LATE infection stages, respectively; the colored-shaded circles at the left of the figure indicate the mean expression value⁶ of the 23 genes over the samples assigned to the respective infection stage (indicated with cyan and dark-pink shaded circles for the SCOV2_EARLY and SCOV2_LATE samples, respectively), including the mean gene-expression value of all genes over all samples (indicated with black shaded circle).

What is important is that all the 23 genes are directly related to IFN-I signaling and especially to ISGs. Table 2 summarizes the functional roles and key antiviral activities of these genes. The findings agree with the discussion made in section 1/′The molecular canvas of immune and defense response during the infection course: the SCOV2 case′ about the special antiviral role of critical ISGs during the whole viral life-cycle process.

In order to validate our findings we collected indicative sets of genes with documented antiviral activity from four relevant studies, namely: (i) a set of 42 ISGs reported in a review about the antiviral functions of ISGs (90); (ii) a set of 33 IFN-I related genes reported in a recent study from the COvid-19 Multi-omics Blood ATlas (COMBAT) consortium (www.combat.ox.ac.uk) (91); (iii) a set of 25 ISGs reported in a recent study to contrast the gene expression profiles of SCOV2 infected samples between adults and children (92); and (iv) a set of 24 key genes interfering the viral life cycle reported in (38). The union of these four gene sets comprises 97 unique genes (consult the provided Supplement File DEGs.xlsx/′1_GSE151513_SCOV2_2-stage′). The fact that 21 of the 23 (91.3%) genes included in the aforementioned two-stage SCOV2 core molecular fingerprint fall within this union supports our findings and provide a strong evidence for the impaired immune/defense host response taking place during the early SCOV2 infection stages.

Here we have to make a special note about three genes, MT1F, MT1G and MT2A which, although not included in the 23 fingerprint genes, they are included in the list of induced 53 induced down-regulated DEGs in the early SCOV2 infection stage (notice also the down-regulation of the enriched ′Metallothioneins bind metals′ enriched REACTOME pathway induced by the enrichment/pathway analysis presented in the sequel; (Figure 4C). These genes belong to the family of Metallothioneins (MTs), a family of small highly conserved cysteine-rich metal-binding proteins. MTs regulate zinc (Zn) (93) that has a beneficial role in physiological and molecular host defense mechanisms during various pathogen infections, including SCOV2 (94). In-vitro experiments on mice reveal a direct and strong increase in the mRNA levels of MTs during acute influenza/A infection, especially at the upper respiratory tract (95). The physiology underlying this increase is attributed to the beneficial antioxidant role of MTs as they are triggered in order to effectively ‘clean-up’ the reactive oxygen species (ROS) generated by the host defense phagocytes during infection. It is also known than Zn contributes to host defense responses by maintaining the membrane barrier structure and function (96) via the modulation of cytokine-induced epithelial cell barrier absorptiveness (97). At the molecular level, recent studies demonstrate that Zn is required for interferon-mediated expression of MTs (98), and helps to enhance IFN-I response during SCOV2 infection, exhibiting an inhibitory ability of SCOV2 RNA polymerase through its failure to be associated with serious clinical outcomes (99). Furthermore, Zn deficiencies are directly linked to anosmia and taste dysfunctions (ageusia), already established as common SCOV2 symptoms (100), especially when decreased levels occurs in the nasopharyngeal tract (101). In addition, it is evidenced that acute viral infection of the nasopharyngeal mucosa lead to a decrease in local Zn levels as part of the normal defense against respiratory pathogens (102).

In order to reveal and highlight the molecular events that take place during the progress of SCOV2 infection, we proceed to the identification of enriched biological processes, pathways and hallmark gene signatures that contrast between EARLY and LATE SCOV2 infection stages. Based on the induced 132 DEGs, a number of GO-biological processes, REACTOME/KEGG pathways and HALLMARK_MSigDB gene signatures were found as significantly enriched and down-regulated in the early infection stage (Figure 4C). The results signify the fact that during the early stage of the infection, key biological processes and pathways engaged in first-line host innate immune/defense responses are ‘blocked’, with the list largely dominated by IFN/cytokine signaling and pro-inflammatory processes. Notably, even the KEGG COVID-19 pathway is found to be down-regulated during the early infection stage. The findings provide strong evidence that down-regulation of key antiviral immune/defense processes during the early course of SCOV2 infection and their up-regulation at the later stages may be the cause for the later emergence of uncontrolled, exaggerated and acute inflammatory COVID-19 clinical outcomes.

We performed the same as above analysis using only the mock-treated samples, with the respective samples assigned to the same infection stages, i.e., samples in hpis 0,1,2 to the early, and samples in hpi 12 to the late stage, respectively. DEG analysis resulted into 119 DEGs; 81 up- and 38 down-regulated in the early stage (data not shown). None of these genes (either up- or down-regulated) belongs to the list of the identified 23 genes that consist the formed two-stage SCOV2 core molecular fingerprint. Performing enrichment analysis of these genes, a list of fifteen significantly enriched biological processes were identified, all up-regulated in the early stage. The list was strongly dominated by cell-cycle events, including ‘Regulation of programmed cell death’, ‘Regulation of apoptotic process’, ‘Apoptotic process’, ‘Cell death’, ‘Programmed cell death’, ‘Negative regulation of cell population proliferation’, and ‘Cellular developmental process’. This highlights the fact that, even if both infected and mock-treated samples are grouped according to the hpi variable (consult Figure 3A), their molecular backgrounds differ drastically. In other words, the infected samples follow the cell-cycle ′norm′ but with overwhelming immune/defense molecular events being under- or over-expressed during the early or late infection stages, respectively. This proves the soundness of our approach to focus only on the SCOV2 samples and contrast between the designated early/late stages in order to explore and assess the two-stage molecular profile that governs the progression of the infection.

Confirmation of the two-stage SCOV2 infection profile was done using the GSE158930 RNAseq dataset (see section 2.1.1/′Task1′). Following the same methodology to identify infection stages, samples in hpi 4 and 24 (n=4) were assigned to the EARLY, and samples in hpi 72 and 96 (n=4) to the LATE SCOV2 infection stage, respectively (see Supplementary Figures 2A–C) for details). Setting FDR≤0.05, FC≥4 (a higher FC threshold was set in order to control the number of differential genes and ease their interpretation), and setting the ′rep′ variable as an extra factor to take in consideration the ′rep′licated/paired samples, a set of 175 DEGs (all protein coding genes) were found to differentiate between the LATE and EARLY stages (consult Supplement File DEGs.xlsx/′1_GSE158930_SCOV2_ 2-stage′); 15 up- and 160 down-regulated in the EARLY infection stage. Applying MCL clustering on the down-regulated genes and following the same network analysis methodology as in the previous experiment, 33 DEGs (all protein coding genes) were coherently grouped into two interconnected clusters (the two-stage SCOV2 core molecular fingerprint for GSE158930, (Figure 5A). Their expression levels are shown in (Figure 5B), and the respective enriched entries (biological processes, pathways and hallmark MSigDB signatures) in (Figure 5C).

The union of the two core two-stage SCOV2 molecular fingerprints identified from the previous experiments (23 genes for GSE151513 and 33 genes for GSE158930) consist of 38 genes being down-regulated in the EARLY SCOV2 infection stage, with 18 of them to be shared between the two fingerprints, namely: IFI27, IFI6 (IFI27-like/ISG12), IFIH1, IFIT1/2/3/5, IFITM1, ISG15, MX1/2, OAS1/2/3, RSAD2, STAT2, UBE2L6 and XAF1 (consult Supplement File DEGs.xlsx/′1_finger.Genes_UNION_COMMON′). The 18 genes and their MCL clustering are shown in Supplementary Figure 3A). As for the respective enriched entries (Supplementary Figure 3B), a set of 12 down-regulated enriched biological processes are shared between the two experiments (80.0% for both experiments); 10 REACTOME pathways (66.7% and 90.0% for GSE151513 and GSE158930, respectively); 7 KEGG pathways (63.6% and 70.0% for GSE151513 and GSE158930, respectively); and 5 hallmark MSigDB signatures (83.3% and 62.5% for GSE151513 and GSE158930, respectively). It can be easily checked that the shared enriched entries are heavily dominated by innate immune/defense processes being down-regulated during the early SCOV2 infection stage.

A specific reference should be made here concerning four genes, namely: HLA-A, HLA-B, HLA-C and HLA-E. Even if these genes are not in the list of 33 genes that compose the core SCOV2 molecular fingerprint for GSE158930, they are in the list of the induced 160 down-regulated genes in the EARLY SCOV2 infection stage. All these genes are harbored at the MHC/HLA (major histocompatibility complex/human leucocyte antigen) genome region, known as the most human polymorphic gene region. It is established that the polymorphic abundance of MCH/HLA genes is the key for their critical role in the regulation of host immune responses to most attaching pathogens (103), including SCOV2 (104). The human MHC/HLA genes are engaged to antigen processing, presentation and immune modulation. Their main functional role is to present antigens to CD8+ cytotoxic T lymphocytes (105), whereas their expression is induced by IFN genes (106). Furthermore, a number of variants in HLA alleles/haplotypes are known to be associated with susceptibility and progression of various infections, including SCOV2, with some of the most indicative summarized in (107), namely: HLA-A*11:01/24:02 with protective role against SCOV2 infection; HLA‐A*24:02 and HLA-B*22 with SCOV2 infection susceptibility; HLA-A*25:01 with moderate disease outcomes; and HLA-A*01:01/02:01, HLA-B*15:03/27:07, HLA-C*05 and HLA-E*01:01 with severe outcomes.

In their original publications the providers of GSE151513 (57) and GSE158930 (58) datasets do not directly contrast between the early and late infection stages, as we done in our analyses. In any case, in the publication of the GSE151513 dataset, the authors report a set of 19 ISGs being up-regulated in the late SCOV2 infection stage. 68.4% (13/19) of them are shared with the set of 18 genes being common between the two SCOV2 core molecular fingerprints identified in the previous experiments. Among other, two genes, STAT2 and UBE2L6, with key antiviral activities (consult Table 2) in the list of the 18 genes are not in the list of the 19 reported genes. As for GSE158930, in their respective original publication the authors report a set of 26 genes being up- and down-regulated in the late (72, 96 hpis) and early (24, 48 hpis) infection stages, respectively, with 42.3% (11/26) of them to be shared with the aforementioned set of 18 genes. Among other, four genes with key antiviral activities in the list of 18 genes are not included in the list of the 26 reported genes, namely: ISG15, RSAD2, STAT2 and UBE2L6. The above indicate, on one hand the competence of our analysis methodology as founded on the designation and direct contrast between early/late infection stages, and on the other, the adequacy of our findings compared to published results.

Our results provide evidence and demonstrate a two-stage core molecular profile that governs the progression of SCOV2 infection from the early to the later stages. The profile is realized by an impaired response of key immune/defense processes during the early infection stage, as materialized by the inhibition of IFN-signaling and the ′blocking′ of key ISGs, with the inverse to hold during the late stages of the infection. The findings are consistent with, and support the results reported in other relevant studies in which, the direct correlation of prolonged activation of high IFN-I levels with disease severity is highlighted (108), and the strong early IFN-I/ISG response to be beneficial for the clinical outcome of the infection (109).