Standardized Echinacea purpurea extract primes an IFN specific innate immune monocyte response via JAK1-STAT1 dependent gene expression and epigenetic control of endogenous retroviral sequences

Background: Herbal remedies of Echinacea purpurea tinctures are widely used today to reduce common cold respiratory tract infections. Methods: A system biology approach involving genomewide transcriptome, epigenome and kinome activity profiles was applied to characterize the immunomodulatory effects of a standardized Echinacea purpurea extract Echinaforce® in THP1 monocytes. Results: Gene expression and DNA methylation analysis revealed activation of innate immunity related antiviral signaling networks, involving IFN, chemotaxis and immunometabolic pathways. Furthermore, phosphopeptide based kinome activity profiling and pharmacological inhibitor experiments with filgotinib confirm a key role for Janus Kinase (JAK)-1 dependent gene expression changes in innate immune signaling. Finally, Echinaforce® treatment induces DNA hypermethylation at intergenic CpG, long/short interspersed nuclear DNA repeat elements (LINE, SINE) or long termininal DNA repeats (LTR). This changes transcription of flanking endogenous retroviral sequences (HERVs), involved in an evolutionary conserved (epi)genomic protective response against viral infections. Conclusions: Altogether, our results suggest that Echinaforce® phytochemicals strengthen antiviral innate immunity by JAK1 responsive gene expression and epigenetic regulation of HERVs in monocytes, which supports its prophylactic use against common cold corona viruses (CoV), Severe Acute Respiratory Syndrome (SARS)-CoV, and new occurring strains such as SARS-CoV-2. . Total RNA was isolated using the RNeasy mini kit (Qiagen, Hilden, Germany) including a DNAse treatment step as suggested by the manufacturer. Then 750 ng RNA was reverse transcribed into cDNA using oligo dT (Invitrogen), M-MLV reverse transcriptase (Promega, Wisconsin USA), 2.5 mM dNTPs and RNaseOUT (Invitrogen). Samples were incubated on 42°C for 60 min and 75°C for 15 min. For the HERV genes, cDNA synthesis was performed using random primers (Invitrogen) and incubation of the samples at 37°C for 60 min and 75°C for 15 min. qPCR was performed using the GoTaq qPCR Master Mix (Promega, Wisconsin USA) on a StepOnePlus Real-Time PCR machine (Applied Biosystems). Following primers were used: MX1 forward primer 5’-GTTTCCGAAGTGGACATCGCA-3’, MX1 reverse primer 5’-CTGCACAGGTTGTTCTCAGC-3’ (NM_001144925), IFITM1 forward primer 5’-CCAAGGTCCACCGTGATTAAC-3’, IFITM1 reverse primer 5’-ACCAGTTCAAGAAGAGGGTGTT-3’ (NM_003641), STAT1 forward primer 5’- CCATCCTTTGGTACAACATGC-3’, STAT1 reverse primer 5’-TGCACATGGTGGAGTCAGG-3’ (NM_007315), IL8 forward primer 5’-GCTCTCTTGGCAGCCTTCCTGA-3’, IL8 reverse primer 5’-ACAATAATTTCTGTGTTGGCGC-3’ (NM_000584), CXCL10 forward primer 5’-GAAAGCAGTTAGCAAGGAAAGGT-3’, CXLC10 reverse primer 5’-GACATATACTCCATGTAGGGAAGTGA-3’ (NM_001565), ACTB forward primer 5’-CTGGAACGGTGAAGGTGACA-3’, and ACTB reverse primer

the safety and efficacy of Echinaforce® to prevent common cold symptoms, showed significantly less cold episodes and of shorter duration as well as lower infection recurrence rate in the Echinaforce® treated versus placebo treated group [24]. Moreover, no differences between placebo and Echinaforce® group were reported in relation to health risk and safety [24]. Despite the promising immune potentiating properties of Echinaforce®, the responsible molecular targets have only partially been identified. For example, multiple studies demonstrate that Echinacea alkylamides exert their action partially through the cannabinoid receptor 2 (CB2) [8,20,25]. Alkylamides are structurally similar to endocannabinoids, and bind CB2 with a higher affinity compared to endogenous cannabinoids [25]. This action of alkylamides on CB2 was the mechanism behind the upregulation of TNF-α in primary monocytes after Echinaforce® treatment [20]. In addition, the cAMP, p38/MAPK and JNK signaling pathways, as well as NF-κB and ATF2/CREB1 transcription factors were found to be involved in the Echinacea-induced TNF-α expression [20].
Furthermore, a recent study showed that ex vivo blood stimulation with LPS after an 8-day oral Echinacea administration resulted in the induction of anti-inflammatory cytokines (IFN-γ, IL-8, IL-10 and MCP-1) only in subjects with low basal levels of these cytokines [11]. Interestingly, individuals with high basal levels of these cytokines didn't show any further induction. The same was seen in subjects with high stress levels or high susceptibility to cold infections, suggesting that Echinaforce® mainly (or only) enhances immune responses in immunocompromised subjects [11]. To further clarify its mode of action, we applied a system biology approach by integrating genomewide transcriptome, epigenome and kinome activity profiles of THP1 monocytes treated with Echinaforce®.

