In this study, we employed four different cell systems: the cell lines HEK-293, SH-SY5Y, and PMA-differentiated THP-1, and primary human epidermal keratinocytes.

The human embryonic kidney cells HEK-293 (Merck, Buchs, Switzerland) were seeded at 60,000 cells per well in a 96-well plate and cultured in Eagle’s minimum essential medium (EMEM; Gibco, Grand Island, NY, United States) supplemented with 10% FBS (Sigma, Burlington, MA, United States), 1% non-essential amino acids (Merck), 1% sodium pyruvate (Gibco), and 1% penicillin/streptomycin (PenStrep; Gibco).

SH-SY5Y is a human bone marrow neuroblastoma cell line that has previously been used to investigate the potential anti-inflammatory effect of anatabine (Paris et al., 2013a). SH-SY5Y cells (CRL-2266, ATCC, Manassas, VA, United States) were seeded at a density of 15,000 cells per well in a 96-well plate and cultured in a medium consisting of 50% EMEM and 50% Ham’s F12 nutrient mix (Gibco), along with 10% FBS, 1% PenStrep, and 1% L-glutamine (Gibco).

The human leukemia-derived monocytic cell line THP-1 is the most widely used in vitro model for primary investigations on human macrophages (Lund et al., 2016). THP-1 cells (Merck) were seeded at 100,000 cells per well in a 96-well plate and cultured in RPMI 1640 medium (Gibco) supplemented with 10% FBS and 1% PenStrep and in the presence of 40 ng/ml phorbol 12-myristate 13-acetate [PMA; Thermo Fisher (Kandel) GmbH, Kandel, Germany]. After 48 h, the culture medium was replaced with fresh medium without PMA and allowed to rest for 24 h before treatment, as previously described (Lund et al., 2016).

Primary human epidermal keratinocytes were included as a relevant cell system for skin. The keratinocytes (Ruwag, Bettlach, Switzerland) of healthy, non-smoking, Caucasian donors aged 20–60 years old, were selected from a donor pool. The keratinocytes were seeded at 17,500 cells per well in a 96-well plate and cultured in KGM-Gold Keratinocyte Growth Medium (Ruwag).

All the four cell systems presented above were seeded at a total volume of 190 μL per well to allow for the addition of 10 μL of 20× anatabine, which was added 24 h after seeding, or the addition of non-PMA media in the case of THP-1 cells.

The same anatabine formulation (WuXi Apptec, Shanghai, China) consisting of racemic anatabine free base (hereafter referred to as “anatabine”; MW: 160.22 g/mol) from a single production batch (Batch B) was used across all experiments.

Anatabine concentrations of less than 400 μM were selected to ensure that at least 80% cell viability was observed after 24 h of exposure. Viability was assessed using CellTiterGlo (Promega, Madison, WI, United States) by measuring the ATP content in each treated sample and comparing them with those of relevant controls, which were typically cells treated only with cell culture medium.

To generate transcriptomics data, cells were exposed to different concentrations of anatabine (0, 100, 200, 300, and 400 μM) for 6 h, and their transcriptome was analyzed using a microarray-based technique.

In detail, all three cell lines were lysed using RLT buffer (Qiagen, Hilden, Germany), and keratinocytes were lysed using the Qiazol lysis reagent (Qiagen). RNA isolation was performed with the RNeasy Micro Kit (Qiagen) on a Qiacube instrument (Qiagen). RNA was quantified with Nanodrop 1000 (Thermo Fisher Scientific, Waltham, MA, United States). The quality of the total RNA, which was required to have an RNA integrity number greater than 6.0, was assessed using 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). Total RNA (50 ng) was processed in the Tecan/Nugen Ovation RNA Amplification system V2 kit (Tecan, Männedorf, Switzerland), according to the manufacturer’s instructions, followed by cDNA fragmentation and labeling using the Encore Biotine module (Tecan). The Human Genome U133 Plus 2.0 microarray was used for hybridization using a Thermo Fisher Oven 645. Furthermore, washing was performed on Affymetrix GeneChip™ Fluidics Stations 450Dx (protocol FS450-0004) and scanning on a ThermoFisher GeneChip™ Scanner 3000 7G.

