A Cytogenetic, Morphological, and Ecological Comparison of Seven Different Species of Achillea.ssp. Accessions in Kurdistan Province, Iran

We conducted the present study on seven important medicinal species of Achillea (4 replications per species) (in a total of 28 samples) in their natural habitats in two consecutive years (2018, 2019) in terms of morphological, cytogenetic, and ecological aspects. This study aimed to examine the environmental variables affecting the morphology, cytogenetics, and evolution of the plant. The results indicated that the populations had a Ploidy base number (x= 9) and the diploid, tetraploid, and hexaploid levels were observed. In addition to the inter-species diversity, there was the intra-species genetic diversity as (4x, 6x) Ac. millefolium (2x, 4x), Ac.vermicularis (2x, 4x), Ac. tenuifolia (2x, 4x), Ac. Alppica(2x) , Ac.talagonica(2x),Ac. biebersteinii, and Ac.wilhelmsii (4x). Further studies also indicated that 11 out of 28 populations had 1A symmetry, 15 populations had 2B symmetry, a population had 2A, and another population had 2B. Principal component analysis (PCA) of cytogenetic variables could not differentiate the species well probably due to the superiority of intra-species diversity of populations to inter-species diversity. Therefore, it seems that the evolution and speciation of this genus are mostly due to the increase or decrease in the amount of chromatin and chromosome length. The examination of principal component analysis in environmental indices also showed that Ac. millefolium hexaploid species was more adapted to the environment with higher percentages of clay and silt while the Ac.tenuifolia tetraploid species preferred a sandy habitat over other environmental factors. Furthermore, Ac.vermicularis tetraploid species indicated the greatest sensitivity to altitude. However, the Ac.biebersteinii tetraploid species reacted to meteorological parameters, such as perception rate and minimum temperature.


