This study, conducted in growing season of 2019, involved the collection of Papaver spp. samples from various species across the western and central regions of Iran, including Kermanshah, Tehran, Khuzestan, and Lorestan ( Table 1 ).

For genetic studies, fresh plant material was collected, identified, and preserved under controlled temperature and light conditions in an herbarium. The samples were then transferred to the laboratory, where DNA extraction was carried out following standard protocols. To assess phytochemical composition, a separate portion of each sample was immediately air-dried under shade at room temperature to minimize compound degradation. The dried samples were stored in airtight containers at -20°C until further chemical analysis. This ensured the stability of key metabolites, including alkaloids, flavonoids, and phenolics.

Ten grams of Papaver leaves were oven-dried at 70°C, followed by extraction in 100 mL of 85% methanol for 48 hours at 25°C. The extracts were then filtered, evaporated under reduced pressure, and freeze-dried. Finally, the samples were stored at -80°C until analysis for antioxidant activity.

The total phenol content was determined using the Folin-Ciocalteu method. First, 125 µL of 10% Folin reagent was added to the extracted methanolic sample, followed by incubation at 35°C for 5 minutes. Subsequently, 200 µL of 7% sodium bicarbonate solution was added, and the mixture was homogenized. Absorbance was measured at 765 nm using a UV-3200 spectrophotometer (MAPADA, Shanghai Mapada Instruments Co., Ltd., China). Results were expressed as mg of gallic acid per gram of dry tissue (Meyers et al., 2003).

The total flavonoid content (TFC) was determined using the aluminum chloride colorimetric method. First, 1.7 mL of 30% ethanol, 150 µL of 0.5 mM sodium nitrite, and 150 µL of 0.3 mM aluminum chloride were added to 150 µL of methanolic extract and homogenized. After 5 minutes, 100 µL of 1 mM sodium hydroxide solution was added. The absorbance was then measured at 510 nm using a UV-3200 spectrophotometer (MAPADA, Shanghai Mapada Instruments Co., Ltd., China). The total flavonoid content was expressed as mg quercetin equivalent (QUE) per gram of dry weight (DW) (Du et al., 2009).

The antioxidant activity of the plant was assessed using the DPPH (2,2-diphenyl-1-picrylhydrazyl) free radical scavenging assay, following the spectrophotometric method described by Brand-Williams et al. (1995). For this, 800 µL of DPPH solution (0.004 g/100 mL of 80% methanol) was added to 200 µL of methanolic extract, and the mixture was incubated in the dark for 30 minutes. The absorbance was then measured at 517 nm using a UV-3200 spectrophotometer (MAPADA, Shanghai Mapada, China). The antioxidant capacity of the extracts was expressed as DPPH inhibition percentage, calculated using the following formula:

where A₁ is the absorbance of the control (DPPH solution without extract), and A₂ is the absorbance of the sample.

Antioxidants reduce Fe³⁺ to Fe²⁺, forming a blue TPTZ–Fe²⁺ complex detected at 593 nm. For the FRAP assay, 250 µg/mL plant extract was added to 2 mL FRAP reagent (10 mM TPTZ in 40 mM HCl, 20 mM FeCl₃, 300 mM acetate buffer, pH 3.6), incubated at 37 °C for 10 min, and absorbance measured at 593 nm. A FeSO₄·7H₂O standard curve (125–1000 µM) was used to express antioxidant capacity as µmol Fe²⁺/g fresh weight (Benzie and Strain, 1996).

To quantify anthocyanin content, 1 g of dried petals was powdered and extracted using 10 mL of acidified methanol (1% HCl, v/v). The mixture was sonicated for 30 min at 20 ± 2°C, followed by storage in the dark at 4°C for 12 h with gentle stirring. After centrifugation at 10,000 rpm for 15 min, the supernatant was collected and stored at -20°C until further analysis.