The Bioconductor affyPLM package in R version 1.64 (Bolstad et al., 2004) was used for quality control checks of all chips. Following quality control procedures, differential gene expression was analyzed using the Bioconductor limma package in R version 3.44.3 (Ritchie et al., 2015). Pairwise comparisons at the gene level, called systems response profiles (SRP), were computed by comparing each concentration-treatment with its respective vehicle control. Genes with a false discovery rate (FDR) (p-value adjusted using the Benjamini and Hochberg method) below 0.05 were considered differentially expressed genes (Benjamini and Hochberg, 1995). SRPs including all genes (∼18,000) were further leveraged in downstream p-value threshold-free gene set enrichment analysis (GSEA), as described in Section 2.6.

Samples for DIA were prepared using the PreOmics iST kit (PreOmics GmbH, Planegg/Martinsried, Germany), according to the manufacturer’s protocol. Briefly, 100 µL of PreOmics lyse buffer was added to the cells, and the mix was incubated at 90°C for 10 min and sonicated for 30 s with a sonifier (Branson, Danbury, CT, United States) at 10% amplitude. Protein concentration was determined using the Pierce 660 nm protein assay (Pierce Biotechnology Inc., Rockford, IL, United States), according to the manufacturer’s protocol. The samples were normalized to 0.5 μg/μL, and 40 µg of each sample was further processed with the PreOmics iST kit with a 3-h long trypsin digestion. Peptides were purified on the cartridge, dried overnight on a vacuum concentrator (Martin Christ, Osterode, Germany), and resuspended in 50 µL of LC Solution (Biognosys AG, Schlieren, Switzerland). iRT reference peptides (1 μL; Biognosys AG) were added to 19 µL of the processed samples, which were analyzed with Easy 1000 nanoLC (Thermo Fisher Scientific) connected online to a Q-Exactive mass spectrometer (Thermo Fisher Scientific). Two microliters of the peptide mixture were separated on a 0.75 × 500 mm, 1.7 μm C18 column (Thermo Fisher Scientific) using solvent A (1% acetonitrile/99% water/0.1 formic acid) over 115 min with a gradient of 5%–35% solvent B (95% acetonitrile/5% water/0.1 formic acid) at 200 nL/min.

Data were acquired on the Q-Exactive system in the DIA mode: MS1 scan at a 140 k resolution was followed by 23 custom MS/MS m/z windows at a 35 k resolution, as previously described (Bruderer et al., 2017) with slight modifications. Data were processed with Spectronaut Pulsar (v. 13.8.190930.43655; Biognosys AG) using the DirectDIA feature.

SRPs were computed by comparing samples at each concentration with their respective vehicle control, as described in Section 2.2.

Cells were exposed to 100, 200, 300, and 400 μM of anatabine for 15 min, 24 min, 6 h, and 24 h. For each time point, the cells were placed on ice, washed with phosphate-buffered saline, and lysed with the addition of lysis buffer (Protavio Ltd, Cambridge, United Kingdom). The obtained lysed samples were frozen at –80°C. Just before the phosphoproteomic measurement, the samples were quickly thawed and centrifuged for 20 min at 2700 g to remove cellular debris. The xMAP assays were performed using a custom-made 18-plex phosphoprotein detection kit (Protavio Ltd), according to the kit instructions, and a Luminex FLEXMAP 3D instrument (Luminex, Austin, TX, United States).

Data were acquired and pre-processed as described previously (Michailidou et al., 2015). SRPs were computed as described in Section 2.2.

An HEK-293-based NRF2/ARE luciferase reporter cell line (Signosis, Santa Clara, CA, United States) was used to evaluate NRF2 activation. The cell line was stably transfected with pTA-ARE-luciferase reporter vector, which contains four repeats of the antioxidant response binding site, a minimal promoter upstream of the firefly luciferase coding region, along with a hygromycin expression vector. Following a 24-h stimulation of the NRF2 reporter cells with various treatment conditions, luminescence was measured using a FLUOstar Omega plate reader (BMG LABTECH, Ortenberg, Germany). Dimethyl fumarate (Alfa Aesar and Sigma) and sulforaphane, two NRF2 activators, were used as positive controls for the assay. SRPs were computed as described in Section 2.2.

Linear model analysis was used to identify genes whose expression changed linearly with an increase in anatabine concentration (0, 100, 200, 300, and 400 μM). The linear model analysis was conducted using the transcriptomics data acquired either from each cell system separately (gene expression ∼ β₀ + β₁*anatabine_concentration + ε) or all cell systems combined (gene expression ∼ β₀ + β₁*anatabine_concentration + β₂*cell_systems + ε) to identify cell system-specific and cell system-independent genes (“core genes”), respectively. Genes with β₁ coefficient-associated FDR <0.05 were considered to have a linear change in their expression with increasing concentrations of anatabine. The limma package (version 3.38.3) in R (version 3.5.1) was used for the analysis.