Cell lines and treatments
Echinaforce® (batch nr. 040070, A. Vogel Bioforce AG, Roggwil, Switzerland) is a standardized preparation obtained by ethanol extraction of freshly harvested Echinacea purpurea herb and roots (95:5). Echinaforce® is marketed as a registered medicinal product and produced under conditions of Good Manufacturing Practice (GMP). Thus, a consistent quality for each produced batch is mandatory and equal to the requirements for an allopathic remedy. The composition of marker compounds like alkylamides (i.e. those compounds known to characterize this species of Echinacea) was described previously [3,7,26]. In contrast to pressed juice extracts, Echinaforce® extract does not contain polysaccharides which are known to stimulate the immune system nonspecifically [27][28][29][30]. The alcohol concentration of Echinaforce® tincture extract was 65% v/v and solvent controls have been included in all experimental in vitro experiments to rule out nonspecific effects. In addition, the preparation was free of detectable endotoxin as determined by means of a commercial assay kit with a lower limit of detection 0.1 unit/ml (Lonza Walkersville Inc., MD). THP1 cells were grown in RPMI-1640 medium supplemented with glutamine, 10 % heat inactivated Fetal Bovine Serum, 50 IU/mL Penicillin, 50 µg/mL Streptomycin, 10 mM HEPES and 0.05 mM β-mercaptoethanol. Cells were treated with 1% Echinaforce® tincture versus ethanol solvent control. Each treatment condition consisted of six biological replicates.

Genome-wide gene expression analysis
Sample preparation and microarray processing THP1 cells were treated for 48h with 1% Echinaforce® or ethanol solvent control. RNA was isolated using the RNeasy mini kit (Qiagen) according to manufacturer's instructions. RNA concentration and purity was measured using the Nanodrop 1000 spectrophotometer (ThermoFischer, CA, USA). RNA integrity of each sample was checked using using the Experion Automated Electrophoresis System (Bio-Rad, MO, USA). Total RNA (500ng) was amplified using the Illumina TotalPrep RNA Amplification kit (Life Technologies, Carlsbad, CA, USA). Briefly, RNA was reverse transcribed using T7 oligo(dT) primers, after which biotinylated complementary or anti-sense RNA (cRNA) was synthesized through an in vitro transcription reaction. Then, 750 ng of amplified cRNA was hybridized to a HumanHT12 beadchip array (Illumina, San Diego, CA, USA) and further incubated for 18 hours at 58 °C in a hybridization oven under continuous rocking. After several consecutive washing steps, bead intensities were read on an Illumina iScan. Microarray data and raw gene expression intensities were preprocessed and analyzed using the beadarray R package [31]. Intensities were quantile normalized and log 2 transformed. Raw and normalized array data were uploaded to the Gene Expression Omnibus (GEO) database and have accession number: GSE117904. Probes with a P-detection value higher than 0.05 in at least six samples were removed. Also, probes annotated as "bad" and "no match" as described before [32] were not kept for further analysis. Differentially gene expression was performed using the limma R package [33]. P-values were corrected for multiple testing using the method of Benjamini and Hochberg. Probes with a log2 fold change higher than 0.4 and an adjusted p-value less than 0.05 were defined as significant and kept for further analysis [34]. The probes were annotated with gene information using the illuminaHumanv4.db annotation dataset [35]. The gene IDs of the significant Illumina expression probes were uploaded into the IPA software to find enriched biological pathways, diseases and networks. Metascape protein-protein-interaction enrichment analysis of differentially expressed genes was performed with the following databases: BioGrid, InWeb_IM, OmniPath [36]. The Molecular Complex Detection (MCODE) algorithm in Metascape identified densely connected network components and applies pathway and process enrichment analysis to each MCODE component independently. The three best-scoring terms by p-value are retained as the functional description of the corresponding MCODE components (as shown in the tables).

