RPMI-1640 (Invitrogen, Gibco, Carlsbad, CA, USA) medium with L-glutamine was used to culture the MCF7 cell line (ATCC^® HTB22™). The medium was supplemented with 10% heat-inactivated fetal bovine serum (Invitrogen/Gibco) and 1 unit penicillin/streptomycin (HyClone, Logan, UT, USA).

The 1 mM stock solution of thymoquinone (Sigma-Aldrich, Saint-Quentin-Fallavier, France) was prepared with DMSO (Sigma-Aldrich). The solution was filtered with a 0.02 μm syringe filter (Hydrophilic Ministart; Sartorius AG, Goettingen, Germany).

The MCF7 cells were seeded at 3–4×10⁶ cells/well in 96-well plates. They were cultured at 0.5% CO₂ in a humidified incubator at 37°C (Thermo Scientific, Waltham, MA, USA) for 24 h. The control cells were treated with 0.05% dimethyl sulfoxide alone (Sigma-Aldrich). Four biological replicates from each sample were prepared in separate culture flasks. The samples were treated with 50 μM thymoquinone for 24 h. EDTA [0.25% (w/v) + Trypsin/0.53% Mm] solution was used to detach the cells from the flask surface. The cells were then centrifuged at 13,000 rpm for 10 min. The supernatant was removed and the cells were washed with PBS solution prior to centrifugation at 13,000 rpm for 5 min.

RNA was isolated from the MCF7 cells using the RNeasy Plus Mini kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. The quantity of RNA was measured using a spectrophotometer (NanoDrop 2000c; Thermo Scientific). Samples with the RNA concentration (A₂₆₀/A₂₈₀ ≥1.8 ng/μl) and purity (A₂₃₀/A₂₆₀ ≥2.0 ng/μl) were selected. An Agilent 2100 Bioanalyser was used to determine the RNA integrity number (RIN). The degradation level was identified using the RNA 6000 Nano LabChip kit (Agilent Technologies, Santa Clara, CA, USA). The samples with RIN >9.8 were selected for further analysis.

The Agilent Low Input Quick Amp Labelling kit (USA) was used to generate cRNA with a sample input of 200 ng total RNA in two-color microarray analysis. The RNA Spike-In kit (Agilent Technologies) was used as an external control (positive control). It monitors and calibrates the linearity, sensitivity and accuracy of the microarray workflow. Spike A Mix with cyanine-3 (cye-3) was used to label the samples (thymoqionone-treated and untreated) and Spike B Mix with cyanine-5 (cye-5) was used to label the internal control (Universal Human Reference RNA, Agilent Technologies). T7 RNA polymerase was used to amplify target material.

The array platform used was 8×60K array SurePrint Technology (Slide Human V1 U252800417900-S01; Interscience, Rockland, MA, USA) and gasket slides with SureHyb Technology. The reference design was selected for the study. Four biological replicates from 50 μM thymoquinone-treated MCF7 and four replicates from untreated MCF7 cells were hybridised against Human Universal Reference RNA using a hybridisation kit (Agilent Technologies). The slide chamber was assembled and placed in rotisserie hybridization oven and rotated at 10 rpm in 65°C for 17 h. The array slide was washed and fluorescent imaging system was used to scan the hybridisation signals (DNA Microarray Scanner, Surescan High-Resolution Technology, Agilent Technologies).

The results were extracted with Feature Extraction Project v10.7.3.1 software and analysed using GeneSpring software v12.1. Statistical analysis was carried out using an unpaired t-test and a fold change with cut-off value >2 with a P-value ≤0.05 was considered to indicate a statistically significant difference.

qRT-PCR was carried otu to validate the results from microarray analysis and performed in two steps. The high capacity RNA-to-cDNA kit protocol (Applied Biosystems, Foster City, CA, USA) was followed to transcribe 2 μg total RNA into single-stranded cDNA. This was followed by amplification using a thermal cycler (Eppendrof, Hamburg, Germany). The ΔC_T of the thymoquinone-treated cells was subtracted from the ΔC_T of the untreated cells to determine the differences (ΔΔC_T) and fold change (2^−ΔΔCT). The human glyceraldehyde-3-phosphate dehydrogenase (GAPDH) reference gene was used as an endogenous control to normalize the fluorescence signals in the untreated and treated cells. Five dilution points were prepared in triplicate for each primer to plot the relative standard curve.