Leveraging these results, GSEA was performed to support biological interpretation (Subramanian et al., 2005). Genes were ranked in decreasing order based on their respective β₁ coefficient-associated t-statistics. The “C2-CP” (Canonical Pathways) and “C3-TFT” (Transcription Factor Targets) MSigDB (version 7.1) gene set collections were used as sources of a priori biological knowledge (Liberzon et al., 2015). Gene sets with normalized enrichment score (NES)-associated FDR <0.05 were considered to be significantly enriched.

To investigate the mode of action of anatabine, gene expression profiles obtained in response to anatabine were compared with those obtained in response to a plethora of perturbagens corresponding mostly to chemical/drug compounds (Messinis et al., 2021). These profiles are part of the publicly available level 5 LINCS L1000 dataset (GEO accession number GSE92742), consisting of 473,647 differential gene expression signatures (profiles), created using various concentrations of 28,927 unique perturbagens for treating multiple cell systems. Each LINCS L1000 gene expression “signature” corresponds to the measurement of 978 genes. Genes from anatabine expression profiles generated using oligonucleotide microarray were matched with those from the L1000 dataset, resulting in 938 genes in common that were used for further comparison analysis.

Spearman correlation was used to compare L1000 and genes from expression signatures of anatabine-treated systems. The ranking of genes for anatabine expression signature was based on their respective β₁ coefficient-associated t-statistics computed using a linear model combining all cell systems. Among the L1000 perturbations, there were several cases where the same compound was used in different experimental conditions. In those cases, the comparison with the highest Spearman correlation coefficient was retained. We ran 100 random transcriptomic signatures against all L1000 signatures in parallel, calculated the significance of the hypothesis that anatabine’s Spearman correlation coefficient values are in the same range as those of the random signatures’, and discarded any observation with a p-value greater than 0.001.

To compare the chemical similarity among compounds, the compounds were encoded using several chemical similarity fingerprints and descriptors. The chemical similarity was compared using the DataWarrior software (Sander et al., 2015). DataWarrior was used to calculate the chemical similarity descriptors “Fragment Fingerprints” (FragFP), “Pathway Fingerprints” (PathFp), “OrgFunctions,” “Spheres Fingerprints” (SphereFp), “SkeletonSpheres” (SkelSpheres), and “Flexophore.”

FragFP (Durant et al., 2002) is a substructure fragment dictionary-based binary fingerprint similar to MDL Keys (Xue and Bajorath, 2000), which contains 512 predefined structure fragments. The FragFP descriptor contains 1 bit for each fragment in the dictionary. A bit is set to 1 if the corresponding fragment is present in the molecule at least one time.