Total anthocyanin content was measured using the modified pH differential method (Virgen-Ortiz et al., 2020). A portion of the extract was diluted (1:10) with 0.025 mol L⁻¹ potassium chloride buffer (pH 1.0) and 0.4 mol L⁻¹ sodium acetate buffer (pH 4.5). The absorbance was recorded at 510 nm and 700 nm at both pH levels using a UV-3200 spectrophotometer (MAPADA, Shanghai Mapada, China). The anthocyanin content was calculated using the following equation (Virgen-Ortiz et al., 2020).

Where A510 - A700 represents the absorbance difference at 510 nm and 700 nm for pH 1.0 and pH 4.5. The MW refers to the molecular weight of cyanidin-3-glucoside (449.2 g/mol), V is the final volume of the extract in milliliters, DF is the dilution factor, ϵ is the molar extinction coefficient of cyanidin-3-glucoside (26,900 L/mol·cm), L is the path length of the cuvette in centimeters, and W is the sample weight in grams. This formula allows for the calculation of the total anthocyanin content in the sample. Finally, the anthocyanin content is quantified as mg of pelargonidin-3-glucoside per 100 gr of dry weight.

The poppy petal samples were ground using a laboratory grinder (IKA, Staufen, Germany) and sieved through a 0.5 mm mesh. The resulting powdered material was then extracted to quantify the codeine and papaverine, via High-Performance Liquid Chromatography (HPLC, Agilent 1290 Infinity II). For analysis, standard solutions with concentrations of 15, 25, and 50 mg/L for both alkaloids were prepared and injected into the HPLC system. Calibration curves for codeine and papaverine were generated by plotting the area under the corresponding peaks as a function of concentration. To prepare the sample extracts, 50 mg of dry powder was mixed with 1 mL of methanol (HPLC-grade) and placed in an ultrasonic bath for 15 minutes to facilitate the dissolution of alkaloid compounds. The resulting mixture was filtered through a 0.45 µm filter before being injected into the HPLC system. The alkaloid concentrations of codeine and papaverine in the sample solutions were determined based on the area under the respective peaks and the corresponding calibration curves, with data analysis performed using the system’s software. Finally, the alkaloid content in the dry extracts was calculated based on the amount of dry material used to prepare the sample solution.

Dried leaves were treated with 70% ethyl alcohol to prevent fungal contamination prior to genomic DNA isolation (Doyle & Doyle, 1987). The samples were homogenized using a pestle, and genomic DNA was extracted from the powdered material using the DNeasy Plant Mini Kit (Qiagen, Germany). The quality and quantity of the isolated DNA were assessed via 1.5% agarose gel electrophoresis. The extracted DNA was then stored at -20°C for future use.

PCR amplification was carried out in a total reaction volume of 50 µL, containing 5 µL of Gold 10X buffer (500 mM KCl, 100 mM Tris–HCl pH = 8.3, 15 mM MgCl2, 1% Triton X-100, 1600 µg/mL BSA, 2 mM dNTP pH = 8.2–8.4), 20 pmol of forward and reverse primers, 2.5 units of AmpliTaq Gold™ DNA polymerase, sterile distilled water, and 10 ng of genomic DNA. The PCR conditions included an initial denaturation step at 95°C for 11 minutes, followed by 28 cycles of 1 minute at 95°C, 1 minute at 50°C, and 1 minute at 72°C, with a final extension at 72°C for 7 minutes. The amplification was performed using an ABI 9700 thermal cycler (Applied BioRad, USA). The ITS1 and ITS4 primer pairs were used with sequences 5’-TCCGTAGGTGAACCTGCGG-3’ and 5’-TCCTCCGCCTTATTGATATGC-3’ (Hsiao et al., 1994), respectively. PCR products were analyzed using 2% agarose gel electrophoresis for qualitative assessment.