Amplification plots and the standard curve of each target gene and endogenous control were created to determine the primer efficiency. Three biological replicates of thymoquinone-treated samples and three biological replicates of untreated samples were examined. The StepOnePlus Real-Time PCR instrument (Applied Biosystems) was used to run this experiment with a 96-well system. The 5′ nuclease assay (TaqMan probe; Life Technologies Corp., Carlsbad, CA, USA) was used to detect the specific PCR product. It was labelled with the reporter fluorescence dye 6-carboxyfluorescein (6-FAM) at the 5′ end and the quencher dye MGBNFQ at the 3′ end. The fast mode was selected to run the experiment for 40 min. The temperature was adjusted to the one recommended in the manufacturer’s instructions (Applied Biosystems,).

The comparative C_T (ΔΔC_T) method was selected to determine the amount of target nucleic acid sequence in each sample relative to the untreated samples. The normalization of the ΔC_T value of the thymoquinone-treated cells to the ΔC_T value of the untreated cells was carried out to determine the ΔΔC_T value.

Filtering was carried out to minimise errors and generated 24,971 entities. There were seven entities out of 577 entities with a fold change of >2.0 and and P-value ≤0.05.

Clustering was carried out to identify entities that were grouped within a cluster which are co-regulated and functionally related. Clustering was carried out on both entities and conditions. This created a two-dimensional dendogram. The untreated samples had a similar expression profile, which linked to form a group. A similar pattern was observed in the thymoquinone-treated samples (Fig. 1).

GO describes gene products in terms of their associated cellular component and molecular function and biological process (Fig. 2). An analysis revealed eight significant GO terms related to biological processes (84%) and molecular functions (16%), as presented in Table I.

SEA identifies matching pathways within the experiment and identifies pathways related to regulated genes. An analysis of the 577 entities with a P≤0.05 and a fold change of >2.0 showed 504 pathways. There were three matched entities in cytochrome p450 (CYP450) (63 entities), phase I and II drug metabolism (49 entities) and interferon type I (54 entities).

Table II presents the downregulated pathways with the highest number of matched entities (16 entities). These were the metapathway biotransformation and interferon signaling pathway (six entities), although these were not significant (P>0.05). The most significant results were observed for the phase I metabolic pathway (P<0.04). The upregulated pathways with the highest number of matched entities were also the metapathway biotransformation (four entities). The most significant results were observed for the GPCR class B secretion pathway (P=0.009).

The most significantly upregulated genes were large intergenic non-coding RNAs (lincRNAs) and the probe was annotated with its genomic location (Table III). The protein tyrosine phosphatase, receptor type, R (PTPRR) gene was upregulated by 2.2-fold. SEA revealed genes in the estrogen pathway that were downregulated following treatment with thymoquinone (Table IV). The genes that were downregulated were UDP glucuronosyltransferase 1 family, polypeptide A8 (UGT1A8) (by −16.26-fold), cytochrome P450, family 1, subfamily A, polypeptide 1 (CYP1A1) (by −2.22-fold), cytochrome P450, family 1, subfamily B, polypeptide 1 (CYP1B1) (by −2.23-fold) and NAD(P)H dehydrogenase, quinone 1 (NQO1) (by −2.07-fold). Genes in the interferon α and β signaling pathways were also downregulated following treatment with thymoquinone; these genes were the interferon-induced protein with tetratricopeptide repeats (IFIT)1 (by −10.65-fold), IFIT2 (by −4.20-fold), IFIT3 (by −5.27-fold), interferon, α-inducible protein (IFI)6 (also known as G1P3) (-by 7.86-fold) and IFI27 (by −3.19-fold). The genes downregulated in the type II interferon signaling pathway, included interferon regulatory factor 9 (IRF9) (−2.01), IFIT2 (−4.20), 2′-5′-oligoadenylate synthetase 1, 40/46 kDa (OAS1) (−3.5) and IFI6 (−7.86).

Expansion and pathway enrichment analysis provided information on the biological functions of the entities. Network A (downregulated genes) showed the DAN domain family member 5, BMP antagonist (DAND5) as the focal point among the matched entities (Fig. 3). In addition, CYP1A1 and CYP1B1 were also noted to be highlighted in the interaction. Network B (upregulated genes) highlighted the involvement of regulator of G-protein signaling 4 (RGS4), proteasome (prosome, macropain) 26S subunit, non-ATPase, 9 (PSMD9), ubiquitin specific peptidase 9, Y-linked (USP9Y) and dystrophin (DMD).

The efficiency of each primer was determined using the relative standard curve method. The optimal concentration based on the efficiency of the primers was selected for relative quantification (RQ). The ΔΔC_T value of CYP1A1, solute carrier family 7 (anionic amino acid transporter light chain, xc-system), member 11 (SLC7A11) and UGT1A8 showed a downregulation of 43-, 15-, and 11-fold, respectively. However, the ΔΔC_T value of the protein tyrosine phosphatase, receptor type, R (PTPRR) and caveolin 1, caveolae protein, 22 kDa (CAV1) genes showed an upregulation of 2.2- and 1.5-fold, respectively.