Kinase activity profiling
Sample preparation THP1 cells were treated with 1% Echinaforce® or ethanol solvent control for 15 min. Cell lysates were prepared according to manufacturer's instructions. In short, cells were washed twice with cold 1X PBS and lysed with lysis buffer (1:100 dilution of Halt Phosphatase Inhibitor Cocktail and Halt Protease Inhibitor Cocktail EDTA free in M-PER Mammalian Extraction Buffer (ThermoFisher Scientific™, Rockford, USA) at a ratio of 100 µl buffer per 1x10 6 cells. Lysates were then incubated on ice for 15 min and centrifuged for 15 minutes at 16000 x g at 4°C. Protein concentration was quantified using the Pierce BCA Protein Assay Kit (ThermoFisher Scientific™, Rockford, USA).

Serine/threonine kinases (STK) and tyrosine kinase (PTK) pamgene assay and data analysis
Kinase activity profiling was performed PamChip® preprocessing and kinase activity profiling was performed according to manufacturer's instructions (PamGene International BV, 's-Hertogenbosch, The Netherlands). The first part of the protocol consisted in the blocking of the arrays with 2% BSA followed by several washing steps. Then 0.5 µg for STK and 5 µg for PTK assays together with the correspondent reaction mixes (purchased from the Pamgene) were loaded onto the arrays and incubated in the microarray system PamStation ® 12 instrument (PamGene International, Den Bosch, The Netherlands). In this step, the ATP contained in the mix leads to the activation of the kinases in the lysate which will result in the phosphorylation of the peptides on the array. Peptide phosphorylation intensities are then detected with the primary STK antibody mix and FITC-labeled antibody for STK assay and with the FITC-labelled PTK antibody (PTK assay). Images are then taken by the CCD camera in the PamStation®12 and processed by the Bionavigator software. Peptide intensities data were log 2 transformed and differences in phosphorylation between Echinaforce® treated and control cultures were determined by using an univariate student t-test analysis corrected for multiple testing using the Benjamini and Hochberg method [34].
To identify potentially activated or inhibited kinases we used the STK or PTK Upstream Kinase analysis PamApp from the Bionavigator Software. The analysis is based on "in silico predictions" for the upstream kinases of phosphorylation sites in the human proteome that are retrieved from the phosphoNET database [38]. In short, a prediction algorithm is derived from known interactions between kinases and phosphorylation sites. The prediction algorithm is then used to predict the strength of undocumented interactions. The Bionavigator application uses PhosphoNet database to map putative kinases upstream of the phospho-peptides (a kinase can have multiple possible phosphosites, and a single site can be phosphorylated by different kinases). For each set of peptides mapped to a specific kinase, a "difference statistics" is calculated (=normalized kinase statistics) using following formula: as the sample mean and variance of the intensity of peptide i in group j, respectively, whereas n is the number of peptides linked with a specific kinase. A positive kinase statistic means that the kinase is activated, while a negative statistic means the kinase is inactivated compared to the control group. The kinases are subsequently ranked based on a specificity and significance score which are calculated using permutation of the peptides and samples, respectively. Following formula is used: where m is the number of times out of M permutations that |τ p | > |τ|, where τ p is the value of the difference statistic obtained after permutation of the samples or peptides. The significance score represents the magnitude of the change represented by the normalized kinase statistic. The specificity score represents the specificity of the of normalized kinase statistic in terms of the set of peptides used for the corresponding kinase. The higher the score the less likely it is that the observed normalized kinase statistics could have been obtained using a random set of peptides from the data set. The sum of the significance and specificity score is used to rank the kinases [39].

Genome-wide DNA methylation analysis
Sample preparation THP1 cells were cultured for 48h with 1% Echinaforce® or ethanol solvent control. Corresponding cellular genomic DNA was isolated using the DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) according to manufacturer's instructions. DNA concentration and purity was measured using the Nanodrop 100 spectrophotomer and 1 µg of DNA was used for bisulfite conversion using the EZ DNA methylation Kit of Zymo Research according to manufacturer's instructions. Successful bisulfite conversion was checked using a methylation-specific PCR in a region of the SALL3 gene (see [40] for primer sequences).