The location of DNA fragments in the gel was determined by staining with ethidium bromide. The PCR product was loaded into the well, and the DNA migration rate was directly proportional to the applied voltage. To estimate fragment sizes, a 50-bp DNA ladder was used as a marker. After electrophoresis, the DNA fragments of interest were excised from the gel, purified, and sent for sequencing. Once the FASTA sequence file was received, the sequences were edited using the corresponding chromatogram file and then subjected to a BLASTN search. The closest matching sequences retrieved from the NCBI database were aligned and used to construct a phylogenetic tree.

Because not all species were collected across all provinces, the dataset does not constitute a fully crossed factorial design. Accordingly, each site-species combination (n=15) was treated as an independent treatment in a one-way ANOVA model. This approach allowed us to compare the means of these combinations without assuming complete species × region interactions. While LSD was used initially for pairwise comparisons, Tukey’s HSD alternatively had applied to control family-wise error in unbalanced designs. Data analysis was conducted using SAS software (PROC GLM, v. 9.4). Moreover, a Mantel test was performed to assess the correlation between genetic distance (based on molecular markers) and biochemical distance (based on phytochemical profiles) among populations. Multivariate analyses, including principal component analysis (PCA) and heat map correlation (HMC), were carried out using GraphPad Prism (v. 10.4.1). Kinship relationships were assessed using MEGA7 software, and neighbor-joining phylograms were constructed in MEGA.

Our analysis of total phenolic and flavonoid content across different Papaver species and geographic regions revealed significant variations. The results indicated that P. hybridum exhibited the highest total phenol content among the studied species, particularly in the Khuzestan region (50.26 mg GAE/g DW), whereas P. rhoeas had the lowest phenolic content (23.5 mg GAE/g DW in Khuzestan) ( Figure 1 ). For total flavonoid content, P. macrostomum exhibited the highest values across multiple regions, including Lorestan (2.66 mg QUE/g DW), Khuzestan (2.58 mg QUE/g DW), and Kermanshah (2.55 mg QUE/g DW). The lowest flavonoid content was recorded in P. rhoeas from Khuzestan (0.91 mg QUE/g DW). Overall, Papaver samples from Lorestan and Kermanshah exhibited higher levels of both phenolic and flavonoid compounds ( Figure 2 ). When analyzing species independent of region, P. hybridum consistently showed the highest phenolic and flavonoid content, while P. rhoeas contained the lowest levels across all tested populations ( Figure 3 ).

The antioxidant activity of Papaver species was assessed using both DPPH and FRAP assays. The highest DPPH radical scavenging activity was recorded in P. hybridum from Khuzestan and Kermanshah provinces, with values of 70.40% and 70.33%, respectively ( Figure 1 ). In contrast, the lowest DPPH activity was observed in P. rhoeas from Khuzestan (25.03%). Similarly, the results obtained from the FRAP assay indicated that P. hybridum from Kermanshah province exhibited the highest ferric reducing antioxidant power (6.91 µmol Fe²⁺ g⁻¹ DW), whereas the lowest FRAP activity was found in P. dubium from Tehran (0.799 µmol Fe²⁺ g⁻¹ DW) ( Figure 1 ). Overall, P. hybridum, which exhibited the highest levels of both phenolic and flavonoid compounds, also demonstrated the strongest antioxidant activity in both assays. The similarity in ranking between these two methods suggests that antioxidant mechanisms in P. hybridum involve both radical scavenging and reducing activities. Interestingly, P. rhoeas from Khuzestan, which had the lowest total phenolic and flavonoid content, also exhibited the lowest DPPH scavenging activity. Overall, Papaver samples Kermanshah exhibited higher levels of both DPPH and FRAP ( Figure 2 ). Moreover, P. hybridum demonstrated the highest DPPH and FRAP concentrations ( Figure 3 ).

As shown in Figure 1 , codeine content was highest in P. hybridum from Khuzestan (68.25 mg/g DW), while P. dubium from the same region exhibited the lowest level (14.00 mg/g DW). In contrast, papaverine was most abundant in P. dubium from Tehran (27.90 mg/g DW) and least abundant in P. bracteatum from the same province (3.32 mg/g DW). The plant samples of Papaver from Lorestan showed the greatest level of codeine and papaverine ( Figure 2 ). Overall, P. hybridum demonstrated the highest codeine concentration, whereas P. dubium from Lorestan contained the highest papaverine levels, indicating significant interspecific and regional variation in alkaloid biosynthesis ( Figure 3 ).