PathFp (O'Boyle et al., 2011) encodes any linear fragment of up to seven atoms into a hashed binary fingerprint of size 512 bits. All possible “paths” of seven or less atoms in the molecule are encoded. A text string that encodes atomic numbers and bond orders is generated from the path in a normalized way. From this text string, a hash value is generated that is used to set the corresponding bit of the fingerprint to 1.

SphereFp (Bremser, 1978) encodes circular spheres of atoms and bonds into a hashed binary fingerprint of size 512 bits. Fragments of increasing size are generated by including 1–5 n layers of atom neighbors for every atom in the molecule. A canonical representation of these circular fragments is obtained considering their aromaticity but not their stereo configurations. A hash code is then generated, which is used to set the respective bit of the fingerprint.

SkelSpheres (Boss et al., 2017) is similar to SphereFp but also encodes stereochemistry. Additionally, it counts duplicate fragments, encodes heteroatom depleted skeletons, and has twice the resolution leading to fewer hash collisions (bit size: 1024).

OrgFunctions (Sander et al., 2015) encodes the functional groups in the molecule and the steric or electronic features of the neighborhood of the functional groups. It encodes 1024 core functions that overlap in some cases. Molecules that contain the same functional groups are considered similar irrespective of the possible presence of carbon skeletons.

Flexophore (von Korff et al., 2008) encodes 3D-pharmacophore features of the molecules. This descriptor can be used to check if two compounds have a similar protein binding mode. A high Flexophore similarity of two molecules indicates that a significant fraction of conformers of both molecules is similar with regard to shape, size, flexibility, and pharmacophore features. In addition to a 3D-pharmacophore model, the Flexophore descriptor matches entire conformer sets rather than comparing individual conformers.

Three-dimensional pharmacophore models of subsets of the 15 compounds were built using the Molecular Operating Environment (MOE) software (2019.01; Chemical Computing Group ULC, Montreal, Canada). The pharmacophores were built using the “Pharmacophore Elucidation” tool in MOE. First, a complete conformational search of each individual compound was performed, and similar conformers were removed. The compounds of interest were then aligned, and the best alignment was selected for implementation of the most relevant pharmacophore features of the selected compounds. The alignment was performed giving emphasis to aromatic, donor, and acceptor atoms. The pharmacophore features were generated using the default list of possible features and the default radius of each feature. Finally, from this set of suggested features, the “Pharmacophore consensus” tool was used to select the most relevant pharmacophore features.

Transcriptomics data were used to infer a network model representative of mechanisms of action of anatabine. As a first step, transcriptomics data corresponding to the outcomes of the linear model including all anatabine concentrations and combined cell systems were reverse-engineered to predict potential active transcription factors using the DoRothEA algorithm (Garcia-Alonso et al., 2019). The top 10 most actively predicted transcription factors were selected for further network construction. Ten is an arbitrarily selected number, so as to obtain a resulting network that does not have too many nodes and therefore is possible to visualize.

As a second step, the CARNIVAL algorithm (Liu et al., 2019) was used to build a network leading to the 10 predicted transcription factors and starting by a hypothesized perturbation, which was anatabine in our case. CARNIVAL requires a prior knowledge network, on which the experimental data will be fitted. Two publicly prior knowledge networks were used: OmniPath and Reactome Functional Interactions (FI). OmniPath (Turei et al., 2016) contains 164,710 interactions, whereas Reactome FI (Wu et al., 2010) consists of 259,151 interactions. Both networks have been assembled utilizing numerous resources and are updated regularly. We used the latest versions available for both (version 2.0.0 of OmniPath R package, and the 2020 version of Reactome FI).

The CARNIVAL algorithm generates an integer linear programming problem based on the provided prior knowledge network and the predicted transcription factors and then solves it with the solver IBM ILOG CPLEX. The solution includes the predicted network topology along with a value for every node of the network ranging from 0 to 100, which corresponds to the node’s predicted level of activation. The network topology was visualized using Cytoscape version 3.8.2 (Shannon et al., 2003).

HEK-293 cells, SH-SY5Y cells, PMA-differentiated THP-1 cells, and primary human keratinocytes were treated with various concentrations of anatabine or vehicle control for 6 h. At the end of the treatment, cell lysates were collected and processed to generate transcriptomic profiles.

A principal component analysis enabled us to investigate the sources of variability in the data (Figure 1A). Principal components (PC) 1 and PC2 explained 19.6% and 17.4% of the variability, respectively, highlighting the combined effects of cell type (prominent effect) and anatabine concentrations.

A linear modeling analysis enabled us to identify genes whose expression changed linearly with increasing concentrations of anatabine, for each individual cell system or all cell systems combined (Figures 1B,C). Figure 1B illustrates, as an example, the linear relationship between anatabine concentrations and changes in the expression of the gene heme oxygenase 1 (HMOX1) for each of the tested cell systems. The slope of each depicted linear regression line portrays the modification in corresponding gene expression per unit of anatabine concentration. Figure 1C shows the top 30 most significantly expressed genes.

To contextualize the processes and pathways represented by these genes, we conducted a GSEA using various gene set collections, namely “C2-CP” and “C3-TFT”. Figure 1D depicts the 30 most significantly enriched gene sets resulting from the GSEA using the gene list generated from comparing all cell systems together. The most significant positively enriched gene set pertains to the nuclear factor-erythroid 2 p45-related factor 2 (NRF2), while the next three are related to the transcription factor heat shock factor 1 (HSF1).

To obtain a more systematic view of anatabine’s activity, we investigated intracellular protein perturbations caused by anatabine, on the PMA-differentiated THP-1 cell system, employing DIA mass spectrometry proteomics.

Figure 2A displays a heatmap of significantly abundant proteins for different concentrations of anatabine. Figure 2B shows the abundance of each protein along with its corresponding gene expression. We observe that the protein abundance of the four most significantly upregulated proteins, HMOX1, glutamate-cysteine ligase modifier subunit (GCLM), thioredoxin reductase 1 (TXNRD1), and NAD(P)H quinone dehydrogenase 1 (NQO1), correlated with the corresponding gene expression changes (Pearson correlation coefficients: 0.97, 0.92, 0.96, and 0.95, respectively).

Based on the observation that these four proteins are NRF2 target genes (Tonelli et al., 2018), we hypothesized that the NRF2 pathway might be involved in anatabine’s mechanism of action.

To test our hypothesis, we used an HEK-293 NRF2/ARE luciferase reporter cell line that assesses whether NRF2 translocates to the cell nucleus and binds to its response element upon treatment with the compound. Along with anatabine, we tested the tobacco alkaloids anabasine, cotinine, nornicotine, and nicotine (to investigate the possibility of a shared effect among tobacco alkaloids) and the NRF2-activating compounds sulforaphane and dimethyl fumarate (to serve as the assay’s positive controls).

The results show a concentration-dependent NRF2 activation, which reached statistical significance (p < 0.05) at 250 μM of anatabine treatment (Figure 3). Interestingly, anabasine, cotinine, nornicotine, and nicotine did not activate NRF2, demonstrating that NRF2 activation is a unique characteristic of anatabine among the tested tobacco alkaloids.

To better understand how anatabine activates NRF2 and the molecules that may act as its functional partners, we leveraged the LINCS L1000 dataset that consists of 473,647 differential gene expression signatures (Subramanian et al., 2017), originating from cell systems treated under various conditions with 28,927 different perturbagens, such as chemical compounds, recombinant proteins, or small interfering RNAs. Upon comparing each of those signatures with anatabine’s gene expression profile, we could rank the LINCS perturbagens based on their similarity to anatabine’s profile.

Of the top 15 compounds with gene expression signatures most similar to those of anatabine (Table 1), six have not yet been referenced in the literature, while two other compounds were mentioned in only four studies. The remaining seven compounds are natural compounds that have been shown to induce oxidative stress response. Among them, four compounds have already been identified as NRF2 activators (piperlongumine, withaferin-A, cucurbitacin-I, and parthenolide).

We observed common structural elements, such as elements of hydropyridine, piperidine, and quinoline, among the 15 compounds. To examine if the structural elements were responsible for the observed similarities in the gene expression induced by these compounds and anatabine, we employed six different chemical similarity fingerprints and descriptors using the DataWarrior software (Sander et al., 2015). Each compound fingerprint was compared to the corresponding fingerprint of anatabine, and the Tanimoto distances of the structural comparison are presented as a heatmap in Figure 4A.

The correlation between gene expression similarity scores and the chemical similarity fingerprints was examined but not found significant. Compound GP-42 exhibited the highest overall Tanimoto distance when compared with anatabine. Therefore, we visualized the common pharmacophore features of GP-42 and anatabine. A 3D pharmacophore model for GP-42 was generated (Figure 4B). Anatabine was fitted in the pharmacophore model, having one aromatic feature, one hydrophobic feature, and two acceptor features in common with GP-42, demonstrating that anatabine and GP-42 could share a similar mode of binding to some degree.

We performed network analysis to extract further insights from the transcriptomics data, leveraging the Reactome and Omnipath biological networks (Wu et al., 2010; Turei et al., 2016), the DoRothEA enrichment algorithm (Garcia-Alonso et al., 2019), and the CARNIVAL optimization algorithm (Liu et al., 2019). The resulting networks indicate potential functional partners of anatabine.

The network analysis revealed p38 MAPK as a central node in the signaling network (Figure 5A). Interestingly, p38 MAPK was predicted to be either positively or negatively regulated depending on the prior knowledge network used in the analysis (Reactome or Omnipath, respectively), even though it was centrally located in both networks. We also observed the presence of dual-specificity phosphatases (DUSP) upstream of the MAPK signaling in both networks. NRF2 was inferred to be activated in both the networks.

To explore the potential involvement of p38 MAPK in facilitating the effects of anatabine, we used multiplexed phosphoproteomic assays. We then applied a linear model approach to capture time- and concentration-dependent phosphoproteomic changes across the four tested cell systems. The phosphoproteomic results after application of the linear model are presented in Figure 5B. The phosphoproteomic data before application of the linear model are detailed in Supplementary Figure S1B.

The phosphoproteomic results, summarized via the linear model across time, concentrations, and cell systems tested, revealed a systematic significant increase in the phosphorylation levels of p53, p38 MAPK, ERK1, MEK1, MARCKS, p90RSK, JNK, c-Jun, HSP27, PRAS40, and p70S6K. Differences among cell systems, concentrations, and treatment times are depicted in Figure 5B and Supplementary Figure S1A. Anatabine treatment also caused a significant decrease in the phosphorylation levels in some of the cell systems tested, specifically, a concentration-dependent inhibition of phosphorylated STAT3 and a time-dependent inhibition of phosphorylated NF-κB (Supplementary Figure S1A).