Sampling

In order to locate P. culiacanum populations, different data bases were consulted: GIFB (Global Information of Biodiversity, 2013), Missouri Botanical Garden, University of Arizona, the Biology School and the Faculty of Agriculture of the Autonomous University of Sinaloa, and the Sonora Regional Station of the National Autonomous University of Mexico (ERNO-UNAM). Areas that comply with the distribution characteristics were also explored. In each population, 10 to 30 individuals were collected (with at least 10 m between them), they were georeferenced, and stored at -80 °C, until DNA purification was performed.

Amplification

For the amplification, the extraction methods proposed by ^{Doyle (1991)}, ^{Porebski et al. (1997)}, ^{Sá et al. (2011)}, ^{Cota-Sánchez et al. (2006)}, ^{Krizman (2006}), ^{Sánchez-Hernández (2006}), and ^{Amani et al. (2011}) were evaluated. The appropriate methodology was based on that described by ^{Azmat et al. (2012}).

Due to the lack of previous information on P. culiacanum’s genetic variability, several chloroplast DNA regions reported in other parasitic plants were evaluated, where sufficient variability to perform population genetics studies were found (^{Taberlet et al., 1991}, ^{de Pamphilis et al., 1997}; ^{Zuber et al., 2000}; ^{Heinze, 2007}; ^{Amico et al., 2009}; ^{Nickrent et al., 2009}). Based on those results, c-Tab F and f-Tab R from ^{Taberlet et al. (1991)} were selected as strand primers, which amplify the intergenic region of Trn L-F. A touchdown-like PCR amplification protocol was used (^{Amico et al., 2007}), with a modification in alignment temperatures: from 48 to 50 °C, and from 52 to 54 °C. Amplification reactions were performed in a Techne Unit Genius^® thermocycler, visualized in a 1 % agarose gel, and stained with ethidium bromide (0.003 %) in an UV transilluminator.

The amplified product was sequenced in the HighThroughput Sequencing Solutions of the Genomic Science Department of the University of Washington, Seattle, U.S. To edit and align the sequences, two software were used: CodonCode Aligner (^{CodonCode Corporation, 2009}), Finch TV 1.4.0 (Geospiza Inc., 2004-2006), and ClustalW multiple alignment from the BioEdit software suite (^{Hall, 1999}).

Genetic variability analysis

DnaSP (Librado and Rosas, 2010) and ARLEQUIN 3.1 (^{Excoffier, 2005}) software were used to calculate genetic diversity indexes: 1) nucleotide diversity (π) (Nei, 1987, ^{Hedrick, 2005}, ^{Hartl and Clark, 2007}); 2) haplotype diversity (Hd), with and without indels as informative sites (Nei, 1987); 3) variation components due to differences between and within populations (AMOVA). In order to evaluate if the populations fit in the isolation by distance model, correlation and linear regression analysis of the genetic and geographic distances were performed, using the Maximum Composite Likelihood model, without taking indels into consideration, in MEGA 4.0, ^{Tamura et al. (2007}), and Microsoft 2007.

Haplotype network

Haplotype sequences were obtained with the DnaSP software, excluding the aligned sites with gaps, missing data, and monomorphic sites. The Network 4.6.1.1 (Fluxus Technology, 2004-2012) software was used to build haplotype networks (^{Bandelt et al., 1999}; ^{Polzin et al., 2003}) with Lennoa sequences as base group.

Haplotype phylogeny

The phylogenetic tree was generated with MEGA 4.0 (^{Tamura et al., 2007}), using the Neighbor-Joining method, with a bootstrap of one thousand iterations. Genetic distances were calculated with the Maximum Composite Likelihood method, without taking into consideration indels. Lennoa sp. sequences were used as an outgroup.

Sampling

Our P. culiacanum distribution area scouting found seven populations: five in Sinaloa and two in Sonora. Out of all the collected individuals, 10 of each population were selected for DNA purification. In some reported distribution sites (GIF, 2013), no populations were found, due to changes in land use - such as primary forests transformed into agricultural or forest exploitation areas (Ranchería Los Álamos, Cerro Siete Gotas, in Culiacán, San Blas)-, to the construction of highways -highway in south Culiacán, Los Mochis-Obregón highway, AltataDautillos highway-, or to tourism and real estate development projects (Nuevo Altata costal area), and the consequent ecosystem fragmentation.

Amplification