Anthocyanins are a subclass of flavonoids responsible for pigmentation in plants and have well-documented antioxidant properties. Among the studied samples, P. dubium from Lorestan exhibited the highest total anthocyanin content ( Figure 1 ). The plant samples from Lorestan province exhibited the higher concentration of TAC. ( Figure 2 ) The remaining Papaver species did not exhibit significant differences in anthocyanin levels ( Figure 3 ).

The PCA analysis ( Figure 4 ) revealed distinct clustering patterns among the biochemical compounds measured in Papaver species, with the first principal component (Dim1) accounting for 39.6% of the total variance and the second principal component (Dim2) explaining 20.6%. Antioxidant-related parameters, including DPPH, FRAP, total phenolic content (TPC) were closely grouped, indicating a strong positive correlation among them. In contrast, alkaloids such as codeine and papaverine were positioned separately, suggesting their biosynthesis follows an independent pathway. Total anthocyanin content (TAC) was distinctly oriented along Dim2, indicating its contribution to a different variance structure than the other measured compounds.

When the provinces were scattered based on the studied variables, the samples from Khuzestan showed a wider spread along Dim1, reflecting higher variability in their phytochemical and antioxidant profiles, while the other regions, especially Tehran, exhibited tighter clustering, indicating more similar compositions. Overall, no clear regional separation was observed in the PCA biplot, suggesting a generally high chemical similarity among the populations, with most of the variation primarily driven by Dim1.

The correlations between all biochemical variables and selected climatic parameters are shown in the Figure 5 . Relative humidity, elevation, and annual precipitation consistently exhibited positive and significant correlations with most biochemical variables, whereas mean maximum temperature (Tmax) and mean minimum temperature (Tmin) generally showed weak and mostly non-significant correlations. Among the biochemical traits, total anthocyanin content demonstrated the strongest and most significant correlations with the climatic variables.

A phylogenetic tree was constructed using the Maximum Likelihood (ML) method with a bootstrap value of 500 and a Gamma Distribution parameter (G = 5), implemented in MEGA5 software. The analysis was based on gene sequences from various Papaver species.

The resulting phylogeny ( Figure 6 ) revealed a close evolutionary relationship among Papaver orientale, P. macrostomum, P. bracteatum, and some accessions of P. armeniacum and P. dubium, all clustering within a major clade. Notably, three sequences from P. argemone formed a distinct subgroup with other accessions of P. dubium, indicating genetic affinity between these species. A few divergent sequences from P. dubium and P. orientale appeared in separate, isolated clades, suggesting intraspecific genetic variability or possible introgression events.

Interestingly, all P. hybridum sequences grouped together and were also associated with two sequences from P. orientale, suggesting a shared genetic background or gene flow between these taxa. The lack of monophyly in the Papaver and Rhoeadium clades, as inferred from the maximum parsimony analysis, further supports the hypothesis that current taxonomic classifications may not fully reflect the evolutionary history of these genera.

A subset of sequences from P. commutatum and P. rhoeas also formed a separate group, reinforcing their distinction within the genus. These findings underscore the complex evolutionary history and taxonomic challenges within Papaver and related genera.

The Mantel test revealed a weak positive correlation between genetic and biochemical distances among the studied populations (r = 0.132) ( Figure 7 ). This low correlation coefficient indicates that there is only a limited association between the patterns of genetic variability and the phytochemical composition of the samples. In other words, differences in genetic makeup do not strongly predict differences in biochemical profiles across populations. This suggests that factors beyond genetic variation—such as environmental conditions, phenotypic plasticity, or local adaptation—may have a considerable influence on the observed chemical diversity.