Low-dose Drosera rotundifolia induces gene expression changes in 16HBE human bronchial epithelial cells

Drosera rotundifolia has been traditionally used for the treatment of respiratory diseases in phytotherapy and homeopathy. The mechanisms of action recognized so far are linked to the known effects of specific components, such as flavonoids, but are not completely understood. In this study, the biological functions of D. rotundifolia were explored in vitro following the treatment of bronchial epithelial cells, which are the potential targets of the pharmacological effects of the herbal medicine. To do so, the whole plant ethanolic extract was 1000-fold diluted in water (D. rotundifolia 3×) and added to a 16HBE human cell line culture for 3 h or 6 h. The effects on gene expression of the treatments and corresponding controls were then investigated by RNA sequencing. The differentially expressed genes were validated through RT-qPCR, and the enriched biological functions involved in the effects of treatment were investigated. D. rotundifolia 3× did not impair cell viability and was shown to be a stimulant of cell functions by regulating the expression of dozens of genes after 3 h, and the effects were amplified after 6 h of treatment. The main differentially expressed genes encoded ligands of epithelial growth factor receptor, proteins involved in xenobiotic detoxification and cytokines, suggesting that D. rotundifolia 3× could stimulate self-repair systems, which are impaired in airway diseases. Furthermore, D. rotundifolia 3× acts on a complex and multifaceted set of genes and may potentially affect different layers of the bronchial mucosa.


Results
Analysis of polyphenols in the D. rotundifolia ethanolic extract. The concentration of polyphenols in the D. rotundifolia ethanolic extract was estimated to be 7.1 × 10 −3 mol/L in terms of the gallic acid equivalents. This concentration corresponds to approximately 1.2 mg/mL of total polyphenols, which in turn corresponds to a concentration of 12 µg/mL in the D. rotundifolia 3× test sample and 1.2 µg/mL (7.1 × 10 −6 mol/L) in the final cell culture. The D. rotundifolia ethanolic extract was also qualitatively analysed by mass spectrometry. In the HPLC 6.590 min retention time peak, the mass spectrum showed peaks corresponding to ellagic acid (m/z = 300.997905) and isoquercitrin (m/z = 463.087935). In addition, at the 8.004 min retention time peak, quercetin (m/z = 301.035282) was detected. Regarding the presence of hyperoside, which has the same molecular mass as isoquercitrin, there is a chance that the isoquercitrin peak represented a mixture of both compounds. The results confirmed the presence of important bioactive components in D. rotundifolia ethanolic extract 13 . Effect of D. rotundifolia 3× on 16HBE cell viability. The viability of the bronchial epithelial cells (16HBE) (Fig. 1) was not impaired after 3 h of exposure to D. rotundifolia 3× or the corresponding control (Ctrl), demonstrating the absence of cytotoxicity caused by the plant preparation at the dilutions used in the experiments. Furthermore, after 24 h of incubation with D. rotundifolia 3×, the cells showed a small (+ 6.9%) but significant (p = 0.019) increase in cell viability. No morphological changes or alteration of cell adherence to the bottom of the well could be observed by optical microscopic inspection (data not shown).

