Naringenin, hesperetin, narirutin, naringin, hesperidin, neohesperidin, nobiletin, neohesperidin dihydrochalcone, and tert-butylhydroquinone were purchased from Yuanye (Shanghai, China). Dimethyl sulfoxide (DMSO), 3-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT), and bovine serum albumin (BSA) were purchased from Fluka (Sigma-Aldrich, Shanghai, China). Fetal bovine serum (FBS) and Dulbecco’s modified Eagle medium (DMEM) were provided by Gibco (Grand Island, NY, USA). Lipofectamine 2000 transfection reagent was obtained from Thermo Fisher Science (San Jose, CA, USA). All common reagents used in this study were of analytical grade.

The reporter plasmid (pARE-Luc) was produced by ligating a synthetic double-stranded oligonucleotide tandem ARE target fragment containing the Nrf2 binding site spanning 5’-TGACTCAGC-3′ into a luciferase reporter gene plasmid vector (25). The renilla luciferase was normalized using the control plasmid pRL-SV40 (Promega, WI, USA).

Human embryonic kidney cell line HEK293T was grown in DMEM enriched with 10% FBS and antibiotics: 100 IU/mL penicillin and 0.1 mg/mL streptomycin. Cells were seeded into 96-well plates and kept in a humidified condition of 5% CO₂ in a 37°C incubator for 12 h. Cell viability was tested via MTT assay as per the method from (21) with modifications (26). Briefly, after the treatment with the tested chemicals, 10% MTT (5 mg/mL) was added and cultured for 4 h at 37°C. After that, the medium was removed, and the formazan crystals were dissolved in 150 μL of DMSO and then quantified by a SpectraMax i3x microlate reader (Molecular Devices, Sunnyvale, CA, USA) at 570 nm.

According to a previous method by Simmons, reporter gene assays were conducted with minor modifications (27). The 1 × 10⁵ cell/well HEK293T cells were plated at the 24-well plate in phenol red-free DMEM added with 2% dxtran-charcoal plus 10% FBS. According to the manufacturer’s instruction of Lipofectamine 2000, the cells were immediately transfected with plasmids with pARE-Luc and pRL-SV40 after 18 h of incubation. Hereafter, the citrus-derived flavonoids at the level of 10, 20 or 40 μM were added and cultured for 24 h, followed by cell lysis using a passive lysis buffer. Meanwhile, the tert-butylhydroquinone, one of the most commonly used phenolic antioxidants, was used as the positive control. Based on the maker’s protocol, both the firefly and Renilla luciferases in cell lysates were tested using a Dual-Luciferase® reporter assay system (Promega, Madison, WI, USA). The luciferase activities were measured using a SpectraMax i3x microlate reader (Molecular Devices).

In study of citrus-derived flavonoids, DFT was chosen in all theoretical computations (28, 29). B3LYP hybrid functional is used with the standard 6-311G (d, p) basis set and with ++ diffuse function for the high-quality theoretical geometry optimization via Gaussian 09 (30). For visualization Gauss view 06 is applied. The important molecular orbitals, such as the lowest unoccupied molecular orbital (LUMO), highest occupied molecular orbital (HOMO), and corresponding energy gap, were calculated to clarify the charge transfer possibility within the citrus-derived flavonoids (31). These computational details also provided valuable details on the binding behavior with the receptor.

To better understanding the chemical reactivity and stability of citrus-derived flavonoids, the molecular electrostatic potential (MSP) and the corresponding quantum descriptors were systematically calculated in detail (32). The visualized demonstrations of MEPs show the sites for electrophilic and nucleophilic location by using Gaussian 09 and Gaussview 06 software. And the auantitative analysis of MSP were performed by multiwfn, which is a powerful program for realizing electronic wavefunction analysis (33).

To investigate the binding behavior of the tested chemicals toward Keap1, molecular docking with grid-based energy assessment was performed with the help of docking program Autodock Vina (34). After modification by adding polar H atoms and Kollman charge, the crystal structure of Keap1 in the file 1X2R was used as a docking receptor (35). According to the crystal structure at the receptor binding site, the specific grid box with grid spacing 0.375 Å was generated to enclose all critical residues of Keap1. For each compound, the lowest binding energy (kcal/mol) was calculated based on the Lamarckian genetic algorithm and chosen for final ranking (36, 37). The binding modes between Keap1 and the tested compounds were visualized on PyMol and Discovery studio 2018 (38).