EPIC DNA methylation array
The Infinium HumanMethylationEPIC BeadChip array (Illumina, San Diego, CA, USA) was used to measure genome-wide DNA methylation. 4 µL of bisulfite-converted DNA from each sample was amplified, fragmented, precipitated, resuspended and subsequently hybridized onto the BeadChips. After overnight incubation of the BeadChips, unhybridized fragments were washed away, while hybridized fragments were extended using fluorescent nucleotide bases. Finally, the BeadChips were scanned using the Illumina iScan system to obtain raw methylation intensities of each probe.

EPIC DNA methylation data preprocessing and analysis
The R package RnBeads was used to preprocess the Illumina 450K methylation data [41]. CpG-probes were filtered before normalization based on following criteria: probes containing a SNP within 3 bp of the analyzed CpG site, bad quality probes based on an iterative greedycut algorithm with a detection p-value threshold of 0.01, and probes with missing values in at least one sample. After filtering these CpG-probes, methylation values were within-array normalized using the beta mixture quantile dilation (BMIQ) method [42]. Another filtering step was performed after normalization based on following criteria: probes measuring methylation not at CpG sites (CC, CAG, CAH, …) and probes on sex chromosomes.
The methylation beta-values were transformed to M-values (M = log 2 (β/(1-β))) prior to further analyses. The moderated t-test incorporated in the limma R package [33] was used to calculate the statistics and p-values of the methylation differences between Echinaforce®-and solvent-treated samples. Significant differentially methylated probes (DMPs) were selected based on a false discovery rate (FDR) < 0.1 and a difference in betavalue of at least 0.05. The DMPs were annotated with gene information using the IlluminaHumanMethylationEPICmanifest R package [43]. Further gene information was retrieved from the UCSC genome browser (human hg19

Assessment of IL-8 and CXCL10 levels
Cell culture supernatants were collected after 3, 6, 12, 24 and 48 hours and assayed for chemokines CXCL10 and IL8 by means of an enzyme-linked immunosorbent assay (ELISA) purchased from Invitrogen (CA, USA) following manufacturer's instructions. The assays have a detection limit of 2pg/ml for CXCL10 and 5 pg/mL for IL-8.

Echinaforce® treatment activates innate immunity genes and pathways
Widespread gene expression changes in monocyte THP1 cells were detected upon 48h 1% Echinaforce® treatment. Based on significance criteria of FDR < 0.05 and absolute log 2 fold change > 0.4, Echinaforce® induced modest upregulation of 205 expression probes (173 genes) while 124 probes (99 genes) were downregulated compared with the ethanol treated solvent controls ( Figure 1A and Supplementary table 1).
Complementary to IPA analysis, protein-protein-interaction enrichment analysis of DEGs by STRING [53] and Metascape [36] algorithms was performed. This revealed strong enrichment of protein-protein interactions responding to a chemical stimulus, which triggers a defensive antiviral innate immune response involving IFN, TLR, NOD, RIG, cytokine, chemokine and NF-κB signaling pathways (Figure 2, Supplementary Table 4). More particularly, Metascape MCODE analysis identified 3 interconnected subnetworks in the antiviral cytokine response: cellular response to interferon, regulation of leukocyte chemotaxis and (mitochondrial) metabolism (Supplementary table 4).
To evaluate and validate the time dependent transcriptional dynamics of different genes involved in IFN and chemotaxis innate immune signaling observed in the microarray data, mRNA expression and protein levels for STAT1, MX1, IFITM1, CXCL8 and CXCL10 were measured using qPCR, ELISA and immunoblotting assays in THP1 monocytes following Echinaforce® treatment. Induction of STAT1 and the interferon-stimulated genes MX1 and IFITM1 expression could clearly be confirmed, with maximal mRNA transcription levels observed after 48h treatment ( Figure 3A). Corresponding changes in STAT1 protein expression levels could also be detected by Western analysis (Figure 3B), whereas we failed to detect significant amounts of MX1 and IFITM1 protein (data not shown). Whether MX1 and IFITM1 protein expression is susceptible to additional posttranscriptional silencing mechanisms (miRNA, lncRNA) or posttranslational modifications which result in high turnover rates and low protein expression levels needs further investigation [54,55]. For the chemokines IL8 and CXCL10, persistent gene induction could be observed until 48h, with peak transcription levels after 3h (Figure 3A). Accumulation of these chemokines in the cell culture supernatants was also detected using ELISA technique and showed that the levels of both chemokines significantly increased with the incubation time of Echinaforce®, as compared to control setups.

Echinaforce® treatment activates JAK1, NF-κB and MAPK kinases
To identify most important upstream kinase pathways responsible for gene expression changes in THP1 monocytes following Echinaforce® treatment, we performed a Pamchip kinome activity profiling assay [39]. This peptide array approach allows characterization of cellular serine/threonine or tyrosine kinome activity profiles following on chip in vitro kinase reaction of 144 conserved kinase consensus peptide motifs in presence of THP1 monocyte lysates left untreated or following Echinaforce® treatment [56][57][58][59]. Using the upstream kinase prediction tool of the Bionavigator PamGene software, the qualitative and quantitative changes in phosphopeptide chip intensities upon Echinaforce® treatment were translated into a pattern of activated or inhibited upstream kinases (Figure 4A and Supplementary table 5). In agreement with the transcriptional activation of the IFN signaling pathway described above (Figure 1C), Pamchip kinome profiling [39] revealed activation of the JAK1 kinase which is important in the phosphorylation of STAT kinases and subsequently downstream regulation of IFN-stimulated genes. Furthermore, in line with pathway analysis of transcriptome data, we also identified activation of the tyrosine kinase TEC (Figure 4B) (Supplementary table 2-4). Surprisingly, our analysis did not detect significant activity changes of early IFN kinases TBK1 and IKK [60].
Besides, we also identified various Echinaforce® activated kinases belonging to the MAPK superfamily of kinases: p38 MAPK (MAPK11, -12, -13, and -14), JNK (MAPK8, -9 and -10) and ERK1 (Figure 4C). This upstream regulators are also predicted by IPA to control various canonical pathways, including pattern recognition receptors in recognition of bacteria and viruses, activation of IRF by cytosolic pattern recognition receptors and role of MAPK signaling in the pathogenesis of influenza among others.
To further verify crucial involvement of JAK kinase activation in downstream gene expression effects upon Echinaforce® treatment, we compared THP1 gene expression changes following Echinaforce® treatment in presence or absence of the pharmacological JAK1 inhibitor filgotinib. We found that filgotinib significantly suppresses the Echinaforce® responsive genes MX1 and IFITM1, whereas STAT1, CXCL10 and IL8 gene expression were less significantly suppressed ( Figure 4D). Altogether, experiments with the JAK1 inhibitor filgotinib strenghten our transcriptome and kinome data analysis, pointing to JAK1-specific regulation of downstream gene expression changes in response to Echinaforce® treatment.

Echinaforce® treatment elicits epigenetic changes in innate immunity gene pathways
Epigenetics seems to be important during monocyte differentiation and in the immunological memory of macrophages [61,62]. Today, various bioactive phytochemicals have been identified which modulate inflammation through epigenetic reprogramming [63,64]. Different phytochemicals and nutrients are known to change DNA methylation and histone modifications by directly influencing epigenetic enzymes or by interfering with the availability of the substrates/cofactors of these enzymes [65][66][67]. To assess whether the Echinaforce® induced changes in transcriptome profiles in THP1 cells are associated with DNA methylation changes, we measured complementary changes in DNA methylation profiles using the Illumina EPIC methylation array. Significant DNA methylation changes were observed following 48h exposure to Echinaforce® (Figure 5A and Supplementary table 6).
A total of 1,875 CpG sites was found differentially methylated (FDR < 0.1) with a methylation difference of at least 5%. Typically, DNA methylation changes after short (24-72h) exposure to phytochemicals and nutrients are much smaller than cancer associated DNA methylation changes in oncogenes or tumor suppressor genes which accumulate for many years in response to the microenvironment [46,68,69]. However, similar DMR effects sizes and cutoff (<5%) were found to be biologically meaningful in various disease etiologies [40,70,71].
From the 1,875 CpG sites identified, only 40 differentially methylated positions (DMPs) were hypomethylated whereas 1,835 DMPs were hypermethylated. DMPs were mainly enriched in gene bodies, intergenic, and CpGpoor regions, while depleted in CpG islands, promoter, and enhancer regions (Figure 5B). Only 1,259 of the 1,875 CpG-probes (67%) were located in a gene or 1,500 bp upstream of a gene. Similarly, DNA methylation variation in the immune system was predominantly found at at CpG islands within gene bodies, which have the properties of cell type-restricted promoters, but infrequently at annotated gene promoters or CGI flanking sequences (CGI "shores") [72]. Subsequent IPA pathway enrichment analysis of the genes containing DMPs revealed inflammation or immunological diseases among others (Supplementary Table 6). Of particular interest, one of the top enriched pathways ('Superpathway of Inositol Phosphate Compounds') controls various epigenetic processes related to the interferon response [73][74][75].
Since both gene expression and kinase profiling both revealed the involvement of interferon signaling pathways, we also checked whether methylation of IFN pathway genes was affected by Echinaforce® treatment. Eight probes located in BCL2, JAK1, STAT1, PIAS1 and TAP1 did show an FDR < 0.1, with small methylation differences (between 1 and 3%) (Figure 5C). Whether these small methylation changes are sufficient to finetune the immune gene response needs further investigation.
Since most of the DMPs were located in intergenic regions and gene bodies, only a small subset of genes containing a DMP also resulted in a significant change in gene expression (Figure 5D). Only seven genes were both differentially methylated and expressed, based on the significance criteria described above: i.e. Calsyntenin 2 (CLSTN2), Enhancer Of Zeste 2 Polycomb Repressive Complex 2 Subunit (EZH2), Growth arrest-specific protein (GAS)-7, neuron navigator (NAV)-3, Thioredoxin Reductase (TXNRD)-1, Tryptophanyl-tRNA synthetase (WARS) and Zinc Finger Transcription Factor (ZNF)-644. When using less stringent significance criteria, leaving out the effect size cutoff (logFC), 574 CpG site -gene pairs were found to be differentially expressed and methylated. Upon further comparing canonical pathways which are significantly enriched for both lists of differentially expressed genes and the list of differentially methylated genes, we identified 10 common biological processes ( Figure 5E). Remarkably, common pathways include NF-κB signaling (NF-κB activation by viruses, NF-κB signaling), MAPK signaling (LPS-stimulated MAPK signaling, UVA-induced MAPK signaling), and immune responses (i.e. Role of pattern recognition receptors in recognition of bacteria and viruses, Role of NFAT in regulation of the immune response, phagosome formation, CD40 signaling, leukocyte extravasation signaling).

Discussion
In this study, we applied for the first time a systems biology approach to characterize a possible mode of action of a standardized medicinal Echinacea purpurea tincture Echinaforce®, which is widely used as a herbal remedy against respiratory tract infections. Microarray, QPCR, Western and ELISA based assays demonstrate that treatment of THP1 monocyte cells with Echinaforce® phytochemicals elicit time dependent gene expression changes in innate immunity signaling networks, involving the IFN (MX1, IFITM1, STAT1, STAT2) chemotaxis (IL8, CXCL10) and immunometabolic (ISG15, PKM2, SQSTM1) pathway.
Most cells express a set of membrane and cytoplasmic receptors to detect viral RNA and DNA molecules: Pattern Recognition Receptors (PRRs). These receptors control innate immune signaling to activate the synthesis of interferons during a viral infection. In addition to pathogens, autophagy, metabolic and chemical stress, DNA damage, unfolded protein response, can also regulate innate immunity through cell-autonomous responses. Either IFN-inducible or constitutive, these processes aim to guarantee cell homeostasis or a biodefense mechanism against (non-self) hazardous molecules [81]. Of importance, these distinct constitutive cellautonomous responses appear to be interconnected and can also be modulated by microbes, viruses and PRRs [82]. Our results suggest that Echinaforce® phytochemicals prime innate immunity pathways via activation of interferon and chemokine gene expression. As such, secondary metabolite phytochemicals involved in plant immunity may train evolutionary conserved innate immune responses across species [83][84][85]. For example, transcriptional upregulation of the protein kinase receptor (PKR, EIF2AK2), a cytoplasmatic pattern-recognition receptor could be observed. PKR is known to transduce RNA helicase (MDA5) dependent virus signals for type I IFN induction [86]. Interferon regulatory factor 7 (IRF7) is another key protein found strongly upregulated. Transcription factors IRF7 together with IRF3 regulate expression of early type I IFN and other proteins involved in the innate antiviral immune response (activation of IRF by cytosolic pattern recognition receptors) [87] (Supplementary tables 1-2-3-4). Signal transduction via PKR occurs mainly via NF-κB and MAPK pathways (Role of PKR in Interferon induction and antiviral response) [88]. Another important intracellular pattern-recognition receptor for viral RNA which was found to be upregulated by Echinaforce® was the RNA helicase MDA5 (IFIH1) [60]. Furthermore, upregulation of the NF-κB subunits RelB and NFKB2/p52 was observed, which can promote downstream production of innate immunity chemokines (NF-κB activation by viruses, NF-κB signaling) [89].
In line with our results which show activation of innate immunity and inhibition of viral infection and replication gene responses, anti-viral effects against influenza infection and activation of IFN pathways have also been reported in vivo following Echinaforce® tincture treatment [11,15,21]. Our in vitro results are also in line with observations in human studies ex vivo/in vivo showing increased immunomodulating as well as chemotactic neutrophil effects following Echinaforce® treatment, especially in immunocompromised people [11,23,24,90]. For example, the antiviral ability of CXCL10 has been attributed to its chemoattractant effects which promote recruitment of neutrophils [91][92][93][94]. The latter illustrates that both neutrophils and inflammatory monocytes are intertwined in the immune system's anti-viral response [95]. Similar results were previously obtained in murine dendritic cells, illustrating that Echinaforce® stimulates cell mobility and chemotaxis and alters expression of cell adhesion and motility genes [96]. Other studies showed that Echinaforce® may reverse the chemokine induction of virus-infected cells [12,[97][98][99]. Paradoxically, Echinaforce® may induce cytokine and chemokine expression in uninfected cells, but suppress their expression upon virus infection or LPS stimulation [30,[97][98][99]. Similarly, Echinaforce® increased the transcription of TNF-α in human monocytes, but reduced the LPSstimulated TNF-α protein production [20]. Although studies suggest that this stimulatory effect may be the result of bacterial-derived LPS and lipoproteins [27][28][29][30], our Echinaforce® tincture contains no polysaccharides, neither endotoxins. Altogether, the latter suggests that its immunomodulatory effects are due to the active compounds present in the formulation [13,20]. Similar activation of IFN innate immunity and viral protection has been observed in presence of avocado and apple extract [100,101]. Interestingly, in the latter case, effects were attributed to oligomeric proanthocyanidins and lost with their monomeric form [101].
With respect to immunometabolism, mitochondrial metabolism shows a remarkable sensitivity to chemokine and IFN signaling [95,102]. For example, ISG15 is an interferon-stimulated, ubiquitin-like protein which regulates mitochondrial homeostasis and targets various proteins involved in catabolic autophagy metabolism in the mitochondria (mitophagy) during infection [103,104]. Moreover, mitochondrial changes in immunometabolism (glycolysis, the tricarboxylic acid (TCA) cycle, the pentose phosphate pathway, fatty acid oxidation, fatty acid synthesis and amino acid metabolism) strongly contribute in (re)shaping immunity and production of neutrophil extracellular traps (NETs) [105][106][107][108].
Finally, genomewide epigenetic analysis of DNA methylation changes following Echinaforce® treatment revealed almost 2000 DMP, enriched for immune disease and immunological pathways. Although the observed methylation changes are relatively small after 72h treatment, cumulative effects can contribute in building an immune memory response by priming chromatin to mount faster and higher innate immune transcription upon re-stimulation of immune cells [84]. Besides the regulation of gene expression, DNA methylation is also involved in regulating alternative splicing, intron retention or promote cryptic transcription of non-annotated TSSs (TINATs) encoding immunogenic peptides which might prime an antiviral innate immune response [124][125][126][127]. As such, it appears that Echinaforce® treatment predominantly promotes epigenetic changes in innate immunity gene pathways and to a less extent of adaptive immune response genes.
Besides, the higher global DNA hypermethylation observed after Echinaforce® treatment in LINE, SINE and LTR transposon repeats flanking endogenous retroviral sequences (HERVs), may be part evolutionary conserved (epi)genomic protective response against retrotransposition and viral infection [128,129]. Similarly, IFN was shown to promote DNA methylation silencing of repeats and noncoding RNAs [37,128,130,131]. Specific HERVs have been proposed to establish a protective effect against exogenous viral infections [132]. HERVs can act as IFN-inducible enhancers and have shaped the evolution of a transcriptional network underlying the IFN response [132][133][134][135]. Of particular interest, the MER41B family of ERV sequences contains a STAT1 binding site and regulates expression of IFN-γ-responsive genes, such as absent in melanoma 2 (AIM2), and IFI6 [136,137]. CRISPR-Cas9 deletion of a subset of these HERV elements in the human genome impaired expression of adjacent IFN-induced genes and revealed their involvement in the regulation of essential immune functions, including activation of the AIM2 inflammasome. Along the same line, DNA methylation inhibitors trigger an IFN response through viral mimicry via transcription of dsRNAs of repetitive elements from HERVs which can activate RIG-I and MDA5 PRRs. RNA transcripts of HERVs can be reverse transcribed to generate ssDNA or expressed to generate proteins with viral signatures, much like the pathogen-associated molecular patterns of exogenous viruses, which allows them to be detected by the innate immune system [138,139]. In another example, silencing of the MLT1C49 HERV decreased expression of CXCL10 and CCL2 chemokines [140]. Finally, transcriptional changes of MLT1B and MER4D HERV transcription and innate immune signaling have also been described upon immunometabolic mitochondrial changes in protein kinase (PK)-M2 activity, which were counteracted by NFkB RelB [141]. From these examples, it appears that HERV regulatory sequences now constitute a dynamic reservoir of IFN-inducible enhancers fueling genetic innovation in mammalian immune defenses [136,142,143].
Previous studies showed that Echinaforce®, besides its immunomodulating activities is also very active as a virucidal agent against viruses with membranes, i.e. HSV-1, respiratory syncytial virus, all tested human and avian strains of influenza A virus, as well as influenza B virus [144]. Along the same line, Echinacea polyphenol quercetin was found to inhibit SARS-CoV entry in the host through binding SARS-CoV-S2 protein [145]. Similarly, protective effects could recently also be observed in a reconstituted nasal epithelium cell culture system by exposing Echinaforce®-treated respiratory epithelium to droplets of HCoV-229E, SARS-or MERS-CoVs, imitating a natural infection [146]. In contrast Echinaforce® was found to be less effective against intracellular virus replication. Consequently, virus already present within a cell could be refractory to the inhibitory effect of Echinaforce®, but virus particles shed into the extracellular fluids would be vulnerable. Therefore the antiviral actions of the Echinaforce® may especially manifest during initial contact with the virus, i.e. at the inception of infection, and also during transmission of virus from infected cells.
In conclusion, our systems biology approach revealed that Echinaforce® phytochemicals trigger multiple antiviral immunemodulating bio-activities. Further studies in preclinical respiratory infection models and placebo controlled intervention studies are needed to identify discrete compounds or mixtures, which could be effective in prophylactic use against common cold corona viruses (CoV), Severe Acute Respiratory Syndrome (SARS)-CoV, and new occurring strains such as SARS-CoV-2.

Figure 1
Echinaforce® induced gene expression activates innate immunity pathways A) Volcano plot showing the upregulated genes (orange color, number of probes: 205), and downregulated genes (blue color, number of probes: 124) upon treatment of THP1 cells for 48h with Echinaforce® tincture (1%) . B) Top enriched IPA canonical pathways. Bars are colored by activation z-score. C) IPA interferon signaling pathway with Echinaforce®-induced upregulated genes colored in red and green, respectively. D) Top enriched IPA infectious diseases and IPA immune trafficking disease and biological function. Bar charts are colored by activation z-score.   Activation of JAK1 and MAPK kinases by Echinaforce®. A) Kinome activity profiling on THP1 cell lysates, following 15 min treatment with Echinaforce® tincture (1%). Showing predicted upstream kinases. Bars are colored by specificity score with red meaning the highest score. The direction of the bars represents the normalized kinase statistics. A positive kinase statistic means a higher activity in Echinaforce® treated samples. B) TEC signaling pathway as predicted by IPA software showing the up-and down-regulated genes (colored in red and green, respectively), after Echinaforce® treatment. C) IPA-enriched P38 MAPK and JNK pathways upstream regulators. Genes colored in orange are predicted to be activated, while genes colored in blue are predicted to be inhibited. D) Effect of JAK1 inhibition on transcript expression of interferon pathway related genes. THP1 cells were either treated with the JAK1 inhibitor Filgotinib alone or in combination with Echinaforce® (n=7). Mean expression LogFC change relative to solvent control is represented together with 95% confidence interval. *: P ≤ 0.05, ** P: ≤ 0.01, *** P: ≤ 0.001.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.