Changes in gene expression after D. rotundifolia 3× treatment. The effects of D. rotundifolia
3× treatment on the global gene expression of 16HBE cells were investigated after 3 h of incubation by comparison with the control. The experiments were reproduced in 4 different biological replicates, and the gene expression profile in response to D. rotundifolia 3× treatment was investigated by RNA-seq analysis.
The reads obtained by the sequencing were quantified as transcripts by Salmon, which showed an average mapping coverage of 92%. Differential gene expression analysis was then performed to identify the significant target genes of D. rotundifolia 3× treatment. The differential expression output obtained from DESeq2 analysis displayed a total of 69 statistically significant (padj < 0.05) differentially expressed genes (DEGs); 44 genes were upregulated (Table 1), and 25 genes were downregulated ( Table 2). The lof2FoldChange ranged from + 0.7 (maximum upregulation, Table 1) to − 0.51 (maximum downregulation, Table 2).
Gene ontology classification of differentially expressed genes. To better understand the functions of the identified genes, functional analysis was performed by submitting the up-and downregulated gene list to gene ontology (GO) enrichment analysis using the Bioconductor/R package gprofiler2. The representative enriched biological processes and the genes associated with these processes are represented in Table 3. The DEGs were mainly enriched in biological processes, including general functions such as "regulation of response to stimulus", which includes 34 of the 69 genes, and more specific biological activities, such as "vasculature development" (13 of 69 genes) and "positive regulation of vascular endothelial growth factor production" (5 of 69 genes). A group of 10 genes is associated with the function "epithelial cell proliferation". This gene set includes the upregulated gene LGR5 (a GPCR receptor protein and member of the Wnt signalling pathway that controls cell proliferation), the downregulated gene IGFBP3 (encoding insulin-like growth factor-binding protein 3, which exerts anti-proliferative effects in many cell types) and 8 other genes with specific GO functions related to epithelial regulation (the upregulated epithelial growth factor-like ligands EREG and EPGN, the downregulated EGFlinked inhibitor ERRFI1), proteins involved in cell membrane repair (DYSF) or in tissue development (NOG) and two transcription factors (HIF1A and JUN). A group of DEGs is related to the production of cytokines, suggesting their possible role in immunity. Other processes that are significantly associated with the target genes are "regulation of cell death" and "positive regulation of the mitotic cell cycle", suggesting their involvement in the recovery of epithelial tissue.
Time course of the D. rotundifolia 3× effect. Some representative and very well-expressed genes were selected to investigate the kinetics of the effects of D. rotundifolia 3× during the time course of the cell treatment.
A new series of four experiments was performed, in which the cells were incubated for 2 h, 6 h, and 24 h with the medicine or the control. Figure 2 displays the mRNA expression of AREG, CTGF, CXCL8, CYP1B1, EGR1, EREG, IL-1α and TIPARP in cells with and without D. rotundifolia 3× at the indicated times. The results of all four experiments are separately shown since one of the experiments resulted in higher levels of mRNA expression for a few genes (e.g., CXCL8, EREG and IL-1α). Despite the heterogeneity in the basal gene expression, the direction of the effect is clear-cut, and it is largely significant for at least one time point based on the statistics for the paired data. The experiments showed that the effect of D. rotundifolia 3× on most genes starts already at two hours, it is expressed in a higher degree after 6 h, while it decreases at 24 h. Exceptions to this trend are CYP1B1 and TIPARP, which showed the maximum activity at 2 h, which was decreased at 6 h and further decreased at 24 h.
Based on these results, a new RNA-seq analysis of the samples incubated for 6 h with D. rotundifolia 3× was performed to describe in detail and under the best conditions all of the effects of the plant on the gene transcription of bronchial cells.  using the same parameters used for the RNA-seq of 16HBE cells treated with D. rotundifolia 3× for 3 h as well as the bioinformatics analysis. All samples passed the quality control tests as described for the RNA-seq of the 3 h-treated samples. The DESeq2 output generated a list of 495 DEGs in 16HBE cells after 6 h of treatment with D. rotundifolia 3× compared to treatment with the control (padj < 0.05), which contains a majority of upregulated genes (n = 334). Figure 3 depicts the DEGs in a Volcano plot, which highlights (blue dots) all the genes (n = 117) that had a statistically significant change (adjusted p-value < 0.05) and a log2FoldChange higher than + 0.4 or less than − 0.4. For this analysis, a log2FoldChange cut-off value was used since the number of DEGs was much higher compared to that identified by the RNA-seq performed for the 3 h treatment. Figure 3 shows that after 6 h of incubation, there were a high number of DEGs, especially upregulated genes, which is in accordance with what is shown in Fig. 2. From a qualitative point of view, a similarity is observed in the genes whose expression is increased at the two different incubation times. Such a similarity is reflected by the upregulation of some genes, such as CYP1B1, EPGN, EREG, SERPINB2, IL-1α, AREG, and PHDLA1, which were also upregulated by cells treated with D. rotundifolia 3× for 3 h (see Table 1). Regarding the downregulated genes (left side of Fig. 3), some new genes associated with an inhibitory effect are highlighted, including IGFBP3, PLAC8, MAN2B2, and BNIP3, which, except for IGFBP3, were identified only after 6 h of incubation and not after a shorter exposure (2 or 3 h). Moreover, genes that were slightly downregulated after 3 h of D. rotundifolia 3× treatment, such as ERRFI1, were not found to be regulated after 6 h of treatment.
Functional analysis. The 117 differentially expressed genes were ranked by the adjusted p-value significance, and GO enrichment analysis was performed using gprofiler2 as described in the "Materials and methods" section. The enrichment results for the 6 h data set (Fig. 4) show a wide range of information concerning the biological processes (blue bars in Fig. 4) and molecular functions (red bars). Similar enriched biological processes were observed after 6 h and 3 h of D. rotundifolia 3× treatment, such as functions related to blood vessel development, regulation of epithelial cell proliferation and cytokine production. Moreover, there was an increase in biological processes related to inflammation and chemotaxis after 6 h of D. rotundifolia 3× treatment due to the To draw a functional picture of the enriched functions and their correlations, a network that associates the genes with the functions in the biological model of bronchial epithelial cells was constructed (Fig. 5). Functions are representative and were chosen based on our experimental model and enrichment p-value. The Log2Fold-Change values from the DESeq2 output are depicted by the coloured nodes, which shows that the enrichment of these biological functions is mostly associated with the upregulated genes. This was expected since the number of upregulated genes observed after D. rotundifolia 3× treatment was almost 3 times higher than the number of downregulated genes. Moreover, genes that were upregulated by D. rotundifolia 3× are associated with many biological processes, while this was not the case for the downregulated genes.