Out of all the sequenced samples, 70 belonged to P.culiacanum and two belonged to Lennoa madreporoides. The size of the amplified fragment was 735 nucleotides, 261 had gaps or were missing data, 151 were monomorphic sites and 323 were polymorphic sites. Two-hundred fifty-one of those polymorphic sites presented two types of substitutions, 69 presented three types, and 3 presented four types. The length of the amplification product, the unamplified regions (rps432-trnL [UAA], trnT [GGU]-trnE [UUC], atpB-rbcL), and the amplified regions with low variation levels (trnT [UGU] - trnL [UAA], rpl32F - trnL [UAG] y trnT [UGU] - trnL [UAA]) all suggest that in P. culiacanum, as in other parasitic plants, the size of the chloroplast was reduced, in comparison with the original size in tobacco plants (1560 nucleotides) (^{Bungard, 2004}).

Variability analysis

The seven analyzed populations are polymorphic, with 11 different haplotypes (2-4 per population), without taking indels into consideration as informative sites (Figure 1).

Figure 1 Haplotype distribution (without taking indels into consideration) of Pholisma culiacanum. Altata (AL), Culiacán (CU), Tosalibampo (TO), San Blas (SB), Heraclio SB (HS), Masiac (MA), and Cerro de Álamos (CA). 

Nucleotide total diversity (π) -without taking indels into consideration- was 0.2737 ± 0.0242; Masiaca and Tosalibampo had the highest and lowest values (0.2990, and 0.0007, respectively). Overall haplotype diversity -without taking indels (Hd) into consideration- was 0.855 ± 0.024; Altata and Cerro de Álamos shared the highest value (0.7330), while Tosalibampo had the lowest (0.4670) (Table 1).

First value (0.1349, p≤0.05) indicates a moderate level of genetic differentiation among studied populations. AMOVA showed that most of the observed genetic variation is mainly due to differences within populations (86.51 %), rather than differences between populations (13.49 %). This result is consistent with those of other parasitic plants studies (^{Jerome et al., 2002}; ^{Amico et al., 2009}) (Table 2). This was to be expected, because the populations’ Hd values are over 0.5 (except in Tosalibampo); high levels of molecular variance within populations result in high haplotipic diversity. High variation within populations can be interpreted as a consequence of gene flow between individuals of the same population (^{García-Franco et al., 1998}). This contrasts with other studies on parasitic and dependent plants, in which most genetic variation was found at population level. These species are characterized by: fragmented or patchy distribution; dependence on their hosts (^{de Vega et al., 2008}); the life history of both the parasite and the host; suffering geographic and stochastic events (bottle necks and founder events). However, there are species that are wind-pollinated or that use the wind to spread their seeds (such as P. culiacanum). These species usually have a lower genetic diversity index between populations (^{Allendorf et al., 2014}).

Although population structure had low levels of genetic differentiation among populations, the Pearson correlation coefficient and the linear regression model for genetic distances showed no evidence of isolation by distance (p>0.05 y r²=0.0026), unlike in other parasitic plants (^{Holzapfel et al., 2002}; ^{Jerome et al., 2002a}; ^{Mutikainen et al., 2002}; ^{deVega et al., 2008}). This allowed us to infer that there is a gene flow amongst the studied populations and that overall they behave as a single evolutionary population unit (^{García-Franco et al., 1998}; ^{Allendorf et al., 2014}), which is large enough to avoid the effects of genetic drift. However, in rare species -such as P. culiacanum-, where variation is equal or larger than in common species, evaluating other non-ecologic factors -such as life cycle strategies, reproductive attributes, demography, interaction with other organisms (the host in this case), and others factors- is advised (Eguiarte et al. 2007).

Haplotype distribution

Haplotype 5 had the highest frequency (19), and it was found in Tosalibampo, Masiaca, San Blas and Cerro de Álamos; haplotype 6 had the highest distribution rate, and it was found in Culiacán, San Blas, Heraclio Sb, Masiaca and Cerro de Álamos (Figure 2).

Figure 2 Haplotype network without taking indels into consideration, generated with Networks software (Bandelt et al., 1998). The Hs and the adjacent numbers indicate the haplotypes that were found. 

