ARTÍCULO

Microsatellite data reveal genetic restructuring of Medicago sinskiae (Fabaceae) in western and southwestern Iran

RAHA ZAREEI¹, ERNEST SMALL², MOSTAFA ASSADI³ & IRAJ MEHREGAN¹

¹ Department of Biology, Science and Research Branch, Islamic Azad University, IR-1477893855 Tehran, Iran
² Agriculture and Agri-Food Canada, Ottawa Research and Development Centre, K.W. Neatby Building, 960 Carling Ave., CA-K1A 0C6 Ottawa, Canada
³ Research Institute of Forests and Rangelands, Agricultural Research Education and Extension Organization (AREEO), IR-13185116 Tehran, Iran

ORCID iD. R. ZAREEI: https://orcid.org/0000-0001-6496-6216, E. SMALL: https://orcid.org/0000-0002-0654-6966, M. ASSADI: https://orcid.org/0000-0001-9014-585X, I. MEHREGAN: https://orcid.org/0000-0002-5108-2558

Author for correspondence: I. Mehregan (iraj@daad-alumni.de; imehregan@srbiau.ac.ir)

Editor: J. López-Alvarado

ABSTRACT
Microsatellite data reveal genetic restructuring of Medicago sinskiae (Fabaceae) in western and southwestern Iran.— Medicago sinskiae appears to be a very rare species in the Iranian flora with only a few records in the last three decades. Eight populations (62 individuals) of M. sinskiae, one population of M. rigidula (seven individuals), and one population of M. constricta (five individuals) from western and southwestern Iran were analyzed for microsatellite data based on newly designed SSR primers using NGS technology. The PCoA, Clustering and Structure analyses showed no geographical pattern of genetically designated clusters. Our results showed that M. sinskiae is mainly an inbreeder. It is assumed that high levels of gene flow (Nm) and generation of genetically homogenous populations seem to be more affected by fast dispersal and not localized gene flow. Extensive collections recently made from the western and southwestern Iran showed that its presence is increasing. Finally, our results indicate that the species is segregated from its very close relatives M. rigidula and M. constricta in Iran.
KEY WORDS: annual medics; Medicago; microsatellites; population genetics.

RESUMEN
Los marcadores microsatélite revelan la reestructuración genética de Medicago sinskiae (Fabaceae) en el oeste y el sudoeste de Irán.— Medicago sinskiae es considerada una especie rara en la flora iraní con únicamente unas pocas citas en las tres últimas décadas. Se han muestreado ocho poblaciones (62 individuos) de M. sinskiae, una población de M. rigidula (siete individuos) y una población de M. constricta (cinco individuos) en el oeste y el suroeste de Irán que han sido analizadas con marcadores microsatélite. Se han utilizado nuevos primers obtenidos con tecnología NGS. Los análisis de PCoA, Clustering y Structure no muestran un patrón geográfico para los clústeres genéticos. Los resultados muestran que M. sinskiae es principalmente una especie autógama. Se asume que los altos niveles de flujo genético (Nm) y la homogeneidad genética poblacional están afectados por una rápida dispersión y un flujo genético no localizado. Recolecciones extensivas realizadas recientemente en el oeste y el suroeste de Iran muestran que el rango de distribución esta especie se está incrementando. Finalmente, nuestros resultados indican que M. sinskiae está diferenciada de las especies M. rigidula y M. constricta en Irán.
PALABRAS CLAVE: genética de poblaciones; marcadores microsatélite; Medicago anuales.

CONTENTS

INTRODUCTIONTop