Introduction
Achillea genus (of family Asteraceae) is of about 130 species that are distributed from southeastern Europe to southeastern Asia and has spread to North America through Eurasia. Different species of this genus have shown signi cant adaptation to different environmental conditions and have spread from deserts and coastal areas to rocky regions. The plants of this genus are perennial and allogamous and they are spread by insects (Mozaffarian, 2005). There are 19 herbaceous species of this genus in Iran. Other species of this genus also grow in Iraq, Anatolia, Syria, Caucasus, Lebanon, Palestine, Central Russia, Transcaucasia, Turkmenistan, Afghanistan, Southwest Asia, and Central Asia in addition to Iran (Ghahreman, 1984). Yarrow is a popular medicinal herb that is widely used in traditional medicine to treat diseases, particularly burns and scars (Muẓ affarīyān, 1996).
The most common species of Achillea, which comprises 12 different species, belongs to a group of plants resulted from natural hybridization between this plant, creating subspecies that are di cult to be separated. Like other plants, polyploidy plays an important role in the evolution and speciation of this plant (as we know, about 30% to 70% of owering plants have autopolyploidy or allopolyploidy phenomenon) (Justin Ramsey & Ramsey, 2014). Polyploidy phenomena are usually found in harsh environmental conditions, such as high altitude and extreme drought. Under these conditions, these species have a higher colonization ability, greater resistance to environmental stress, and greater ecological distribution than their diploid species. Unlike allopolyploidy, autopolyploidy is known as an evolutionary dead-end over time since normally, more than two sets of homologous chromosomes are paired together during meiosis, resulting in irregular meiosis and fertility decline(J. Ramsey, 2007 Millefolium species have shown different levels of ploidies (2x, 4x, 6x and 8x) with 2A symmetrical Stebbins index (Afshari et al., 2013). Meanwhile, study of nine populations biebersteinii species demonstrated a diploid (2x); however, chromosomal interspecies variation has been observed and there has been a symmetrical Stebbins's index 1A and 2A (Chehregani Rad et al., 2017). Kai et al. (2019) reported that a high level Barley ploidy (4x) has a higher ability to dry environmental stress compared to diploid accessions. These high abilities have been detected in many biochemical pathways and morphological diversity aspects (Zhou et al., 2019).
Traditional models of allopolyploid creation predict a uniform genetic polyploid, which is due to the genetic contribution of a plant of any parent. In view of this formation along with an understandable buffering capacity of a duplicated genome, the evolutionary biologists of plants such as Stebbins and Wagner consider polyploidy as an evolutionary dead-end and a genetic uniformity that cannot respond to a changing environment(Van De Peer et al., 2017). Thus, these species will extinct as circumstances change. However, the current comments on the polyploid species sometimes show the genetic diversity among populations due to the multiple shares of parents. Furthermore, crossing between genetically distinct individuals of separate origins may create greater genetic diversity through independent sets. The result could be a set of genetically different polyploid individuals who may respond to different selection pressures in different ways and provide further opportunities for polyploid species in order to survive in a changing environment(Van De Peer et al., 2017).
As previously thought, autopolyploidy is not also a rare phenomenon. On the contrary, auto polyploidy populations could be widely dispersed in colonization in some species, for instance Biscutella laevigata, Tolmiea menziesii, and Hordeum marinum ssp. Additionally, increasing the cell size, gene diversity, and pressure resistance may make the autopolyploidy an important factor in the plant evolution (Zhou et al., 2019). However, the adaptive value of polyploidy remains unclear. A study on various species of Achillea on the coast of California indicated that hexaploid species were ve times more superior to tetraploids(Justin Ramsey, 2011;Weiss-Schneeweiss et al., 2013). The hexaploid and tetraploid species have been respectively replaced in different habitats of meadows and sand hills. These results implied that the genome replication changes features of an A.borealis species to adapt to a new environment(Justin Ramsey, 2011).
Tetraploid and hexaploid Achillea cytotypes do not coexist. Hexaploid populations mostly prefer xeric habitats over their tetraploid competitors. Since most studies have not provided the differences in the comparison of alternative cytotypes habitats, Ramsey analyzed the association of soil texture and the resident species(J. Ramsey, 2007;Justin Ramsey, 2011). Soils in tetraploid habitats have a signi cantly higher organic matter and less sand than soils in hexaploid habitats. The hexaploid species have become compatible to live in sandy hills while tetraploid species have become compatible with species of grassland and forest margins. Ramsey also indicated that a 70% superiority was observed in hexaploid species over their tetraploid competitors(Justin Ramsey & Ramsey, 2014).
The present study sought to examine seven important species of yarrow (with different ploidy levels) in their natural habitats in Sanandaj and its suburbs in terms of morphology and cytogenetics. Furthermore, efforts were made to examine the possible associations between environmental data of each habitat with the data obtained from morphological and cytogenetic analyses in order to identify the relationships of environmental parameters affecting the evolution of such species.

Plant materials
All of the 28 samples in this study, including 7 species of Achillea (Ac. millefolium, Ac.vermicularis, Ac. tenuifolia, Ac. alppica, Ac. biebersteinii, Ac.wilhelmsii, and Ac.talagonica), with 4 replicates in each species, were collected in west of Iran, Kurdistan, Sanadaj. This region is located at a longitude of 46° 59' 45" E and latitude of 35° 19' 00" N. This study was conducted in two consecutive years of 2018 and 2019. To this end, we detected 4 replicates of each species and carefully recorded the position of their geological information. To identify every species, a sample was collected from each point. It was accurately detected in the laboratory using the morphological characteristics listed in the ora of Iran (Ahmadi et al., 2013). Figure 1 represents the exact position of each location. Furthermore, we recorded the geographical position of each location using the GPS, and a soil sample was taken from a 1-35-cm pro le in each location and sent to the laboratory to estimate the pedology variables. The pedology parameters included the electrical conductivity, soil pH, soil saturation coe cient, organic carbon, total soil nitrogen, soil texture (percentage of clay, sand, and silt), soil phosphorus, and potassium. Table (1) presents the latitude and longitude of each site.

Plant morphology
All the traits continued by activating the plant in March and at the beginning of the growing season until the plant's yellowing (full maturity) and harvesting their seeds. The traits included the plant height, plant base thickness, plant stem thickness, third leaf length of in orescence, number of leaves, in orescence length, number of stems leading to ower, seed weight obtained from the harvested mass (seed yield), owering time, and nally, the amount of essential oil obtained from 40 grams of the plant. It should be noted that all the parts of the plant except for the roots were harvested in the owering stage (50%), and the essential oils were extracted utilizing a Clevenger apparatus after the plants dried in the shade. We performed all the measurements using a caliper and ruler.

Cytogenetic study of species
The seeds obtained from every point were disinfected employing the solution of the Sodium hypochlorite 2%, under sterile conditions, inside a petri dish, and on the lter paper. Afterwards, the seeds germinated at room temperature. Following two to ve days, their roots reached the proper size for sampling (roots with a length of 0.5-1 cm are appropriate for sampling). After applying the pre-treatment, the root samples were exposed to a 0.5% α-Bromonaphthalene solution for 4 hours, and running water for 30 min to remove the remains of the solution. Subsequently, we performed the xation to kill cells and maintain the internal structure of cells, and x the cell division. Thus, we used Levitsky solution as a suitable xator for karyotypic studies, and the samples were in the xation solution for 16 hours (Levitsky, 1931;Levitus et al., 2010). Following the xation, the samples were rinsed with running water for 3 hours to eliminate the residuals of the xation solution. We then used squash (softening the tissue) at an optimal level to separate the cells and put them at the same level, and make staining better. To this end, we removed the roots from 70% ethyl alcohol solution, rinsed them with running water for 30 minutes, put them in a hydrolyzer (1 N NaOH), and placed them in the oven at 60 °C for 8 minutes. After hydrolysis, the samples were dried with lter paper and placed in dye solution (hematoxylin) for 3-4 hours to stain the chromosomes (Abbaszade et al., 2017). Chromosomal images were transferred to the monitor and saved with a Digital Color CCD Camera mounted on a light microscope. The chromosomes of each cell were cut in Photoshop and arranged in a separate le. Using Micro Measure software and specifying the beginning and end of chromosomes and their centromere locations, certain characteristics such as short and long arm length, the total chromosome length, and relative chromosome length were calculated. The results were stored in Excel. In the present study, ve cells (replications) were selected and evaluated from each slide to measure chromosomal parameters. The parameters calculated for the karyotypes were as follows: CI, AR, %RL, TL,%SA,%LA, SC, A1, A2, %TF, DRL, VRC, and chromosome length range (Abbaszade et al., 2017).
According to the number of the replications, we calculated the standard deviation for the traits and the con dence interval for some of them. The chromosome form was determined using a method by Levan (Levan, 1964). After measuring the chromosomes, we drew the ideogram associated with the karyotype of the populations based on the lengths of short and long arms, in which the order of chromosomes was considered based on the length of the short arm (from large to small). We utilized the Stebbins method for comparing the karyotypic symmetry in the species (Stebbins, 1971).

Statistical data analysis
We performed all the statistical analyses employing R software. The principal component analysis (PCA) was performed separately for each environmental, morphological, and cytogenetic data series using the statistical packages, factoextra, FactoMineR, and devtools (Kassambara, 2017). We then used the results obtained from each section to interpret other sections. To extract geological and meteorological data, we extracted the data of each geographical point using a scienti c website (https://www.worldclim.org/),R software, and raster and MapTool packages(Bivand & Lewin-Koh, 2013)( Figure 1). Therefore, we obtained the data such as the minimum and maximum temperature, average temperature, perception level, radiation intensity, and slope for each point.

Results
Analysis of karyotype of the species The cytogenetic examination indicated that the accessions of the species were diploid, tetraploid, or hexaploid and had 2x, 4x, and 5x chromosomes, respectively. Therefore, we found diversity in species and also inside the accessions in rare cases in terms of ploidy levels. AL1(alppica): This sample was collected from Koodak Park in Sanandaj. The karyotype formula of the species was as 2n = 2x = 18 = 8m +1sm ( Table 2). The length of its chromosomes varied from at least 2.63 to 4.48 microns (Figure 2 for ideogram of AL1, and Figure 3). The percentage of relative length of these chromosomes was 11.11% and %TF= 39.97%, and its karyotypic type (KA) belonged to the 1A Stebbins group. The asymmetry indices, A1 and A2, had values of 0.334 and 0.17, respectively (Tables 2 and 3). AL2: This population was collected from an altitude of 2095 meters in Salavat Abad, Sanandaj. Its karyotypic formula was as 2n = 2x = 18 = 9m and all the chromosomes were metacentric (Figure 2 for L2 ideogram, and Figure 3). Base on Figure 4, the lengths of its chromosomes ranged from 2.53-4.39 microns. The relative percentage of chromosomes (%RL) was about 11.11, %TF= 39.78%, and its karyotypic type belonged to the 1A Stebbins group. The A1 and A2 asymmetry indices were 0.335 and 0.162, respectively (Tables 2 and 3).
AL3: The seeds of this population were harvested from an altitude of 1829 meters in Hasanabad of Sanandaj. Its karyotypic formula was as 2n=2x=18= 6m+3sm ( Table 2). The lengths of its chromosomes ( Figure 2 for ideogram of AL3, and Figure 3) were from 3.18-5.39 microns. %RL= 11.11% and %TF=37.6%, and it was in class 1A Stebbins in terms of karyotypic symmetry. The asymmetry A1 and A2 were equal to 0.395 and 0.169, respectively (Tables 2 and 3).
AL4: This sample was obtained from an altitude of 1607 meters in Goyran village. Its karyotypic formula was as 2n= 2x= 18= 8m+1sm ( Table 2). The variation of chromosomes ranged from a minimum of 2.58 to a maximum of 4.33 ( Figure 2 for ideogram of AL4, and Figure 3). %RL=11.11% and %TF=39.27, and it was in 1A Stebbins group in terms of karyotypic symmetry. The asymmetry A1 and A2 indices were equal to 0.351 and 0.166, respectively (Tables 2 and 3). BI1(biebersteinii): This sample was collected from an altitude of 1334 meters in Savarian village. Its karyotypic formula was as 2n=4x= 36= 14m+ 4sm ( Table 2). The variation in lengths of its chromosomes ranged from 1.95-4.81 microns ( Figure 2 for BI1 ideogram, and Figure 3). %RL= 5.59, %TF= 38.42%, and it was put in 1B class of Stebbins in terms of karyotypic symmetry. The intra-chromosomal asymmetry A1 and inter-chromosomal asymmetry A2 were respectively equal to 0.356 and 0.228 (Tables 2 and 3).
BI2: This sample was collected from an altitude of 1436 meters in Pichon village, Sanandaj. The sample had a karyotype formula of 2n=4x=36=15m+3sm. The total length of the chromosome from a minimum to a maximum was 1.59-3.79 microns ( Figure 2 for the ideogram of BI2, and Figure 3). %RL=5.5, %TF=39.43, and it was in class 1B of Stebbins in terms of karyotypic symmetry. The asymmetry A1 and A2 indices were equal to 0.343 and 0.19, respectively (Tables 2 and 3).
BI3: This sample was collected from an altitude of 1436 meters in Danikesh village. Its karyotypic formula was as 2n=2x=4x=36=14m+4sm ( Table 2). The variation of chromosome length ranged from 2.53-5.38 microns ( Figure 2 for the ideogram of BI3, and Figure 3). %RL= 5.56, %TF=38.44, and it was put in Stebbins 1B group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were equal to 0.376 and 0.195, respectively (Tables 2 and 3). BI4: This sample was collected from an altitude of 1593 meters in Chehel Gazi village. Its karyotypic formula was 2n=4x=36=14m+4sm ( Table 2). The variation of its chromosome length ranged from 2.01-4.35 microns ( Figure 2 for the ideogram of BI4, and Figure 3). %RL and %TF were 56.5 and 38.62, respectively. In terms of symmetry, it was put in the 2B group of Stebbins. Other characteristics of the population were the intra-and inter-chromosomal indices, respectively equal to 0.343 and 0.217 (Tables 2   and 3). MI1(millefolium): This sample was collected from an altitude of 1985 meters in Salavat Abad village. Its karyotypic formula was 2n=6x=54=20+7sm ( Table 2). The variation of its chromosome length ranged from 2.22-4.29 microns ( Figure 2 for the ideogram of MI1, and Figure 3). %RL and %TF were 3.7 and 38.8, respectively. In terms of symmetry, it was put in the 2A group of Stebbins. The asymmetry A1 and A2 indices were equal to 0.355 and 0.153, respectively (Tables 2 and 3).
MI2: The second sample on this sub-population was collected from an altitude of 1498 meters in Baba Riz village of Sanandaj. Its karyotypic formula was as 2n=6x=54=17m+10sm ( Table 2). The variation of chromosome length ranged from 2.02-4.24 microns ( Figure 2 for the ideogram of MI2, and Figure 3). %RL= 3.7, %TF=37.89, and it was put in Stebbins 1B group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were equal to 0.378 and 0.167, respectively (Tables 2 and 3). MI3: The third sample on this sub-population was collected from an altitude of 1993 meters in Dolbandi village of Sanandaj. Its karyotypic formula was as 2n=4x=36=15m+3sm ( Table 2). The variation of chromosome length ranged from 2.15-4.18 microns ( Figure 2 for the ideogram of MI3, and Figure 3). %RL= 5.56, %TF=39.29, and it was put in the Stebbins 1A group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were respectively equal to 0.337 and 0.183 (Tables 2 and 3). MI4: The fourth sample on this sub-population was collected from an altitude of 1642 meters in Jebreilian village of Sanandaj. Its karyotypic formula was as 2n=6x=54=20m+7sm ( Table 2). The variation of chromosome length ranged from 2.24-5.12 microns ( Figure 2 for the ideogram of MI4, and Figure 3). %RL= 3.7, %TF=27.85, and it was put in Stebbins 1B group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were equal to 0.375 and 0.194, respectively (Tables 2 and 3).

TA1(talagonica):
The rst sub-sample of this species was collected from an altitude of 1993 meters in Charandu village 20 km outside of Sanandaj. Its karyotypic formula was as 2n=2x=18=9m ( Table 2). The variation of chromosome length ranged from 2.24-5.12 microns ( Figure 2 for the ideogram of TA1, and Figure 3). %RL= 11.11, %TF=42.02, and it was put in Stebbins 1B group in terms of karyotypic symmetry.
Asymmetry A1 and A2 indices were equal to 0.277 and 0.236, respectively (Tables 2 and 3). TA2: Another sub-sample of this species was collected from an altitude of 1980 meters in Baynchub village 54 km outside of Sanandaj. Its karyotypic formula was as 2n=2x=18=7m+2sm ( Table 2). The variation of chromosome length ranged from 2.34-5.14 microns ( Figure 2 for the ideogram of TA2, and Figure 3). %RL= 11.11, %TF=40.37, and it was put in Stebbins 1B group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were respectively equal to 0.224 and 0.311 (Tables 2 and 3). TA3: This sub-sample was collected from an altitude of 1919 meters in Mamukh-e So a village. Its karyotypic formula was as 2n=2x=18=7m+2sm ( Table 2). The variation of chromosome length ranged from 2.34-5.14 microns ( Figure 2 for the ideogram of TA3, and Figure 3). %RL= 11.11, %TF=40.68, and it was put in Stebbins 1B group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were equal to 0.315 and 0.223, respectively (Tables 2 and 3). TA4: The fourth sample of this group was collected from an altitude of 1577 meters in Sarab Qamish village. Its karyotypic formula was as 2n=2x=18=6m+3sm ( Table 2). The variation of chromosome length ranged from 1.8-3.54 microns ( Figure 2 for the ideogram of TA4, and Figure 3). %RL= 11.11, %TF=38.86, and it was put in the Stebbins 1A group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were equal to 0.363 and 0.209, respectively (Tables 2 and 3).

VE1(vermicularis):
The rst sub-population of this group was collected from an altitude of 2334 meters in Qalvaz village. Its karyotypic formula was as 2n=2x=18=9m ( Table 2). The variation of chromosome length ranged from 4.92-9.85 microns ( Figure 2 for the ideogram of VE1, and Figure 3). %RL= 11.11, %TF=39.74, and it was put in the Stebbins 1A group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were equal to 0.339 and 0.21, respectively (Tables 2 and 3). VE2: The second sub-population of this group was collected from an altitude of 2145 meters in Mamukh-e So a village. Its karyotypic formula was as 2n=4x=15m+3sm ( Table 2). The variation of chromosome length ranged from 2.35-3.29 microns ( Figure 2 for the ideogram of VE2, and Figure 3). %RL= 5.56, %TF=39.84, and it was put in Stebbins 1B group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were equal to 0.347 and 0.223, respectively (Tables 2 and 3). VE3: The third sub-population of this group was collected from an altitude of 2152 meters in Sangse d village. Its karyotypic formula was as 2n=4x=36=17m+1sm (Table 2)  VE4: The fourth sub-population of this group was collected from an altitude of 1838 meters in Dul rahman village. Its karyotypic formula was as 2n=4x=36=17m+1sm ( Table 2). The variation of chromosome length ranged from 2.36-5.3 microns ( Figure 2 for the ideogram of VE4, and Figure 3). %RL= 5.56, %TF=40.52, and it was put in Stebbins 1B group in terms of karyotypic symmetry. Asymmetry A1 and A2 indices were respectively equal to 0.317 and 0.218 (Tables 2 and 3). TE1(tenifolia): This sample was taken from Mamukh mountain pass and an altitude of 1886 meters. Its karyotypic formula was as 2n=4x=36=18m and it was put in the Stebbins 1A group in terms of karyotypic symmetry. The variation of chromosome length ranged from 1.94-3.37 microns ( Figure 2 for the ideogram of TE1, and Figure 3). %RL= 5.56 and %TF=42.5. Asymmetry A1 and A2 indices were equal to 0.257 and 0.159, respectively (Tables 2 and 3). TE2: This sample was taken from Charandu village and an altitude of 1841 meters. Its karyotypic formula was as 2n=4x=36=14m and it was put in Stebbins 1B group in terms of karyotypic symmetry ( Table 2). The variation of chromosome length ranged from 11.82-4.1 microns ( Figure 2 for the ideogram of TE2, and Figure 3). %RL= 5.41 and %TF=39.36. Asymmetry A1 and A2 indices were equal to 0.347 and 0.223, respectively (Tables 2 and 3).
TE3: This sample was taken from Bazi Rabab village and an altitude of 1934 meters. Its karyotypic formula was as 2n=4x=36=11m+7sm and it was put in Stebbins 1B group in terms of karyotypic symmetry ( Table 2). The variation of chromosome length ranged from 1.6-5.93 microns ( Figure 2 for the ideogram of TE3, and Figure 3). %RL= 5.39 and %TF=38.56. Asymmetry A1 and A2 indices were equal to 0.369 and 0.199, respectively (Tables 2 and 3).
TE4: This sample was taken from Sarab Qamish village and an altitude of 1577 meters. Its karyotypic formula was as 2n=4x=36=17m+7sm and it was put in Stebbins 1B group in terms of karyotypic symmetry ( Table 2). The variation of chromosome length ranged from 2.16-3.43 microns ( Figure 2 for the ideogram of TE3, and Figure 3). %RL= 11.11 and %TF=46.11. Asymmetry A1 and A2 indices were respectively equal to 0.141 and 0.147 (Tables 2 and 3). Wl1(wilhelmsii): This sample was taken from Mamukh-e So a village and an altitude of 1999 meters. Its karyotypic formula was as 2n=4x=36=15m+3sm and it was put in Stebbins 1B group in terms of karyotypic symmetry ( Table 2). The variation of chromosome length ranged from 2.25-4.85 microns ( Figure 2 for the ideogram of WI1, and Figure 3). %RL=5.56 and %TF=39. Asymmetry A1 and A2 indices were equal to 0.346 and 0.19, respectively (Tables 2 and 3).
Wl2: The second sample of this species was taken from Arandan village and an altitude of 1974 meters. Its karyotypic formula was as 2n=4x=36=14m+4sm and it was put in Stebbins 1B group in terms of karyotypic symmetry ( Table 2). The variation of chromosome length ranged from 2.1-4.65 microns ( Figure   2 for the ideogram of WI2, and Figure 3). %RL=5.56 and %TF=39.13. Asymmetry A1 and A2 indices were equal to 0.371 and 0.205, respectively (Tables 2 and 3).
Wl3: The third sample of this species was taken from Klatei village and an altitude of 2208 meters. Its karyotypic formula was as 2n=4x=36=16m+2sm and it was put in the Stebbins 1A group in terms of karyotypic symmetry ( Table 2). The variation of chromosome length ranged from 2.46-4.84 microns ( Figure 2 for the ideogram of WI3, and Figure 3). %RL=5.56 and %TF=39.13. Asymmetry A1 and A2 indices were equal to 0.342 and 0.168, respectively (Tables 2 and 3).
Wl4: Ultimately, the fourth sample of this species was taken from Gavdarreh village and an altitude of 2026 meters. Its karyotypic formula was as 2n=4x=36=14m+4sm and it was put in the Stebbins 1A group in terms of karyotypic symmetry ( Table 2). The variation of chromosome length ranged from 2.45-4.84 microns ( Figure 2 for the ideogram of WI4, and Figure 3). %RL=5.56 and %TF=37.75. Asymmetry A1 and A2 indices were equal to 0.38 and 0.183, respectively (Tables 2 and 3).

Comparison of cytogenetic parameters between species
The comparison of the results of karyotypic characteristics in the populations indicated that the base chromosome number was x=9 in all the populations and there were hexa-, tetra-and diploid levels for the populations. Regarding the ploidy level, there was diversity not only among the species, but also among the populations of the three species, Ac. millefolium (tetra-and hexa-ploidy), Ac.vermicularis (tetra and diploidy), and Ac.tenuifolia (tetra and diploidy). Ac. alpine and Ac. talagonica species were diploid and Ac.biebersteinii and Ac.willhelmsii species were tetraploid. Karyotype formulas of inter-species and intraspecies populations were different and all the chromosomes were metacentric only in populations AL2, TA1, VE1, TE1, and TE4; the karyotype consisted of a large number of metacentric chromosomes and a small number of chromosomes were submetacentric in other populations. According to the Stebbins bilateral table, most of the populations were in Classes 1A and 1B; only population BI4 was in class 2B, and MI1 in the class 2A. Therefore, a symmetrical karyotype was observed for the species of this genus. The highest relative amount of chromatin belonged to population AL3 with an average of 4.15 microns whereas the lowest relative amount of chromatin belonged to population TA2 with an average of 2.55. Except for population AL3, the relative chromatin levels of the populations were less than 4 and more than 2 microns. Since the relative difference in the lengths of chromosomes had an inverse relationship with intra-species ploidy levels, the most asymmetric chromosomes among the hexaploid populations, based on DRL index, belonged to population MI4 with an average of 38.3% (Tables 2 and 3). For diploid populations, BI1 population had the highest rate of chromosomal asymmetry with the highest DRL (5.32 percent). Among the diploid populations, VE1 population and four populations of Ac. talagonica species had the highest rate of DRL, and the most asymmetric chromosomes. The percentage of overall chromosomes form ranged from 37.6 to 46.11, and the highest percentage of overall chromosome form belonged to populations TE4, TE1, and TA1; thus, they had a more symmetrical karyotype compared to the other populations. On the contrary, AL3, MI2, MI4, and WI4 had the lowest percentage of overall chromosome form; therefore, they had the most asymmetric karyotypes (Tables 2 and 3). The lowest intra-chromosomal asymmetry index belonged to populations TE4, TE1, and TA1; consequently, they had more symmetrical karyotypes than the other populations. Based on index A1, AL3 and Wl4 had the highest chromosomal asymmetry. Hence, it was found that the intra-species diversity was high for A1 and TF%, and the species were indistinguishable based on the parameters. For the inter-chromosomal asymmetry, the intra-species diversity was somewhat lower, and the species could be divided into three categories; the species of the rst class included Ac.alppica and Ac.millefolium, whose populations had an inter-chromosomal index of less than 0.2 and symmetrical chromosomes based on the index. On the contrary, the populations of two species, Ac. talagonica and Ac.vermicularis, had an index A2 of over 0.2 and asymmetric chromosomes based on the index. However, the populations of other species had higher intra-species diversity compared to the abovementioned four species populations and also had populations with low inter-chromosomal asymmetry and high inter-chromosomal asymmetry ( Table 2).
In terms of the centromeric index (Table 3), populations TE4, MI1, TE1, VE3, and TA1 had a centromeric index between 0.42 and 0.46 and they had symmetrical chromosomes based on the index. Meanwhile, populations AL3, BI3, MI2, MI4, WI2, and WI4 with a centromeric index of 0.38 had the most asymmetric chromosomes based on the index. The lowest ratio of long to short arm belonged to TE4, TE1, and TA1 populations with average values of 1.19, 1.37, and 1.4, respectively, and had symmetrical chromosomes based on the index. On the contrary, the highest value for the index with long to short arm ratio between 1.6 and 1.6 belonged to populations AL3, BI3, MI2, MI4, WI4, TE3, and WI2, and thus, they had asymmetric chromosomes. The highest average total chromosome length belonged to population AL3 with an average of 4.15 microns, and other populations had an average total chromosome length between 2.55 and 3.64, among which populations TA2, TE1, TA4, TE4, BI2, VE2, MI1, VE1, BI4, and BI1 had an average total chromosome length of less than 3 microns. The other populations had an average total chromosome length between 3 and 3.64 microns. Therefore, it was found that intra-species diversity was high for AR, CI, and TL indices ( Table 3). Based on the parameters, the species were indistinguishable. The range of the total chromosome length varied widely from a minimum range of 1.59 microns in the population BI2 to a maximum of 9.85 microns in the population VE1; hence, the longest chromosome was 6.19 times higher than the shortest chromosome.
Environmental data analysis  (Fig. 4b). In the rst zone for the coordinates of components, SP, Silt, and Clay variables indicated the highest impact (the variables such as soil carbon and potassium were in the second degree of importance. In the second zone, there were two variables, namely soil phosphorus and regional altitude, among which the altitude was much more involved. In the third zone, two variables, namely the regional slope and the soil sand amount, were more involved in explaining the variance. Finally, the parameters such as the average temperature, minimum temperature, maximum temperature, perception, soil pH, and soil nitrogen levels were put in the fourth zone. Furthermore, meteorological variables are often located in this zone. It should be noted that the components, namely average temperature, minimum temperature, and precipitation rate have greater effects than the other parameters. Figure 4c represents the involvement of variables on the rst component (28.6%). As shown in the gure, variables such as Tmin, perception, and soil sand level had the greatest effect on the variance of the rst component. Figure 4d shows the involvement of variables on the second component. SP, soil clay amount, sand amount, carbon, potassium, soil pH, and light intensity variables had the greatest impact on the variance of the second component. As mentioned earlier, these variables are often in the rst coordinate region. Figure 4c, d illustrates that the variables such as SP, Tmin, Temp, alts, and amounts of sand, clay, and silt, nitrogen, Tmax, and amount of precipitation play roles in explaining the variance of the rst and second components. Furthermore, the variables had the highest environmental impact on the distribution of the species. Figure 4e shows the distribution of seven species with 4 replications in the zone. It explains 28.6% in the rst dimension, 25.8% in the second dimension, and 54% of the total variance. In the rst zone, there were three replications of species MI2, MI3, and MI4 (blue dotted line), and species TE2, Wl2, and VE4. There was Millifolume species in this zone. In the second zone, there were species VE (1, 2, 3), WI (1,4), and MI1. In other words, these three species had a more prominent presence in this zone. The size of each point perfectly indicated the degree of species participation in explaining the total variance of the components. The more the color inclined to red and the larger the size of each dot was, the greater its participation in the data variance would be. In the third zone, there were species TE (1, 2, 3, 4), AL (2, 3), TA3, and WI3. There were also species BI (1, 2, 3, 4), AL (1,4), and TA1 in the fourth zone. In other words, the presence of two species, BI and AL, was more prominent in this zone. Figure 5a depicts the results of PCA for cytogenetic variables of different species. A total of 11 variables were studied. According to the gure, the rst two components explained 74.8% of the total variance ( Figure 5b). The components highlighted in red had a very high contribution to the variance explanation, and as the color intensity tends to be turquoise, the variable contribution to the variance explanation decreased.

Results of analysis of cytogenetic variables
In the rst zone, there were SA, VRC, TL, and LA variables with almost equal weights probably due to the alignment of these variables with each other. These variables essentially had the same nature. In the second zone, there were A2, DRL, %TF, and CI variables among which %TF had a higher contribution rate compared to A2 and DRL. There were A1, AR, and X variables in the fourth zone.   3), and BI3 species whose size and color represent their contribution to variance. There were TA3, TA1, VE3, and TE4 in the second zone, and also TE (1, 2), BI2, MI3 VE (1, 2), and TA (2, 4) species in the third zone. Eventually, there were BI (1, 4), MI (1,2,4), and WI (2,4) in the fourth zone. Three groups, AL, MI, and WI were separated from the other species in the analysis, and the cytogenetic variables could not separate other species from each other. Figure 6a illustrates the results of PCA on the morphological variables of different species in the rst year of the study. The data analysis indicated that the rst component explained 36.4%, the second component 18.6%, and a total of 55% of the whole variance is explained (Figure 6b). The owering time, SD, BD, PHT, LTLin , and x variables were in the rst zone. There were two variables, StkinFlo and NLeaf in the second zone. Finally, there were two variables, seed yield, and essential oil yield of the species in the fourth zone. Figure 6c shows the degree of contribution of the variables in the explanation of the rst component. Five variables, BD, SD, PHT, LTLin , and NLeaf had the largest contribution to the explanation of the variance of the rst component. In Figure 6d, four variables, namely stkin o, LinFlo, Yld, and Nleaf had the largest contribution to the explanation of the variance of the second component. Figure 6e demonstrates the results of the species based on morphological data in the rst year. Species TE (1,2,3,4) was in the rst zone and was well separated from the other species. In the second zone, two species, VE (1,2,3,4) and AL (1,2,3,4) were completely different from each other and the other species. In the third zone, there were other two species, WI and TA; nevertheless, the dispersion rate of WI was higher than that of species BI. Finally, there were two species, MI and BI, in the fourth zone. As could be seen, the morphological data could clearly distinguish all the species from each other.   No speci c trends were observed in the rst zone. In the second zone, two species, AL and VE, were clearly distinguished from other species and each other. The third zone included two species, WI and TA, which were completely on the left bottom corner of the diagram. Ultimately, there were three species in the fourth zone. MI species was almost in the zone close to the axis, and other species were away from it, on the right top of the graph. These species included BI and TE.

Discussion
There are high diversity and differences in chromosomal length characteristics of the inter and intraspecies of this genus. Given that the existence of diversity and difference in chromosome length indicates an advanced karyotype and has chromosomes in different sizes (Afshari et al., 2013) , the species of this genus have advanced karyotypes. The existence of x=9 as the base chromosome number on the yarrow genus has been proven in several reports, yet the number of chromosomes and ploidy levels vary among different species of this genus, which could range from 2n=2x= 18 to 2n=8x=72 even though most species are diploid (Baltisberger & Widmer, 2016;Guo et al., 2005) . In addition to inter-species diversity in ploidy levels, there are numerous reports of ploidy level diversity in populations within a species. In other words, different ploidy levels are reported for populations of a species (Ebrahim et al., 2012; Hoshi et al., 2010) . Accordingly, a range between diploid to hexaploid has been reported for Ac.aleppica species (Chehregani & Javaheri, 2014); however, all the accessions of the species were diploid in the present study. The tetraploid level was reported for A. bieberestini species (Afshari et al., 2013), which was consistent with the present results Afshari and et al., reported diploid and tetraploid levels for Ac.millefolium species (Afshari et al., 2013). In another study, Hexa and Octa-ploidy levels were reported for the species (Ebrahim et al., 2012). The two reports were consistent with the present study in terms of Ac.millefolium species. For four populations of Ac. talagonica species, the diploid level was in accordance with results of studies by Sahin et al. (Sahin et al., 2006). Finally, the results obtained for ploidy levels of two species, Ac. vermicularis and Ac. Tenuifolia, were in agreement with other reports ( (Table 3); hence, there were more symmetrical chromosomes in the present research than other reports since the karyotype A2 was mostly reported in other reports, and fewer cases had A1 and B1 symmetries (Kiran et al., 2012;Sahin et al., 2006). Accordingly, no obvious differences were reported in karyotype asymmetry between yarrow species; all the species had symmetrical karyotype structures because most chromosomes were metacentric and sub-metacentric (Kiran et al., 2012).
Satellites were observed more in populations with tetra-and hexaploid levels and on chromosome 1 ( Figure  2). No satellites were observed in diploid populations, and there was only a satellite for each population. Our results were consistent with those of a report by Sahin et al. (Sahin et al., 2006). On the contrary, no satellites were reported in certain studies (Afshari et  There were chromosomes B in two populations of Ac. tenuifoli species (Table 2). A chromosome B was also reported for the species in some populations (Chehregani Rad et al., 2017), and there were some reports on the existence of chromosome B in other species on the genus (Baltisberger & Widmer, 2016). Nevertheless, there was no chromosome B in some reports (Kiran et al., 2012).
There was no inter-species diversity for the chromatin content, arm length, and chromosome length.
Furthermore, the intra-species populations showed more diversity (Table 3), but there was inter-species diversity for ratios to arms (centromeric index and large to small arm ratio); however, the intra-species populations had diversity. Therefore, the evolution and speciation of the genus was through intra-chromosomal asymmetry rather than increasing or decreasing the chromatin content and chromosome length. The average length of each chromosome ranged from 2.93 to 3.55 microns for the species, which was consistent with other reports (Afshari et al., 2013;Sahin et al., 2006). Meanwhile, the chromosome length range was higher in certain reports than that in the results of the present study, and longer chromosomes were reported for the species (Aksu et al., 2013). Based on the karyotypic characteristics, the Ac.alppicaIn species had more karyotypic evolution in terms of chromatin content, and three species, Ac.biebersteinii, Ac.wilhemsii, and Ac.millefolium, had more complete karyotypes due to the intrachromosomal asymmetry, and a higher evolution in terms of chromosome length characteristics and chromatin content. Ac.talangonica, Ac.tenifolia, and Ac.vermicularis had karyotypic evolution due to the chromosomal asymmetry; thus, Ac.biebersteinii, Ac.wilhemsii, and Ac.millefolium had more evolved Karyotypes than the other species. According to the results, the karyotypic characteristics could not separate the populations of Yarrow species due to the intra-species diversity, and the populations of different species were in the same group in several cases. Despite the lack of comprehensive reports on the study of inter and intra-species relationships for the yarrow genus based on karyotypic characteristics, the few available reports indicated that the karyotypic characteristics were unable to completely separate the According to the results, the seven species had ploidy levels of 2x, 4x, and 6x with the intra-species diversity in some samples. The soil and ecological data also indicated the high diversity in the natural habitat of the species. Even though the Achillea genus performs pollination through the allogamy by insects, its cytotypes are located in the far distance in nature. The analysis of environmental indices has shown that the species could be largely distinguished. Millifolume species (with a hexa ploidy level) was in the rst zone of the PCA on the environmental data, indicating that the environmental indices, such as soil permeability, silt percentage, and soil clay, are more critical for the growth of this species(Justin Ramsey, 2011;Weiss-Schneeweiss et al., 2013). In other words, this species has become more adapted to these environmental conditions and evolved for this purpose. In the second zone, there is the environmental parameter, height, and it has the vermicularis tetraploid, indicating the response of species to the habitat height. In the third zone, there is the soil sand percentage component with which the tetraploid species, Ac.tenuifolia, is consistent. In other words, the Ac.tenuifolia species prefers the sandy habitat over other environmental factors(Justin Ramsey, 2011). In the fourth zone, there are meteorological variables, such as minimum and medium temperatures, and precipitation against which the tetraploid species, Ac.biebersteinii, has had the most reaction (Figures 4a and 4e). This might answer the question of why cytotypes of this genus cannot be found together in one place. It seems that different evolutionary processes, autopolyploidy of each species for instance, have led to special habitats and thus, they cannot be found in one place(Justin Ramsey & Ramsey, 2014).
The results of the PCA on cytogenetic indices (Figure 5e) can only separate the three species from other species. For instance, the diploid species, Ac. Alppica (marked with a light green label) in the rst zone is consistent with indices such as SA, VRC, TL, and LA. In the fourth zone, the hexaploid species, Millifolume (marked with a turquoise label), is more prominent with two indices, A1 and AR, compared to the rest of the species. The tetraploid species, Ac.wilhelmsii (marked with the mustard label), is located exactly at the interface between these two groups; hence, it has the intermediate characteristics of the two species. The cytogenetic indices could not appropriately isolate the species so that they would have insu cient diversity to distinguish species from each other. It is noteworthy that the diversity of the intra-species cytogenetic indices, as previously mentioned, is probably a reason why species cannot be separated with the help of the principal component analysis.
On the contrary, the PCA results for both years of morphological data could perfectly distinguish all different species, and almost the same trend is followed in both years. Many morphological variables are located in the rst zone of the coordinate axis, and there are few indices in other zones. The morphological parameters could well distinguish the species in the two years probably due to the high diversity of the indices among the species.

Conclusion
The highest mode of speciation in the next 500 years is polyploidy (auto or allo), which considering climate change, leads to an increase in plant polyploidy from 35% to 50%. In addition to its complexity, there are several variables involved in the creation and evolution of these species. Identifying these effective factors requires further research from different aspects and perspectives. These factors differ from one genus to another, and the intra-species determinants may be also different. This could pave the way for the evolution of a new species. The present study aimed to identify the important parameters affecting yarrow in terms of oil compounds, including the pedology variables (clay, sand, and penetration), meteorological factors (minimum temperature and precipitation), and altitude, which were consistent with the natural habitats of the species and triggered their proliferation, adaptation, and separation at ecological niche layers. However, the analysis models of the present study, the PCA for instance, could not explain more than 54% of the total variance of the environmental variables, indicating that we were far from identifying the factors (probably due to other unknown parameters not considered in the present study). Furthermore, the cytogenetic indices could not separate different species in the present study. Therefore, there is a need for precise genomic indices n top of reviewing the phenotyping methods and recording the environmental data to enter further parameters in the model and pursue the polyploidy codes of plants, particularly Yarrow and its factors, with a novel perspective.  Table 3. Karyotypes and chromosomal parameters in this study for each species.
Positions samples were collected on the map. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 3
The morphological diversity of chromosomes between and within species in the metaphase stage. which explained about 54.4% of the total variance. The share of each species is determined by the intensity of the color; the more reddish the color is, the higher the share would be and vice versa. Drawn ellipses distinguish the groups (text description).