Drug-likeness and toxicity can forecast if a specific chemical has properties as an orally active compound (39, 40). To determine the Rule of Five for drug-likeness and toxicity, the canonical simplified molecular input line entry (SMILE) of ligand extracted from PubChem was fed into a freely accessible online tool SwissADME through a login-free website (41).¹

The data were reported as mean ± standard deviation (SD). All assays were carried out independently at least three times, and statistical analysis was conducted using one-way SPSS. All the data were from three independent repeats each treatment (*/#/$, **/##/$$, and ***/###/$$$: p < 0.05, p < 0.01 and p < 0.001).

In order to performed the luciferase reporter gene assays, the concentration range of citrus-derived flavonoids with no significant inhibited effect of HEK293T should be determined first. This confirmed the changes in ARE activation activity were due to xenobiotic compound exposure, rather than the side effect of cell death. The tested compounds at a concentration of 50 μM did not exhibit cytotoxic effect (Figure 1), therefore, 10, 20, and 40 μM were chosen for the subsequent experiments.

To assess the ARE-driven luciferase activity of citrus-derived flavonoids, ARE-driven reporter gene assay in vitro was constructed. The luciferase activity of the transient transfection system in HEK293T cells was measured after 24 h of incubation. The luciferase activity of HEK293T cells was significantly induced by 12.5 μM tert-butylhydroquinone, as shown in Figure 2. The citrus-derived flavonoids also induced transactivation of ARE in a dose-dependent manner. Two categories were classified: strong agonists (naringenin, narirutin, naringin, hesperetin) and weak agonists (hesperidin, neohesperidin, nobiletin, neohesperidin dihydrochalcone).

As a reliable approach for related organic molecules, DFT/B3LYP hybrid functional with 6-311G (d, p) was broadly applied to calculate and analyze the quantum chemistry parameters in this study (42, 43). Specifically, LUMO+1, LUMO, HOMO, and HOMO-1, which collectively called FMOs, play a crucial role in the chemical characteristics of the bioactive compounds (44, 45). Among them, the HOMOs and LUMOs represent the electron-providing and -accepting abilities of the tested chemicals (46). The three-dimensional plots and energy results of FMOs were tabulated in Figure 3 and Table 1. The negative and positive phases were represented in green and red, respectively. The computations found that the HOMOs and LUMOs of nobiletin were mainly localized over the entire molecule. The LUMO+1, LUMO, HOMO, and HOMO-1 energies are −0.02560, −0.6004, −0.21727, and − 0.22707, respectively. While the HOMO of naringenin is positioned on nearly the whole structure, its LUMO is particularly located on the phenol group. The energy gap represents the chemical reactivity and intramolecular charge of the compounds. A larger energy gap means lower kinetic stability, and vice versa. The DFT calculated energy gaps of naringenin, hesperetin, naringin, and narirutin were 0.1794, 0.17517, 0.17388, and 0.16772 eV respectively, which are associated with the decreasing energy gap and stand for the reactivity of the structures.

Based on the DFT computation at the same level of theory, the MEP of chemicals could visually represent the electronic distribution (47). The calculated MEPs in Figure 4 show a surface analysis throughout the tested compounds calculated with the B3LYP density functional method. Potential ranks as red < green < blue. The red region represents negative electrostatic potential, and indicate the proton is attracted by the aggregate electron density in the molecule and related to electrophilic reactivity. MEP surface plots have been traced with the help of the different colors as the given significance of electrostatic: the blue regions are the positive electrostatic potential, and suggest the proton is repelled by the atomic nuclei and associated with nucleophilic reactivity. The most electrophilic attacks occur on the C=O (Figure 4). Moreover, the blue-green region corresponds to the hydroxyl of benzene ring in the tested chemicals. The MEP results imply that the hydrogen and oxygen atoms suggest the strongest attraction and repulsion, respectively, which may affect the electrostatic interactions between ligands and surrounding residues of Keap1 and subsequently influence their binding potency.

Since the electrostatic potentials greatly account for the interaction behavior between the receptor and ligands, quantitative study of molecular surface is necessary. The global minimum on the surface of the tested compound occurred at C=O nearly −50 kcal/mol, which is the largest minimum value on the surface (Table 2). The global minimum arising from the cationic hydroxyl varies largely from 302.19627 to 671.84562 kcal/mol. In comparison, more electronegative O atoms connected to the H atoms attract the electron cloud, making the H atoms more positive. The internal charge separations of strong agonists are generally higher than weak agonist groups.

Docking calculations between the compounds and Nrf2 reproduced the tested binding conformations, as found in Keap1-Nrf2 complex (PDB: 1X2R), implying the Keap1 model is feasible for calculating binding free energy. Docking calculations showed the binding affinity of the tested compounds with Keap1 ranked as narirutin > naringenin > hesperetin > hesperidin > naringin = neohesperidin dihydrochalcone > neohesperidin > nobiletin (Table 3). In the 2D and 3D representations in Figures 5, 6, the dark green, light green orange, pink, red, and cyan represent hydrogen bond, van der Waals, carbon hydrogen bond, and Pi-cation, respectively (35). Based on the experiments above, these compounds can be divided into weak agonists (naringin, neohesperidin dihydrochalcone, neohesperidin, nobiletin) and strong agonists (narirutin, naringenin, hesperetin, hesperidin). Basically, the order of activation potency for compounds agrees well with the binding affinity with Keap1. The docking results showed the top three molecules (naringenin, hesperetin, narirutin) were well accommodated in the binding pocket of Keap1 with high hydrophobicity. Except the two unfavorable donor-donor interactions, the hydroxyl groups and carbonyl of naringenin form three hydrogen bonds with VAL⁵¹² and VAL⁴⁶⁵, which contribute to the high binding energy between Keap1 and naringenin. Similarly, in addition to the three hydrogen bonds with VAL⁵¹² and VAL⁴⁶⁵, the hydroxyl groups of hesperetin also form other hydrogen bonds with GLY³⁶⁷, VAL⁶⁰⁶, and ARG⁴¹⁵. Different from naringin, narirutin enters the binding pocket of Keap1 with the highest binding affinity and exhibits different interactions via hydrogen bonds, van der Waals, and Pi-Alkyl interactions (48). Additionally, the binding mode of weak agonists is strikingly different from that of the strong agonists. The weak agonists form several hydrogen bonds with the surrounding residues of Keap1, such as ARG³⁸⁰, SER^431, ARG⁴⁸³, ASN⁴¹⁴, ARG⁴¹⁵, LEU³⁶⁵, and TYR^572. The multiple hydrogen bonds as well as Pi-Alkyl and Pi-Sigma interactions account for the stability and closeness between the ligand and Keap1.

With the rapid increase in biological and chemical information, evaluation of physicochemical and ADME properties is vital in discovery of phytocompounds from traditional plants. To detect bioavailability potential of phytocompounds, the pharmacokinetics, bioavailability, drug-likeness and medicinal chemistry friendliness of the citrus-derived flavonoids were characterized primarily based on “Lipinski’s rule of five,” which are the molecular mass shall not surpass 500 Dalton, there shall not be over five hydrogen bond donors or over ten acceptors, and the partition index LogP shall not exceed 5. The Swiss ADME results (Table 4) show that naringenin, hesperetin, and Nobiletin satisfy Lipinski’s rule of five with zero violations. Others showing a violation of more than two rules depict a molecule as not orally active. Hence, the three tested compounds can be candidates for further studies.

The predictive data of individual ADME behaviors of the tested compounds were summarized in Table 5, such as water solubility, pharmacokinetics, bioavailability, drug-likeness, and medicinal chemistry. The skin absorption of the tested chemicals was reflected by skin permeability (Kp, cm/s). A more negative logKp means less skin absorption. The Kp of all the compounds is within −6.17 and −10.54, which presumes low skin permeability. Additionally, layout and absorption of drug molecules were detected using blood brain barrier (BBB) and gastrointestinal (GI) permeation. The SwissADME forecasting variables show all compounds have GI absorption, especially naringenin, hesperetin, and nobiletin. Hence, ADME prediction and physicochemical analysis suggest naringenin, hesperetin, and nobiletin may be better candidates than other tested chemicals.