The genus Medicago L. (Fabaceae) has more than 80 annual or perennial species mainly distributed around the Mediterranean Basin (Small, 2011Small, E. 2011. Alfalfa and relatives. Evolution and classification of Medicago. NRC Research Press, Ottawa.). Medicago sinskiae Uljanova, the subject of this study, is a very poorly known annual herb, with the following key characteristics. The mainly simple-hair pubescent stems are usually 15–25 cm (rarely up to 40 cm) long, prostrate to ascending, branched from base. The stipules are 2–4 mm, rarely 5 mm long, dentate to laciniate. Each peduncle bears 1–3 (rarely up to 5) flowers. The flowers are 3–6 mm in length, with pubescent calyx and yellow or orange-yellow corolla. The mature pods are ovoid, cylindrical, or discoid, pubescent with both simple and gland-tipped hairs, 4–8 mm long, 4–6 mm wide, with 2.5–5 coils, spineless or with spines up to 4 mm long, hardened at maturity with some gaps often present between coils (Fig. 1). Seeds are 1.75–2.5 mm long, 1–2 per coil, separated by spongy fruit partitions, smooth in surface, yellow to yellow-brown in color, with the radicle about half as long as the length of the seed. Flowering starts in early April and the fruits mature in May and June (Small, 2011Small, E. 2011. Alfalfa and relatives. Evolution and classification of Medicago. NRC Research Press, Ottawa.). Only plants with spineless fruits have been reported from Turkmenistan (Uljanova, 1964Uljanova, T. 1964. A new species of medic from Turkmenia. Novosti Sistematiki Vysshikh Rastenii Moscow, Leningrad 175: 175–177.; Small & Brookes, 1991Small, E. & Brookes, B. 1991. A clarification of Medicago sinskiae Uljan. Canadian Journal of Botany 69: 100–106. https://doi.org/10.1139/b91-014). Mehregan et al. (2002Mehregan, I., Rahiminejad, M. R. & Azizian, D. 2002. A taxonomic revision of the genus Medicago L. (Fabaceae) in Iran. Iranian Journal of Botany 9: 207–221.) reported both spineless and spiny plants in Iran, with the spines up to 3.5 mm long and usually hooked at the apex.

Figure 1. Variation of pods of Medicago sinskiae as seen in selected samples from some populations. See Table 1 for abbreviations. [View full size] [Descargar tamaño completo]

Medicago sinskiae was first described from the western Kopet-Dagh (Turkmenistan) by Uljanova (1964Uljanova, T. 1964. A new species of medic from Turkmenia. Novosti Sistematiki Vysshikh Rastenii Moscow, Leningrad 175: 175–177.) based on very limited material from a locality in Turkmenistan (southwestern Turkmenia, Karal-Kalin, 10th western Kopet-Dag, Kuraty Canyon of the Sumbar-Chandyrskii watershed ridge, river slope, debris cone, 700 m, 1961). It was rarely accepted by botanists for nearly three decades until it was tentatively recognized by Small & Brookes (1991Small, E. & Brookes, B. 1991. A clarification of Medicago sinskiae Uljan. Canadian Journal of Botany 69: 100–106. https://doi.org/10.1139/b91-014) based on the specimens grown from the seeds of the type collection. They suggested that M. sinskiae was a distinctive species, derived from M. rigidula (L.) All–M. rigiduloides E. Small complex and was not related to other species of Medicago (Small & Brookes, 1991Small, E. & Brookes, B. 1991. A clarification of Medicago sinskiae Uljan. Canadian Journal of Botany 69: 100–106. https://doi.org/10.1139/b91-014). The chromosome number of both M. rigidula and M. sinskiae are n = 7 and 8 (Heyn, 1963Heyn, C. C. 1963. The annual species of Medicago. Magnes Press, The Hebrew University, Jerusalem.; Small & Brookes, 1991Small, E. & Brookes, B. 1991. A clarification of Medicago sinskiae Uljan. Canadian Journal of Botany 69: 100–106. https://doi.org/10.1139/b91-014). Like many other annual species of Medicago with small flowers, M. sinskiae seems to be largely an inbreeder (Novoselova, 2003Novoselova, L. V. 2003. Bud cleistogamy at annual species of genus Medicago. Czech Journal of Genetics and Plant Breeding 39: 285–287.; Small, 2011Small, E. 2011. Alfalfa and relatives. Evolution and classification of Medicago. NRC Research Press, Ottawa.).

For nearly four decades no further samples of M. sinskiae were collected again. Mehregan et al. (2002Mehregan, I., Rahiminejad, M. R. & Azizian, D. 2002. A taxonomic revision of the genus Medicago L. (Fabaceae) in Iran. Iranian Journal of Botany 9: 207–221.) reported some material from western Iran, hundreds of kilometers away from the type locality. Until collections made from the Zagrosian region of western Iran by Mehregan et al. (2002Mehregan, I., Rahiminejad, M. R. & Azizian, D. 2002. A taxonomic revision of the genus Medicago L. (Fabaceae) in Iran. Iranian Journal of Botany 9: 207–221.), the only known population of M. sinskiae was the type collection, and it was thought the species was only endemic to Kopet-Dagh. Based on morphological similarities, Mehregan et al. (2002Mehregan, I., Rahiminejad, M. R. & Azizian, D. 2002. A taxonomic revision of the genus Medicago L. (Fabaceae) in Iran. Iranian Journal of Botany 9: 207–221.) treated M. sinskiae, M. constricta Durieu, M. rigidula and M. rigiduloides as one species. Using ITS marker, Zareei et al. (2020Zareei, R., Small, E., Assadi, M. & Mehregan, I. 2020. The taxonomic status of Medicago sinskiae: insights from morphological and molecular data. Journal of Taxonomy and Biosystematics 12: 1–14. https://10.22108/tbj.2020.125443.1132) proposed that M. sinskiae is a separate species, sister to M. rigidula and M. rigiduloides. Furthermore, in many regional floras, both M. rigidula and M. rigiduloides are collectively treated as M. rigidula (Bayat et al., 2021Bayat, M., Assadi, M., Small, E. & Mehregan, I. 2021. Molecular studies of Iranian populations support the morphology-based taxonomic separation of Medicago rigidula and M. rigiduloides. Phytotaxa 518: 281–299. https://doi.org/10.11646/phytotaxa.518.4.5).

Different markers are available for studying populations of intra and inter species. SSR (simple sequence repeats) or microsatellites are widely used in studying the structure of plant populations and genetic diversity (Chabane et al., 2008Chabane, K., Varshney, R. K., Graner, A. & Valkoun, J. 2008. Generation and exploitation of EST-derived SSR markers for assaying molecular diversity in durum wheat populations. Genetic Resources and Crop Evolution 55: 869–881. https://doi.org/10.1007/s10722-007-9292-8; Enayat Avval, 2017Enayat Avval, S. 2017. Assessing polymorphism information content (PIC) using SSR markers on local species of Citrullus colocynthis. Case study: Iran, Sistan-Balouchestan Province. Journal of Molecular Biology Research 7: 42–49. https://doi.org/10.5539/jmbr.v7n1p42 ). The codominant SSRs are among the most reliable markers in population genetics studies (Freeland, 2020Freeland, J. R. 2020. Molecular ecology (3rd ed.). Wiley-Blackwell, London.). Next generation sequencing (NGS) technology is frequently used to identify microsatellite regions and develop SSR primers (Shendure & Ji, 2008Shendure, J. & Ji, H. 2008. Next-generation DNA sequencing. Nature Biotechnology 26: 1135–1145. https://doi.org/10.1038/nbt1486; Yang et al., 2015Yang, T., Fang, L., Zhang, X., Hu, J., Bao, S., Hao, J., Li, L., He, Y., Jiang, J., Wang, F., Tian, S. & Zong, X. 2015. High-throughput development of SSR markers from Pea (Pisum sativum L.) based on next generation sequencing of a purified Chinese commercial variety. PLoS ONE 10: e0139775. https://doi.org/10.1371/journal.pone.0139775; Emami-Tabatabaei et al., 2021Emami-Tabatabaei, S. S., Small, E., Assadi, M., Dehshiri, M. M. & Mehregan, I. 2021. Genetic variation among Iranian Medicago polymorpha L. populations based on SSR markers. Genetic Resources and Crop Evolution 68: 1411–1424. https://doi.org/10.1007/s10722-020-01071-7).

Our close examination of material recently collected from Iran showed that M. sinskiae has a wider distribution in Iran. This study aims to use SSR markers identified by NGS technologies to clarify the genetic structure of this species at population level.

MATERIALS AND METHODSTop

Sample collection

In total, pods of 74 individuals including 62 individuals from eight populations of M. sinskiae, five individuals from a single population of M. constricta and seven individuals from a single population of M. rigidula were collected from western and southwestern Iran in July 2017 (Table 1). Considering the limited occurrence of M. sinskiae, the restricted population size, and the criterion of individuals of each population being sampled at 20 m minimum intervals, no more effective number of individuals could be gathered. We could not find any material from Turkmenistan. All collected samples were identified and labeled based on authoritative identification keys (Heyn, 1984Heyn, C. C. 1984. Papilionaceae. In: Rechinger, K. H. (Eds.), Flora Iranica. Akademischer Druck-u. Verlagsanstalt, Graz: 253–271.; Small & Jomphe, 1989Small, E. & Jomphe, M. 1989. A synopsis of the genus Medicago (Leguminosae). Canadian Journal of Botany 67: 3260–3294. https://doi.org/10.1139/b89-405; Mehregan et al., 2002Mehregan, I., Rahiminejad, M. R. & Azizian, D. 2002. A taxonomic revision of the genus Medicago L. (Fabaceae) in Iran. Iranian Journal of Botany 9: 207–221.; Small, 2011Small, E. 2011. Alfalfa and relatives. Evolution and classification of Medicago. NRC Research Press, Ottawa.). Six to eight seeds of each individual were cultivated in separate pots on a research farm in southwestern Iran with similar ecological conditions to their natural habitats. Samples were taken at different growth stages up until fully ripened pods were developed. Total DNA was extracted from young leaves. Morphological examinations were performed on the fully grown plants and pods. Voucher specimens were deposited at IAUH (Islamic Azad University Herbarium).

DNA extraction

Total genomic DNA was extracted from young leaves dried in silica gel using CTAB (cetyltrimethylammonium bromide) method of Doyle & Doyle (1987Doyle, J. J. & Doyle, J. L. 1987. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochemical Bulletin 19: 11–15.) employing Nucleospin© Plants kits (Machery-Nagel, Germany) after manufacturer’s instructions. Success of DNA extraction was initially checked on 1% agarose gel. Density and purity of extracted DNA were examined on a NanoDrop™ 2000 (Thermo Scientific).

Identification of microsatellite loci via Next Generation Sequencing

Method of Yang et al. (2015Yang, T., Fang, L., Zhang, X., Hu, J., Bao, S., Hao, J., Li, L., He, Y., Jiang, J., Wang, F., Tian, S. & Zong, X. 2015. High-throughput development of SSR markers from Pea (Pisum sativum L.) based on next generation sequencing of a purified Chinese commercial variety. PLoS ONE 10: e0139775. https://doi.org/10.1371/journal.pone.0139775) was used to identify and develop SSR markers. 100 ng of the genomic DNA of a single sample was used to generate an Illumina DNA library. After DNA was fragmented, repaired at the ends, “A-tailed”, and ligated to the TruSeq adapters, the library was amplified in eight cycles. The average size of the library was 670 bp, corresponding to an average integral length of 500 bp. Sequencing of the “library” was carried out in an Illumina-MiSeq system (Illumina, San Diego, CA) with 300 bp each in the “paired-end” mode. In order to remove residues of adapter sequences, the overlapping “paired-end-reads” were first trimmed at the ends. The quality score was set to at least 20. FLASH software (Magoč & Salzberg, 2011Magoč, T. & Salzberg, S. L. 2011. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27: 2957–2963. https://doi.org/10.1093/bioinformatics/btr507) was used to assemble the reads. The resulting sequences were bioinformatically analyzed for existing microsatellites (Faircloth, 2008Faircloth, B. C. 2008. MSATCOMMANDER: Detection of microsatellite repeat arrays and automated, locus-specific primer design. Molecular Ecology Resources 8: 92–94. https://doi.org/10.1111/j.1471-8286.2007.01884.x). Examination of 817,510 potential loci with any tandem repeat were performed with following criteria: mismatch = 0, motive-length = 3–5 and length of repeats = 40–185.

Primer design and test of polymorphism

Primer3 software (Koressaar et al., 2018Koressaar, T., Lepamets, M., Kaplinski, L., Raime, K., Andreson, R. & Remm, M. 2018. Primer3_masker: integrating masking of template sequence with primer design software. Bioinformatics 34: 1937–1938. https://doi.org/10.1093/bioinformatics/bty036) was used to design primers for 513 out of 2410 eligible loci suitable for primer designing with the following criteria: PCR-product size = 150–350 bp, primer length = 18–22 bp, and TM-value = 58–62°C. Primer pairs were synthesized for 124 loci and tested on 12 samples for suitability. Those primer pairs which yielded unique PCR fragments were tested on a batch of four individuals. In the best case those primers should result in different unique fragments in each individual. PCR products from individuals were sequenced to confirm the tandem repeat pattern. Out of the 124 primer pairs only six polymorphic loci could be detected. Primers with different labeling were synthesized (Table 2). In order to avoid the possible problems with null alleles, we tried to improve our methods by using labelled primers, finding the best annealing temperature, and sequencing the PCR products for estimating the product size (Dakin & Avise, 2004Dakin, E. & Avise, J. 2004. Microsatellite null alleles in parentage analysis. Heredity 93: 504–509. https://doi.org/10.1038/sj.hdy.6800545).

PCR and fragment analysis

Polymerase chain reactions (PCR) were carried out in 25 μl reaction volumes containing 5 μl of 5× PCR-buffer, 2 μl of MgCl₂ (25mM), 2 μl of dNTPs (2.5 mM), 0.1 μl of Taq-Polymerase (1.25 U), 1 μl of forward primer (10 pmol/μl), 1 μl of reverse primer (10 pmol/μl), 1 μl of genomic DNA (~15 ng), and 12.9 μl of ddH₂O. The PCR reactions were performed on a Labcycler Gradient (SensoQuest GmbH, Göttingen, Germany) under the following conditions: initial denaturation at 95°C for 180 s, 34 cycles of denaturation at 95°C for 30 s, annealing at primer specific temperature for 30 s and extension at 72°C for 30 s, and a final extension at 72°C for 300 s. PCR products of MED-01/MED-02/MED-03 and MED-04/MED-05/MED-06 primer pairs were pooled separately. Two μl of each pool was mixed with 7.75 μl of HiDi formamide (Applied Biosystems) and 0.25 μl ROX-500 internal size standard (Applied Biosystems) and then injected to an 3730xl Applied Biosystems capillary sequencer. Raw data were visualized with GeneMarker v4.0 (Applied Biosystems, Foster City, CA, USA). Output files were aligned with the ROX-500 size standard using GeneMarker v2.4.2 (GeneMarker, SoftGenetics, State College, PA, USA). Each peak with a signal intensity of more than 1000 was scored as present. The binary matrices (1: presence, 0: absence) of each two primer pools were combined and prepared for further analyses. Some analyses need data to be entered as co-dominant. To do so, allele sizes for each locus were entered in GenAlEx software after publisher’s tutorials.

Multivariate analyses

UPGMA (Unweighted Pair Group Method with Arithmetic mean) algorithm of clustering with Dice similarity index as well as PCoA (Principal Coordinate Analysis) analyses of dataset were performed with PAST3 software package (Hammer et al., 2001Hammer, Ø., Harper, D. A. T. & Ryan, P. D. 2001. PAST: Paleontological statistics software package for education and data analysis. Palaeontologia Electronica 4: 4. Retrieved May 20, 2021, from http://palaeo-electronica.org/2001_1/past/issue1_01.htm). After circumscribing M. sinskiae with UPGMA and PCoA, further analysis was performed on M. sinskiae populations only. In POPTREEW software genetic distances were measured and phylogenetic tree was constructed.

Analysis of population structure

Genetic structure of M. sinkiae populations was estimated using a Bayesian Markov Chain Monte Carlo model (MCMC) implemented in Structure v2.3.4 (Pritchard et al., 2000Pritchard, J. K., Stephens, M. & Donnelly, P. 2000. Inference of population structure using multilocus genotype data. Genetics 155: 945–959. https://doi.org/10.1093/genetics/155.2.945). The true number of subpopulations (K) was calculated using Evanno et al. (2005Evanno, G., Regnaut, S. & Goudet, J. 2005. Detecting the number of clusters of individuals using the software structure: a simulation study. Molecular Ecology 14: 2611–2620. https://doi.org/101111/j.1365-294X.2005.02553.x ) method summarized in CLUMPP_Windows v1.1.2 software (Jakobsson & Rosenberg, 2007Jakobsson, M. & Rosenberg, N. A. 2007. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23: 1801–1806. https://doi.org/10.1093/bioinformatics/btm233) on the Structure Harvester website (Earl & vonHoldt, 2012Earl, D. A. & vonHoldt, B. M. 2012. Structure Harvester: a website and program for visualizing Structure output and implementing the Evanno method. Conservation Genetics Resources 4: 359–361. https://doi.org/10.1007/s12686-011-9548-7). It was tested for K = 1 to K = 8 with 20 independent simulations at 60,000 samplings with a burn-in period of 10,000 first iterations. The final analysis with the resulting K (= 3) was conducted by 1,250,000 repetitions after a burn-in of first 500,000 replications. Individuals with at least 80% probability of membership in a cluster were considered to belong to that cluster. Individuals with probabilities of membership below 80% were interpreted as a hybrid genotype.

Estimating frequencies, diversity and population structure

Different parameters of M. sinskiae populations including the haploid number of migrants (Nm), number of different alleles (N_a), number of effective alleles (N_e), number of private alleles (N_p), expected heterozygosity (H_e), unbiased expected heterozygosity (uH_e), and Shannon information Index (I) for each population were calculated using GenAlEx v6.503 (Peakall & Smouse, 2012Peakall, R. & Smouse, P. E. 2012. GenAlEx 6.5: genetic analysis in excel. Population genetic software for teaching and research an-update. Bioinformatics 28: 2537–2539. https://doi.org/10.1093/bioinformatics/bts460) and POPGENE (Yeh et al., 1999Yeh, F. C., Yang, R. C., Boyle, T. B. J., Ye, Z. H. & Mao, J. X. 1999. POPGENE, version 1.32: the user friendly shareware for population genetic analysis. Biotechnology and Molecular Biology 10: 295–301.). F-statistics employ inbreeding coefficients to describe the partitioning of genetic variation within and among populations and can be calculated at three different levels (F_IS, F_ST, F_IT). GenAlEx was also used to calculate the mean of haploid number of migrants or gene flow (Nm), genetic differentiation between subpopulations (F_ST), inbreeding coefficient of an individual relative to the subpopulation (F_IS), and inbreeding coefficient of an individual relative to the total population (F_IT and G_ST) for each primer pair (locus). G_ST is assumed to be an analogue of F_ST, and G_ST is equivalent to F_ST when there are only two alleles per locus and is the weighted average of F_ST for all alleles in case of multiple alleles per locus (Freeland, 2020Freeland, J. R. 2020. Molecular ecology (3rd ed.). Wiley-Blackwell, London.). F_IS was estimated by dividing (H_e – H_o)/H_e (Pagnotta, 2018Pagnotta, M. A. 2018. Comparison among methods and statistical software package to analyze germplasm genetic diversity by means of codominant markers. J 1: 197–215. https://doi.org/10.3390/j1010018). The selfing rates (S) were calculated by dividing 2F_IS/ (1+ F_IS) in our populations (Burkil et al., 2017Burkil, A., Sieber, N., Seppala, K. & Jokela, J. 2017. Comparing direct and indirect selfing rate estimates: when are population-structure estimates reliable. Heredity 118: 525–533. https://doi.org/10.1038/hdy.2017.1).

The polymorphism information content (PIC) was also calculated manually in spreadsheet as PIC = 1 – Σⁿ_i = 1p_i², where p_i (i = 1, 2, 3, … I) are the frequencies of i^th alleles for the given locus and “i” is the number of distinct alleles at a locus (Luo et al., 2019Luo, Z., Brock, J., Dyer, J. M., Kutchan, T., Schachtman, D., Augustin, M., Ge, Y., Fahlgren, N. & Abdel-Haleem, H. 2019. Genetic diversity and population structure of a Camelina sativa spring panel. Frontiers in Plant Science 10: 184. https://doi.org/10.3389/fpls.2019.00184). Genetic distances and genetic diversity with pairwise test and Nei’s Genetic Distance (GD) were calculated using POPGENE and GenAlEx. Data were analyzed with POPTREE2 (Takezaki et al., 2009Takezaki, N., Nei, M. & Tamura, K. 2009. POPTREE2: Software for constructing population trees from allele frequency data and computing other population statistics with Windows interface. Molecular Biology and Evolution 27: 747–752. https://doi.org/10.1093/molbev/msp312). A neighbor-joining (NJ) method based on Nei’s genetic distances using POPTREE2 with 1000 replicates of bootstrapping was used to illustrate the relationships between populations. The SplitsTree v4.15.1 software (Huson & Bryant, 2006Huson, D. H. & Bryant, D. 2006. Application of phylogenetic networks in evolutionary studies. Molecular Biology and Evolution 23: 254–267. https://doi.org/10.1093/molbev/msj030) was used to calculate the genetic distances for a split neighbor net. The genetic difference within and among the studied populations (Table 1) as well as regions (Region 1, western Iran: populations ABD, KHW, KHR, SPC, SPD, and PLS; Region 2, southwestern Iran: populations FTH, BSN) was tested by AMOVA (Analysis of Molecular Variance) with 1000 permutations using GenAlEx. In order to evaluate the impact of genetic distance and geographical distance on population differentiation, a Mantel test (Mantel, 1967Mantel, N. 1967. The detection of disease clustering and a generalized regression approach. Cancer Research 27: 209–220. ) was performed to correlate two matrixes of genetic distance based on F_ST and Nei’s genetic distance using GenAlEx.

RESULTSTop

Multivariate analyses

The UPGMA dendrogram of SSR data with Dice similarity index showed that the 74 individuals sampled fell into three distinct groups (Fig. 2). Each of these groups represented the individuals of a different species. Similar results were also observed in the PCoA analysis of the same dataset (Fig. 3). In both analyses, all individuals of M. sinskiae were circumscribed as a single distinct group and therefore is in accordance with the proposal of Small & Brookes (1991Small, E. & Brookes, B. 1991. A clarification of Medicago sinskiae Uljan. Canadian Journal of Botany 69: 100–106. https://doi.org/10.1139/b91-014) and the results of Zareei et al. (2020Zareei, R., Small, E., Assadi, M. & Mehregan, I. 2020. The taxonomic status of Medicago sinskiae: insights from morphological and molecular data. Journal of Taxonomy and Biosystematics 12: 1–14. https://10.22108/tbj.2020.125443.1132). Further analyses were performed on the 62 individuals of M. sinskiae (see Figs. 4–10 and Tables 4–6).

Figure 2. UPGMA tree based on the SSR analysis of 74 individuals of Medicago sinskiae, M. constricta and M. rigidula. [View full size] [Descargar tamaño completo]

Figure 3. Principal Coordinate Analysis (PCoA) plot of 74 individuals of Medicago sinskiae, M. constricta and M. rigidula. [View full size] [Descargar tamaño completo]

Genetic diversity

A total of 19 alleles were detected from the analysis of six loci (Table 3). The number of alleles per locus ranged from 2 to 6. The mean of Shannon Index (I) for eight populations was 0.458 (Table 4). The Shannon Index (I) values were higher than expected heterozygosity (H_e) values. Population ABD showed the maximum values of N_a (2.667), N_e (1.870), I (0.701), H_e (0.418) and uH_e (0.440), and population KHR showed the minimum values of N_a (1.333), N_e (1.308), I (0.224), H_o (0), H_e (0.160) and uH_e (0.178). The mean of observed heterozygosity (H_o), expected heterozygosity (H_e), unbiased expected heterozygosity (uH_e), number of different alleles (N_a), number of effective alleles (N_e) and fixation Index (F) for all eight populations were 0.040, 0.296, 0.317, 1.917, 1.535, and 0.849 respectively (Table 4). The F_IS value ranged from 0.582 to 1 in eight populations and population BSN showed the lowest level of F_IS (0.582), compared to the highest levels of F_IS (1) observed in populations PLS and KHR. Selfing rate (S) ranged from 0.735 to 1 with a mean value of 0.928 (Table 4). Nei’s genetic distance (GD) results (Table 5) showed the populations with more similar alleles having smaller genetic distances. Populations KHW and SPC both from Lurestan province were the most similar (Nei’s GD = 0.013). The largest amount of genetic distance was observed between populations BSN and ABD (Nei’s GD = 0.223).

Mean of polymorphism percentage was 77.08% (Table 4). Although populations ABD, KHW, SPC and FTH showed to be 100% polymorphic, population KHR showed the lowest polymorphism (33.3%) as well as lowest values of genetic parameters (Table 4). Our SSR markers displayed a high level of polymorphism, and this species seems to be a polymorphic plant. Both PIC and H_e (= gene diversity) values are measures of genetic diversity, although PIC values are not useful in linkage analyses when determining the inheritance between offspring and parental genotypes, and expected heterozygosity (H_e) is useful for haploid markers (Luo et al., 2019Luo, Z., Brock, J., Dyer, J. M., Kutchan, T., Schachtman, D., Augustin, M., Ge, Y., Fahlgren, N. & Abdel-Haleem, H. 2019. Genetic diversity and population structure of a Camelina sativa spring panel. Frontiers in Plant Science 10: 184. https://doi.org/10.3389/fpls.2019.00184). The mean of following parameters was observed for six loci: haploid number of migrants or gene flow (Nm) = 1.999, coefficient of genetic differentiation among populations (F_ST) = 0.157, inbreeding coefficient of an individual relative to the subpopulation (F_IS) = 0.820, inbreeding coefficient of an individual relative to the total population (F_IT) = 0.847, and among-population genetic differentiation (G_ST) = 0.402. Locus MED-02 showed the highest value of F_IS, G_ST and Nm and the highest value of F_ST and F_IT were observed in locus MED-04 (Table 3). The SSR markers showed the PIC values ranging from 0.196 (locus MED-01) to 0.999 (locus MED-05) (Table 3). Among the primers, locus MED-06 with six alleles had the highest number of polymorphic bands and locus MED-01, MED-02, MED-03 with two alleles had the lowest number of polymorphic bands (Table 3). The PIC value was used to measure the informativeness of primers. Having the highest PIC value (0.999), locus MED-05 showed higher polymorphism and had more impact in differentiation of individuals. The minimum amount of PIC value (0.196) was observed in the monomorphic locus MED-01, which was uniform in all individuals. All other primers were polymorphic (Table 3).

Population structure and genetic relationships

The true number of subpopulations (K = 3, Fig. 4) was obtained using the method of Evanno et al. (2005Evanno, G., Regnaut, S. & Goudet, J. 2005. Detecting the number of clusters of individuals using the software structure: a simulation study. Molecular Ecology 14: 2611–2620. https://doi.org/101111/j.1365-294X.2005.02553.x ). As seen in Fig. 4, none of the populations seems to be uniformly consisting of a single cluster. All populations had individuals from different genetic clusters. Results of UPGMA clustering analysis showed that there was no major cluster formed by individuals solely from a single population or individuals with geographical proximity (Fig. 5). Different clusters included individuals from different populations with no geographical proximity. For example, individuals of population SPC were present in all four major clusters A1a, A1b, A2 and B. As seen in Figs. 6 and 7, individuals of different populations were scattered all over the PCoA plot and Neighbor-net network. In accordance with the results of Structure analysis (Fig. 4), UPGMA clustering dendrogram (Fig. 5), the PCoA plot (Fig. 6) and Neighbor-net network (Fig. 7), none of populations studied were unmixed. Despite being clearly circumscribed when analyzed alongside with M. rigidula and M. constricta, populations of M. sinskiae showed no grouping based on geographical proximity. In POPTREE software genetic distances measured for constructing phylogenetic trees of eight populations showed no relationship based on geographical proximity (Fig. 8). The AMOVA test was performed (Table 6) to study population differentiation and to estimate the percentage of intrapopulation and interpopulation genetic variation. Most of the genetic variation occurred among individuals (95%). The calculated genetic variation among populations was 5%, and the number for variation among regions was 0% (Table 6). The Mantel test indicated no meaningful correlation between genetic distance and geographical distribution (R = -0.019, P = 0.659).

Figure 4. Geographical distribution and population structure of Medicago sinskiae populations in Iran based on K = 3. Magnitude of Delta K as a function of K = 2–8 is shown at the upper corner on left. [View full size] [Descargar tamaño completo]

Figure 5. Unweighted Pair Group Method with Arithmetic mean (UPGMA) tree based on the SSR analysis of 62 individuals of Medicago sinskiae in Iran. [View full size] [Descargar tamaño completo]

Figure 6. Principal Coordinate Analysis (PCoA) of Iranian populations of Medicago sinskiae. [View full size] [Descargar tamaño completo]

Figure 7. Neighbor-net network of individuals of Medicago sinskiae generated from the complement of Dice similarity coefficient. [View full size] [Descargar tamaño completo]

Figure 8. Dendrogram based on Nei’s Genetic Distance among the studied populations of Medicago sinskiae. [View full size] [Descargar tamaño completo]

DISCUSSIONTop

We showed that SSR markers have the potential to identify and separate the closely related species M. constricta, M. rigidula and M. sinskiae. SSR markers have been widely used to evaluate genetic diversity and polymorphism in plants. In our study of six loci in Medicago sinskiae, the mean of polymorphism in genotypes was remarkably high (77.08%) and was 100% in some populations. This indicates that M. sinskiae is a polymorphic species. Higher genetic diversity is usual in outcrossing species (Szczecińska et al., 2016Szczecińska, M., Sramko, G., Wołosz, K. & Sawicki, J. 2016. Genetic diversity and population structure of rare and endangered plant species Pulsatilla patens (L.) Mill. in East Central Europe. PLoS ONE 11: 1–24. https://doi.org/10.1371/journal.pone.0151730). In a population genetic study of outcrossing Marrubium L., Salehi et al. (2018Salehi, N., Kharazian, N. & Shiran, B. 2018. Genetic diversity of Marrubium species from Zagros region (Iran), using inter simple sequence repeat molecular marker. Journal of Sciences, Islamic Republic of Iran 29: 7–19.) reported a high level of polymorphism (100% vs. our 77.08%), genetic diversity (G_ST = 0.99 vs. our 0.4) and Shannon information index (I = 0.51 vs our 0.45). Our obtained values are quite high for species like annual Medicago that are are known to be inbreeders (Small, 2011Small, E. 2011. Alfalfa and relatives. Evolution and classification of Medicago. NRC Research Press, Ottawa.). A lower heterozygosity (H_e = 0.348–0.479) was observed for self-pollinated annual species M. truncatula Gaertn. in the French Mediterranean region. Although in self-pollinated species pollen dispersal is very infrequent compared to outcrossing species, it is important because it generates novel genetic variability via recombinant lines (Bonnin et al., 2001Bonnin, I., Ronfort, J., Wozniak, F. & Olivieri, I. 2001. Spatial effects and rare outcrossing events in Medicago truncatula (Fabaceae). Molecular Ecology 10: 1371–1383. https://doi.org/10.1046/j.1365-294x.2001.01278.x). Riday et al. (2015Riday, H., Reisen, P., Raasch, J. A., Santa-Martinez, E. & Brunet, J. 2015. Selfing rate in Alfalfa seed production field pollinated with leafcutter bees. Crop Science 55: 1087–1095. https://doi.org/10.2135/cropsci2014.04.0295) reported a variable selfing rate (0 to 52.2%; mean 11.8%) in populations of perennial Medicago sativa L. Besides detecting a very high (99%) rate of selfing and a small seed dispersal distance in Medicago truncatula, Siol et al. (2008Siol, M., Prosperi, J. M., Bonnin, I. & Ronfort, J. 2008. How multilocus genotypic pattern helps to understand the history of selfing populations: a case study in Medicago truncatula. Heredity 100: 517–525. https://doi.org/10.1038/hdy.2008.5) reported high genotypic diversity and polymorphism in some few inbred lines per population. The mean value of F_IS in M. sinskiae is 0.923 (ranging from 0.582 to 1.000). Populations PLS and KHR showed higher levels of selfing rates (S = 1.000). Higher levels of gene flow (Nm) and population structure seems to have been more affected by dispersal patterns and not localized gene flow (Bayat et al., 2021Bayat, M., Assadi, M., Small, E. & Mehregan, I. 2021. Molecular studies of Iranian populations support the morphology-based taxonomic separation of Medicago rigidula and M. rigiduloides. Phytotaxa 518: 281–299. https://doi.org/10.11646/phytotaxa.518.4.5; Emami-Tabatabaei et al., 2021Emami-Tabatabaei, S. S., Small, E., Assadi, M., Dehshiri, M. M. & Mehregan, I. 2021. Genetic variation among Iranian Medicago polymorpha L. populations based on SSR markers. Genetic Resources and Crop Evolution 68: 1411–1424. https://doi.org/10.1007/s10722-020-01071-7; Bagheri et al., 2022Bagheri, Z., Assadi, M., Small, E. & Mehregan, I. 2022. Cryptic molecular-geographical divergence within Medicago minima revealed by SSR markers. Iranian Journal of Science and Technology, Transactions A: Science 46: 49–60. https://doi.org/10.1007/s40995-021-01236-8). Small increases in gene flow (Nm) will reduce population differentiation (F_ST). Each individual in a selfing population differs from others. The accessions preserve their traits, and it will continue in the coming generation. It was shown in M. truncatula, when selfing rates are very large, the genetic and genotypic diversity can be high, while selfers are composed of a few inbred lines per population (Bataillon & Ronfort, 2006Bataillon, T. & Ronfort, J. 2006. Evolutionary and ecological genetics of Medicago truncatula. In: Mathesius, U., Journet, E. P. & Sumner, L. W. (Eds.), The Medicago truncatula handbook. Noble Research Institute, Ardmore: 1–13.; Siol et al., 2008Siol, M., Prosperi, J. M., Bonnin, I. & Ronfort, J. 2008. How multilocus genotypic pattern helps to understand the history of selfing populations: a case study in Medicago truncatula. Heredity 100: 517–525. https://doi.org/10.1038/hdy.2008.5). Yan et al. (2009Yan, J., Chu, H.-J., Wang, H.-C., Li, J.-Q. & Sang, T. 2009. Population genetic structure of two Medicago species shaped by distinct life form, mating system and seed dispersal. Annals of Botany 103: 825–834. https://doi.org/10.1093/aob/mcp006) showed that self-pollination and dispersal mechanisms shaped the population genetic structure and geographical distribution of Medicago lupulina L. (F_ST = 0.535) and Medicago ruthenica (L.) Trautv. (F_ST = 0.130), and F_ST in self-pollinated annual species is higher than outcrossing perennials. Selfing rates estimated from F_IS values were more than 95% for M. lupulina but much lower (ca. 30%) for M. ruthenica (calculated from Yan et al., 2009Yan, J., Chu, H.-J., Wang, H.-C., Li, J.-Q. & Sang, T. 2009. Population genetic structure of two Medicago species shaped by distinct life form, mating system and seed dispersal. Annals of Botany 103: 825–834. https://doi.org/10.1093/aob/mcp006). It is shown that selfing species of Zingiber Mill. have less genetic diversity at the population and species levels compared to outcrossing ones (Huang et al., 2019Huang, R., Chu, Q. H., Lu, G. H. & Wang, Y. Q. 2019. Comparative studies on population genetic structure of two closely related selfing and outcrossing Zingiber species in Hainan Island. Scientific Reports 9: 17997. https://doi.org/10.1038/s41598-019-54526-y).

Medicago’s close relative genera Trigonella L., Melilotus Mill. and Trifolium L. have a passive floral pollination mechanism allowing flowers to be pollinated frequently. In contrast, the explosive tripping mechanism of pollination in the genus Medicago allows the flowers to be visited by pollinators only once (Small, 2011Small, E. 2011. Alfalfa and relatives. Evolution and classification of Medicago. NRC Research Press, Ottawa.). The genus Medicago includes both perennial and annual species. The perennial Medicago sativa is chiefly an outcrossing species with some populations benefitting from self-pollination. The floral structure of most of the annual species of Medicago is related to their mostly self-pollination nature. The annuals have flowers that may be closed (cleistogamous), although they are usually opened (chasmogamous) and auto tripping (Novoselova, 2003Novoselova, L. V. 2003. Bud cleistogamy at annual species of genus Medicago. Czech Journal of Genetics and Plant Breeding 39: 285–287.). Outcrossing, at least partly, is present in nearly all perennial species of Medicago. In contrast, all the annual species of the genus Medicago seem to be strongly self-pollinated, with limited association with pollinators (Small, 2011Small, E. 2011. Alfalfa and relatives. Evolution and classification of Medicago. NRC Research Press, Ottawa.). Self-pollination, occurring in different ways, can have lasting impacts on genetic diversity. The morphological and phenological characteristics of flowers have impact on each mode of self-pollination (Lloyd & Schoen, 1992Lloyd, D. G. & Schoen, D. J. 1992. Self-fertilization and cross-fertilization in plants. I. Functional dimensions. International Journal of Plant Sciences 153: 358–369. https://doi.org/10.1086/297040). In the absence of disturbance, migration events can partition populations into several independent recombinant lines, allowing a high level of genetic diversity to be sustained (Bonnin et al., 2001Bonnin, I., Ronfort, J., Wozniak, F. & Olivieri, I. 2001. Spatial effects and rare outcrossing events in Medicago truncatula (Fabaceae). Molecular Ecology 10: 1371–1383. https://doi.org/10.1046/j.1365-294x.2001.01278.x). The special genetic structure of M. sinskiae populations unrelated to geographical proximity is unlike many other annual self-pollinated medics (Bayat et al., 2021Bayat, M., Assadi, M., Small, E. & Mehregan, I. 2021. Molecular studies of Iranian populations support the morphology-based taxonomic separation of Medicago rigidula and M. rigiduloides. Phytotaxa 518: 281–299. https://doi.org/10.11646/phytotaxa.518.4.5; Emami-Tabatabaei et al., 2021Emami-Tabatabaei, S. S., Small, E., Assadi, M., Dehshiri, M. M. & Mehregan, I. 2021. Genetic variation among Iranian Medicago polymorpha L. populations based on SSR markers. Genetic Resources and Crop Evolution 68: 1411–1424. https://doi.org/10.1007/s10722-020-01071-7; Bagheri et al., 2022Bagheri, Z., Assadi, M., Small, E. & Mehregan, I. 2022. Cryptic molecular-geographical divergence within Medicago minima revealed by SSR markers. Iranian Journal of Science and Technology, Transactions A: Science 46: 49–60. https://doi.org/10.1007/s40995-021-01236-8) and should be explained differently.

The presence of annual Medicago in Iran is well documented (Heyn, 1963Heyn, C. C. 1963. The annual species of Medicago. Magnes Press, The Hebrew University, Jerusalem.; Mehregan et al., 2002Mehregan, I., Rahiminejad, M. R. & Azizian, D. 2002. A taxonomic revision of the genus Medicago L. (Fabaceae) in Iran. Iranian Journal of Botany 9: 207–221.). There is no record of M. sinskiae in Iran in the literature published before 2002 and searching for M. sinskiae in major Iranian herbaria were unsuccessful. First presence of M. sinskiae in Iran was spotted in 1999 by collecting some pods from western regions represented in this study by populations ABD, KHW, KHR, SPC, SPD and PLS. Localities BSN and FTH are among the regions searched for medics by the authors between 1994 and 1997, where no material matching the description of M. sinskiae was collected. This study is based on the new material collected from western and southwestern Iran in 2016–2017. Once limited to western Iran, M. sinskiae started to appear in southwestern Iran. The historical relationships of the Iranian and Turkmen populations are unclear. They could be the remainder of a species that was once widespread, or a relatively newly generated species that is now expanding its range; in any case, M. sinskiae is clearly expanding in Iran. Inbreeding species such as annual medics often have enhanced ability to rapidly expand because they do not require pollinators (Barrett et al., 2008Barrett, S. C. H., Colautti, R. I. & Eckert, C. G. 2008. Plant reproductive system and evolution during biological invasion. Molecular Ecology 17: 373–383. https://doi.org/10.1111/j.1365-294x.2007.03503.x; Kalisz et al., 2004Kalisz, S., Vogler, D. W. & Hanley, K. M. 2004. Context-dependent autonomous self-fertilization yield reproductive assurance and mixed mating. Nature 430: 884–887. https://doi.org/10.1038/nature02776). Lower genetic diversity with self-fertilization combined with human-mediated dispersal boosted rapid expansion of Brassica tournefortii Gouan in the United States (Winkler et al., 2019Winkler, D. E., Chapin, K. J., François, O., Garmon, J. D., Gaut, B. S. & Huxman, T. E. 2019. Multiple introductions and population structure during the rapid expansion of invasive Sahara mustard (Brassica tournefortii). Ecology and Evolution 9: 7928–7941. https://doi.org/10.1002/ece3.5239). When a species expands rapidly, distinct clusters based on geographical proximity cannot be distinguished. We suggest that similar structures observed in populations of M. sinskiae (Fig. 5) would not be the consequence of gene flow, migration, and connectivity of populations; rather, they were originated by the rapid expansion. This hypothesis agrees with AMOVA results (Table 6), which show very low genetic diversity among populations (5%) and regions (0%). Annual medics such as M. sinskiae are excellent plants for feeding livestock (Khassanov, 1972Khassanov, O. Kh. 1972. Wild alfalfas of central Asia. Academy of Sciences of the Uzbek SSR, Tashkent.). They can disperse across relatively large distances via animal fur because of their spines and can dispersed with wind, rivers, and human activities. Medicago sinskiae is rapidly expanding in the western part of Zagrosian regions of Iran, a region dominated by oak forests (Zohary, 1973Zohary, M. 1973. Geobotanical foundations of the Middle East. Gustav Fischer Verlag, Stuttgart.). This area shows a rich diversity of wild and domestic animal life including sheep and goats. Life of many rural and nomad people of this area are dependent on grazing. The indehiscent spiny fruits of the annual medics are well adapted to dispersal in animal fur (Small, 2011Small, E. 2011. Alfalfa and relatives. Evolution and classification of Medicago. NRC Research Press, Ottawa.). This would explain why M. sinskiae in western Iran is expanding so fast. Given that climate change is now rapidly changing the distribution range of many species (Kelly & Goulden, 2008Kelly, A. E. & Goulden, M. L. 2008. Rapid shifts in plant distribution with recent climate change. Proceedings of the National Academy of Sciences of the United States of America 105: 11823–11826. https://doi.org/10.1073/pnas.0802891105; Gómez-Ruiz & Lacher Jr., 2019Gómez-Ruiz, E. P. & Lacher Jr., T. E. 2019. Climate change, range shifts, and the disruption of a pollinator-plant complex. Scientific Reports 9: 14048. https://doi.org/10.1038/s41598-019-50059-6), it will be interesting to follow the future distribution pattern of the species.

CONCLUSIONSTop

The genetic structure of M. sinskae is consistent with inbreeding, at least in Iranian populations, and it seems to be a species expanding its range. Medicago sinskiae may continue to expand into different regions of western and southwestern Iran, northern Iraq, and even southeastern Turkey, where its potential presence should be monitored. Furthermore, and despite further sampling is needed, results of our clustering and PCoA analyses suggest that M. sinskiae can be recognized as a separate species as it is differentiated from Iranian populations of M. rigidula and M. constricta.

ACKNOWLEDGMENTSTop

The authors would like to thank Dr. Sven Bikar, Dr. Bettina Ebner, and Dr. Tilmann Laufs for their helps in part of the bioinformatic analyses. We also would like to thank Dr. Javier López-Alvarado and two anonym reviewers for their valuable comments on an earlier draft of the manuscript.

REFERENCESTop