Populations with unique haplotypes were Altata (haplotype 1), Culiacán (haplotypes 3, 9, and 11), Tosalibampo (haplotype 10), and Cerro de Álamos (haplotypes 2 and 8). The other three populations (Masiaca, San Blas, and Heraclio Sb) were made up entirely by haplotypes present in other populations (haplotypes 4, 5, 6, and 7). Unshared haplotypes can be the result of mutations related to the colonization of new ecological niches (hosts, ecosystems, or both). This P. culiacanum ability to parasitize different plant genera in different ecosystems was documented by ^{Yatskievych (1985)}. It can also be a result of the expansion of hosts to new territories and the subsequent colonization of the parasite as proposed by ^{deVega et al. (2008)}.

Haplotype networks

In Figure 2, each haplotype is represented by a black circle the size of which is proportional to the frequency of the sequences that it is made up by. Mean vectors are represented in white and were generated by the software to link distant haplotypes. Mean vectors can be mutational steps, extinct, or not sampled subspecies. The analysis detected three major haplotype groups: 1) haplotypes 1 and 4, which are 51 mutational steps away from Lennoa (H12); 2) a group formed by haplotypes 2, 3, 5, 9, 10 and 11 is 160 steps away ; and 3) another group formed by haplotypes 6, 7 and 8 is 207 steps away. The methods used in haplotype networks are designed to model the evolutionary past via a stochastic process (coalescence), which is based on the idea that the alleles in a population can be traced back in time to the point where they coalesce with an ancestral allele. Therefore, the haplotypes furthest from the root sequence (6, 7, 8, 3, 11 y 9) are the most recent and the haplotypes closest to the root sequence (1 and 4) are the oldest (Eguiarte et al., 2007).

Haplotype phylogeny

The phylogenetic tree of the 11 haplotypes showed four groups: 1) the outgroup group (Lennoa); 2) haplotypes 1 and 4, as the basal group of P. culiacanum and the closest to Lennoa; 3) a branch which bifurcates into two major groups (80 % support), on the one hand, haplotypes close to the basal group 6, 8 and 7 (100 % support); and (4) a branch (100 % support) with haplotypes 2, 3, 5, 9, 10, 11. The general distribution matched the haplotype network diagram, including the proximity or distance from the root sequences (Figure 3).

Figure 3 Pholisma culiacanum haplotype phylogeny, without taking indels into consideration and using the Neighbor Joining method with Mega software (^{Tamura et al. 2007}). Support bootstrap values used 1000 iterations. 

Populations’ morphology varies geographically. There were two major differentiated groups. The first one includes the populations of Altata and Heraclio Sb, and the second is made up by the populations of Culiacán, Tosalibampo, San Blas, Masiaca, and Cerro de Álamos. The former has characteristic thin or branched stems originating from a single haustorium, similar to those of Lennoa. The second group has one single thick stem, over 1 cm in diameter. These observed characteristics matched the molecular evidence that resulted from the phylogenetic analysis and the haplotype network, which located most of group 1 haplotypes far apart from group 2 haplotypes.

The haplotype networks showed several mutational steps between those haplotypes that are mostly formed by individuals from Altata and Heraclio Sb (haplotypes 1, 4, 6, 7 and 8) and those from all other populations. Haplotype phylogeny also located these morphologically different populations’ haplotypes close to each other and to Lennoa. Other relevant characteristics were the different kind of soils in which the populations were found, ranging from coastal dunes (Altata) and sandy soil (Heraclio Sb), to clay soil in deciduous forests and the ecosystems (thorn forest, xeric shrublands, deciduous and semideciduous forests). Pholisma arenarium and Pholisma sonorae share these types of soil (Nabhan, 1979; Felger, 1980; ^{Yatskievych, 1985}; ^{Yatskievych and Mason, 1986}; ^{AGFD Plant Abstract, 2004}), which might suggests that these P. culiacanum populations belong to a different species of the Pholisma genus or to a P. culiacanum subspecies. For this reason, a comparative molecular analysis of the three species of the genus is necessary to determine the reason for these differences.

Given the high levels of haplotype genetic diversity, it is advisable to establish genetic resources conservation areas in the Culiacán, Cerro de Álamos and Altata ranges. These would include almost all the genetic diversity of the studied populations and most of the unique haplotypes would be found in the preserved areas. This is important because the aforementioned ranges are the most affected by anthropogenic activities, such as urban and touristic developments, agriculture, forestry, and land-use changes.

The Tosalibampo range is recommended for in situ exploitation, because the genetic diversity and haplotype indexes were the lowest; only two haplotypes were found and it had the most abundant population.