Discussion
In this study, functional genomic research was performed to explore the effects of a medicinal plant traditionally used for respiratory diseases on its potential target cells in vitro, i.e., bronchial epithelial cells. This study shows that a diluted herbal extract used in a traditional homeopathic pharmacopoeia 5,[9][10][11][12] can induce a change in mRNA expression that can be measured reliably and reproducibly. Low doses of traditional plants and homeopathic dilutions are now being investigated by means of molecular biology techniques for differential gene expression analysis and bioinformatics [24][25][26][27][28] .
Screening in bronchial cells with RNA-seq analysis showed that low amounts of D. rotundifolia at a 3× dilution, which represents a dose commonly used in homeopathic and phytotherapeutic syrups, changed the expression of dozens of genes after 3 h, and this effect was amplified after 6 h of treatment. Validation with RT-qPCR confirmed the differential expression of the genes of interest and showed that after 6 h of D. rotundifolia 3× treatment, there was an increase in mRNA expression. Since RT-qPCR sensitivity is higher than that of RNA-seq, such results support the findings of the RNA-seq analysis, suggesting that treatment with D. rotundifolia × can be followed-up by the investigation of the expression of these genes.
In our experiments, D. rotundifolia 3× did not impair cell viability or adherence, suggesting that this dilution is safe when tested directly in cell culture. The same dilution exerted a stimulatory effect on the expression of several genes, including those of inflammatory cytokines. Furthermore, the transcriptomic analysis and RT-qPCR confirmation suggest that, at the dilution used, D. rotundifolia 3× works mainly as a stimulant and not as an inhibitor of cell functions. This result is different from the data reported by others in HMC-1 human mast cells 29 , suggesting the anti-inflammatory effects of the plant. This apparent discrepancy could be due to the difference in the biological model, since mast cells are typical inflammatory cells and were stimulated before treatment with D. rotundifolia, while our epithelial cells were treated without previous inflammatory stimuli. However, even more importantly, different doses may be utilized: Fukushima et al. 29 obtained an inhibitory effect by using nondiluted fractions of the whole extract of D. rotundifolia, while in the present study, the cells were treated with a 3× dilution (1000 times) of the ethanolic extract. Table 3. Representative biological processes and associated differentially expressed genes (DEGs) in 16HBE cells treated with D. rotundifolia 3 × for 3 h. Table reports enriched biological processes identification (GO term and Biological process name), the number of differentially expressed genes belonging to each biological processes and their gene symbols. p-values indicate the enrichment significance obtained with g_SCS algorithm of gprofiler2 analysis, as described in the "Materials and methods" section.  CYP1B1, IL1A, SERPINB2, PPP1R10, EREG, ERRFI1, AREG, EGR1,  PPP1R15A, NCOA7, EPGN, DDIT4, HLA-F, IER3, LGR5, GIPR,  DYSF, NOG, AJUBA, JUN, UBR5, ARRDC3, CXADR, NLGN2,  HIF1A, GPAM, ATF4, LTA, IGFBP3, CTDP1, ADM, 30 reported that extracts of D. rotundifolia and other species, tested in hen's egg model, show efficacy as anti-inflammatory, antispasmodic and antiangiogenic agents. In contrast, the data of our study suggest that treatment with D. rotundifolia 3× can trigger a mild inflammatory response in 16HBE cells. The latter effect can be attributed to the increase in the mRNA expression of pro-inflammatory cytokines such as IL-1α, IL-1β and IL-6 and chemokines such as CXCL1, CXCL2 and CXCL8. Discrepancies in the data from Paper et al. 30 could be due to the different models or, again, to the doses applied, since we used a homeopathic dilution of the plant extract, which represents a biological stimulus and does not inhibit key cell functions 24 , as conventional antiinflammatory agents do. Moreover, it is well known that a mild inflammatory response could be beneficial to the immune system by increasing the apoptosis and clearance of inflammatory cells 31,32 . This concept is fully in agreement with homeopathic theory and tradition, in which low doses of pathogenic substances and/or minimally stressful stimuli trigger endogenous responses of "vital energy", which eventually lead to healing at the cell, tissue or systemic levels [33][34][35] . The presence of CXCL1, CXCL2 and CXCL8, which are mostly chemoattractants for neutrophils, suggests that treatment with D. rotundifolia 3× can help to establish a ready-state immune barrier against pathogens in case of bronchial epithelial tissue damage 36 .
Interestingly, D. rotundifolia 3× treatment increased the expression of specific epidermal growth factors, such as AREG, EREG and EPGN, by 16HBE cells. The upregulation of these factors suggests that D. rotundifolia 3× plays a positive role in the regulation of cell survival, cell proliferation and wound healing. Moreover, all these factors bind to epidermal growth factor receptor (EGFR) to accomplish their functions. In this context, EGFR is widely expressed on the cells present in the bronchial environment, suggesting that the possible release of these growth factors affects not only the proliferation of bronchial epithelial cells but also that of endothelial cells, fibroblasts and vascular cells. In the case of damage to bronchial epithelial tissue, D. rotundifolia 3× could contribute to faster healing of the bronchial microenvironment.
CYP1B1 was the most highly expressed gene after treatment with D. rotundifolia 3× at both time points shown in this study. CYP1B1 is probably activated as a response to treatment with D. rotundifolia 3x, since CYP1B1 Table 4. mRNA expression in 16HBE cells incubated with D. rotundifolia 3 × or Ctrl for 3 h. 16HBE cells were cultured with D. rotundifolia 3 × or the control for 3 h, and their mRNA expression of the indicated genes was evaluated by RT-qPCR. Gene expression is depicted as the mean normalized expression (MNE) after GAPDH mRNA normalization (mean ± SD of the indicated number of experiments www.nature.com/scientificreports/ belongs to the cytochrome family and is involved in the metabolism of xenobiotics. Moreover, as shown in Fig. 5, CYP1B1 is connected to the regulation of angiogenesis and the positive regulation of cytokine production because of its role as an antioxidant and NF-kB regulator 37,38 , respectively. This study has some limitations since some of the findings are based on bioinformatics analysis, and mRNA expression and functional experiments are required to verify the mechanisms of action and effects of the treatment with D. rotundifolia 3×. In addition, it is important to observe that after 6 h of D. rotundifolia 3× treatment, there was a number of differentially expressed genes with log2FoldChange values ranging from + 0.4 to − 0.4, which were not included in the analysis for this manuscript but could be considered for a future study. Indeed, the effects on gene regulation represent only a first functional step of the action of the plant that we have highlighted, and our hypotheses will have to be confirmed with adequate studies in laboratory animals to investigate the therapeutic potential and the mechanisms of action of this plant. In conclusion, the data highlight the complex and multifaceted action of the plant on the different layers of the bronchial mucosa, as summarized in Fig. 6. The mRNA expression changes in 16HBE cells treated with D. rotundifolia 3× suggest its direct action on the epithelial cell, which protects its integrity with respect to toxic substances (CYP1B1) and stimulates its reparative capacity (AREG, EREG and EPGN). In addition, epithelial cells transmit molecular signals that activate mild inflammation (pro-inflammatory cytokines) and recruit innate defence-and angiogenesis-associated cells. Before the experiments, D. rotundifolia 1× ethanolic extract was serially diluted 1:10 (0.5 mL + 4.5 mL) twice in sterile pyrogen-free water (B-Braun, Melsungen, Germany) in a 14-mL clean glass tube, which immediately vigorously shaken (succussed) with a Dyna A mechanical shaker that delivered 20 strokes/s for 7.5 s with an 11-mm travel distance. The final solution corresponds to a 3× dilution in 0.45% ethanol solution. The control solution was prepared starting from 45% ethanol (AppliChem, Darmstadt, Germany) fresh solution and succussed as previously described. This was considered the 1× control solution. Two serial decimal dilutions/succussions in ultrapure water were applied as reported above to obtain the control 3× dilution (Ctrl). No filtration was applied at any step. These solutions were used in the culture media at a 1:10 ratio (0.1 mL test solution in 0.9 mL culture medium). Therefore, the final dilution of D. rotundifolia was 1000 times greater than that of the ethanolic extract and the ethanol concentration was 0.045%.

Phenolic quantification.
The total phenolic content of the D. rotundifolia ethanolic extract was determined by the Folin-Ciocalteu assay 39 . Briefly, 50 μL of extract at an appropriate dilution was mixed with 155 μL of Folin-Ciocalteu reagent diluted 1:10 v/v with water. After 1 min, 40 μL of 20% Na 2 CO 3 solution was added, and the samples were incubated for 30 min in the dark at 37 °C. The absorbance of each sample was measured at 765 nm. Gallic acid was used as a standard for the calibration curve, and the phenolic content was expressed as gallic acid equivalents. Each determination was repeated three times, and the results are expressed as the mean ± SD.   Evaluation of cell viability. Cell viability was checked by the cell proliferation reagent WST-1 (Gibco ThermoFisher Scientific). The WST assay evaluates the cell metabolic activities (NADH reductase) by measuring the chemical modification of the WST tetrazolium salt. 16HBE cells were seeded at a density of 50,000 cells/well in 96-well plates and treated for 3 or 24 h with the D. rotundifolia 3× or the Ctrl solutions, which were added to the cell culture at a 1/10 volume ratio. After treatment, 1:10 (v/v) prewarmed WST-1 solution was added to the cells, and the plate was incubated for 90 min. The absorbance (OD) of the samples was measured using a Vic- Figure 5. Functional network of enriched biological processes and associated genes. The network was constructed using the selected significant GO terms shown in Fig. 4, which were uploaded into Cytoscape software as described in "Materials and method" section. Diamonds indicate biological processes and circles indicate associated genes. Transcript quantification was performed with the high-performance computing server at the Computational Unit of the Technologic Platform Centre (University of Verona, Italy). Transcript quantification was determined from the FASTQ reads using the mapping-based mode of Salmon (Version 0.13.1) 43 with the following parameters: "-i index -l SR -r name -validateMappings -gcBias". The reference sequence and annotation files were downloaded from the GENCODE repository (Human Release 32, GRCh38.p13) at the following website: https ://www.genco degen es.org/human /. The output of Salmon was then imported into the RStudio environment (R Version 3.5.3; R studio version 1.1.463) using Bioconductor/R package tximport (Version 1.10.10) 44 , which converted data from the transcript level to the gene level during the importing process. Genes with less than 10 reads were discarded from the posterior analysis. Differential expression analysis was performed using the Bioconductor/R package DESeq2 (Version 1.22.2) 45 . After differential analysis, the lfcShrink function in DESeq2 was applied to shrink the log2FoldChanges. Significance values were based on a Wald significance test, and differences in gene expression with an adjusted p-value < 0.05 (corrected by the Benjamini-Hochberg method) were considered significant. For the RNA-seq analysis of the experiments, in which the cells were treated for 6 h with Those genes, which bind to EGFR to accomplish their functions, are regulated by ERFFI1, which is in turn downregulated by D. rotundifolia 3× in the first 3 h. CYP1B1 is highly expressed and is related to the metabolism of xenobiotics and could be involved in the regulation of angiogenesis. Furthermore, inflammatory cytokines (IL-1α and IL-6) and chemokines (CXCL1, CXCL2 and CXCL8) are induced by D. rotundifolia 3× and can trigger mild inflammation, increasing chemotaxis and angiogenesis, helping the system in fighting infections. www.nature.com/scientificreports/ D. rotundifolia 3×, a log2FoldChange cut-off value (< 0.4 and > 0.4) was applied. To associate a function with the differentially expressed genes, functional enrichment analysis was performed using the Bioconductor/R package gprofiler2 (Version 0.1.7) 46 . The differentially expressed genes were ranked by the adjusted p-value significance, and gene ontology (GO) analysis was queried using the "gost" function in gprofiler2 in "ordered" mode and "g_SCS" as correction/statistical method. The list of significant GO terms was used for the construction of the network using Cytoscape software (Version 3.7.2) 47 .

Reverse transcription quantitative real-time PCR (RT-qPCR). Total RNA extracted from 16HBE
cells was reverse transcribed into cDNA using the PrimeScript RT reagent Kit (Takara Bio, Kusatsu, Japan), while qPCR was carried out using TB Green Premix Ex Taq (Tli RNase H Plus) (Takara Bio) 42 . The sequences of the gene-specific primers (Thermo Fisher Scientific, Waltham, Massachusetts, USA) used in this study are listed in Table 5. qPCR was performed using a Viia7 Real-time PCR system (Thermo Fisher Scientific). Data were calculated by Q-Gene software (http://www.gene-quant ifica tion.de/downl oad.html) and expressed as the mean normalized expression (MNE) units after GAPDH normalization 42 . Statistical evaluation was performed by twoway ANOVA followed by Bonferroni's post hoc test. Values of p < 0.05 were considered statistically significant.
Ethics approval and consent to participate. Not applicable, as the study did not involve patients, volunteers or animals.

Data availability
The raw and processed RNA-seq data are available in the public GEO database (https ://www.ncbi.nlm.nih.gov/ geo/) under the accession number GSE144215.