Microscale Adaptive Response of Charophytes of the Negev Desert, Israel: Species Divergences by AFLP

Genetic diversity of Chara populations found in Israeli Negev Desert streams was analyzed using amplified fragment length polymorphism (AFLP) in 31 individuals. The adjacent streams Avdat and Aqev were used as study sites in the Ein Avdat National Park. A Jaccard phylogenetic tree, based on 468 loci, showed clear separation between four Chara species: C. vulgaris, C. contraria, C. gymnophylla and Chara sp. The last two species served as out groups from the Northern Israel. In both C. vulgaris and C. contraria, genetic differences were found between the populations originating in the two streams: higher private loci number (80 vs. 24 in C. vulgaris and 58 vs. 25 in C. contraria) and polymorphic loci level (45% vs. 23% in C. vulgaris and 39% vs. 21% in C. contraria) were found in Avdat compared to Aqev. The genetic divergence revealed within and between the two streams is presumably adaptive and determined by natural selection associated with ecological stress. Sunlight intensity, water level and pH were found to be the main ecological variables associated with species clustering, through selection, in the varying ecology of the sampling stations along each of the streams. Evidence for gene flow between and within the streams was found using the structure analysis, suggesting that sampling sites condition is the regulating factor for oospores establishing, and hence, gene flow occurs more often between sampling sites with similar ecological conditions.


Introduction
The submerged green algae Charophytes are globally distributed, found predominately in standing and running water [1]. They are considered the closest living relatives to land plants [2], and include living algae in the family Characeae, and a number of extinct fossil families, mainly recognized by Lime-shells (gyrogonites) which are produced around the oospores [3]. Chara species' identification has been controversial among taxonomists, who must contend with limited morphological characteristics, intermediate species, and the unclear extent to which habitat conditions or developmental levels lead to morphological variation [4]. Chara vulgaris and Chara contraria are considered by most authors as separate species of the Chara genus [5]. Wood [6] characterized C. contraria as a conspecific variety of C. vulgaris, due to their similar morphology. C. vulgaris and C. contraria differ mainly in their cortex arrangement and spine cell position, where C. vulgaris cortex is defined as aulacanthous ('the secondary rows are more prominent, the spine-cells appear to sit in furrows') and C. contraria as tylacanthous ('the primary rows are prominent, spine cells appear to sit on the ridges', Bryant and Stewart 2002, p. 595) [7]. Both species are characterized as haplobiontic, monoecious, and reproduce mostly by selfing. Consequently, these species exhibit extremely limited genetic variation, compared to diploid forms. Nevertheless, the variation in their chromosome number (ploidy level) has been suggested to contribute to their genetic diversity [5]. Experimental crosses between C. vulgaris and C. contraria clearly indicated that they were reproductively isolated, although their "bridger" form could fertilize oospores that failed to germinate [5].
Molecular genetic studies using AFLP technology Vos [8] improved Chara taxonomy, which was based on Wood's definitions [6]. Mannschreck [9] revealed that C. intermedia and C. hispida could be regarded as separate species, in contrast to Wood [6], who classified C. intermedia as a variety of C. hispida. Boegle [1] found that C. intermedia and C. baltica, though morphologically similar, differ in salinity preferences and are indeed genetically different. In another study, C. curta and C. aspera were found to be morphologically different, although no genetic separation was found between them; as a result, the two species were treated synonymously. In this case, phenotypic plasticity resulted from variability in environmental conditions [4].
Studies of C. vulgrais and C. contraria in Israel were first described by Grant and Proctor [5]. Charophytes were recently found in central and northern Israel, mainly in the Upper Jordan and in the Oren River basin [10][11][12]. The first Chara species in Israel was found in the Hadera River as well as near the Dead Sea in Wadi Zerka and Ghuweira habitats [13]; unfortunately, these habitats were destroyed.
The current study focused on the genetic diversity among the Negev Desert Chara populations by AFLP markers. The analysis included mainly C. vulgaris and C. contraria species adaptation caused by natural selection associated with ecological local stresses of sunlight intensity and water level.

Description of the Study Site
A The Negev is a desert region in southern Israel. Geographically, it covers over 13,000 km², forming an inverted triangle whose western side is contiguous with the Sinai peninsula desert, and whose eastern and southern borders are the Arava valley and the Gulf of Eilat, respectively. The Negev Desert is arid, with an annual average precipitation of 25-200 mm [14]. Most of the rain water percolates through the ground and is stored as groundwater, throughout the Negev. One of the largest groundwater aquifers is found in the Negev mountain area, and contains brackish water used for agriculture. Surface runoff, occasionally in the form of floods that wash away soil and rocks, is commonly utilized for irrigation. The Zin stream is one of the largest in the Negev. Measured from the northern edge of Makhtesh Ramon in the central Negev to the Sodom south of the Dead Sea, the stream is 110 km long. Its drainage basin covers 1,400 km2. It crosses several areas with various topography, geomorphology, and geological substrates, which include the Ramon ridge, Avdat highland, Zin valley, Zinim ridge, the north wilderness, and Sodom wilderness. About a third of the Negev seasonal springs are usually dry and located in the Zin basin. Springs flowing from the rim of the Avdat highland are En Kadish and En Koriate (Kadesh-barnea) in the west, and En Avdat, En Aqev, and En Zik in the north ( Figure 1). These springs flow from the cliff of Avdat to the Zin basin.  [14].
The Avdat (30º49'54''N, 34º46'12''E) and the Aqev (30º48'57''N, 34º48'47''E) streams are usually characterized as a unique ecology habitat which can be distinguish by water level, sunlight intensity etc., although some of the sampling sites, especially adjacent ones, may have similar characters. Some of the sampling stations in Avdat from which samples were taken are relatively deep (2-4 m), some are intermediate (0.5-1.0 m) and some are shallow (~10 cm in some cases); all sampled stations in the Aqev stream were below 1m deep. Charophytes in the Dan River (33º14'28''N; 35º44'21''E) and the El Qinia pool in Northern Israel (33º13'47''N; 35º43'48''E) where taken for comparative studies.

Ramat Avdat Geology
Ramat Avdat is a table mountain stage located between Sede-Boqer and Makhtesh Ramon. It is a flat highland elevated about 800 m (average) above sea level. At its north-east it borders with the Zinim cliff, which rises 150 m (average) above sea level. The north-east of Ramat Avdat is drained by two main rivers: Nahal Nizana which flows to the west and Nahal Zin which flows to the east. Israel's watershed line passes between these two river systems. By the end of the Eocene epoch, the entire area of today's Israel was uplifted as a single plate [14]. Consequently, the ancient Tethys Sea receded, and the rocky layers at its seabed became exposed to erosion. In the Negev Mountain, the horizontal layers from the Eocene epoch covered most of the area. During the Neogene period, river systems were developed on the vast plains created by these layers, which led to depositional basins. Evidence for the vast extent of these rivers can be found in existing Neogene conglomerate outcrops in the area. The formation of the Arava rift valley, in the Neogene period, created elevation differences of hundreds of meters between the fault and the eastbound flowing rivers. This change diverted the rivers' flow to the Arava. In the Sede Boqer region, the river incised intensively while removing the Neogene conglomerate that covered the area, and later, also removing the Eocene limestone and chalk layers, and exposed the soft marl rocks underneath it. The river widened its channel over these rocks, and formed a wide valley, delimited at its northern side by the slopes of the Halukim anticline and the Hatzer anticline, and at its southern side, by the Zinim cliff, for 20 km. Today the river flows at a level about 90 m lower than its ancient level. The cliffs surrounding the Avdat highland were formed by the intensive incision of the rivers in the soft marl (Pliocene) that constitutes the basis of the plateau. The hard limestone and chalk rocks laid over them (Eocene epoch) stop this erosion process and form the cliffs. Judging by the concentration of travertine and from remnants of prehistoric sites in the area, it is estimated that the canyon in its current form was created about 45,000 years ago.

Sampling and Cleaning
Twenty nine Chara individuals (C. vulgaris, C. contraria and Chara sp.) from Avdat and Aqev streams were sampled. The Avdat stream individuals were sampled from six sampling sites along 2 km; in the Aqev stream, individuals were sampled from five sampling sites along 0.5 km ( Figure  1).
Two additional individuals of C. gymnophylla were sampled from the Dan River and the El Qinia pool in Northern Israel, and were used as an out group for comparison. In each of the sampling sites, one or two individuals were collected from a Chara patch of 0.20 m2, and each sample was examined separately. Distance between patches at the same stream was 5-2000 m.
Thalluses were cleaned with distilled water, epiphytes (such as filamentous green algae and diatoms) were mechanically removed under microscopic observation (x40 and x100), and the samples were kept in -80°C for further investigation. Taxonomy comes from Cavalier-Smith and algaebase.org [15].

Environmental Variables Measurements
Temperature, conductivity, and pH were measured mainly during 2008-2010, at least four times in a year (in four seasons) in each of the stream stations. The samples for chemical analyses of major elements (Cl -, Na + , Ca +2 , Mg +2 and HCO 3-, PO4 -3 , NO 3-, and Fe 2 O 3 ) were sampled twice a year (summer/winter). Chemical analyses were doing in the Newe Yar laboratory. Sunlight intensity was measured in the sites in klux, in order to characterize sunlight differences among the stations in each stream and between the streams.

DNA Extraction
DNA was extracted from frozen material that was fractured in liquid nitrogen (LN2). Genomic DNA was isolated using the DNeasy Plant Mini Kit (Qiagen, Hilden, Germany).

AFLP Analysis
AFLP analysis was performed according to Vos [8] with slight modification for a non-radioactive method as described by Satish [16]. Selective amplification was carried out using EcoRI and MseI primers, each with three additional selective nucleotides. Polymerase chain reaction (PCR) was performed using MseI primers and an EcoRI primer with the following selective nucleotides: ACA, ACG, ACT, or ACC, which were labeled by Fam, Ned, Vic, and Pet. The unlabeled MseI primers included the following selective nucleotides: CAC, CTT, CAA, and CTG were used as the second primer in the PCR reaction. Altogether, 12 primer combinations were used for PCR, and the most informative primer combinations were selected for further analysis.
Fragment analysis was performed using automated sequencer (ABI3130XL; Applied Biosystem, Foster City, STATE, USA). PCR amplification products (1ul) were added to 0.35 µl internal size standard (Liz-600) and 8.5 µl formamide for the polyacrylamide gel electrophoresis. Semi-automated AFLP fragment analysis was performed with GeneMapper® DNA fragment analysis software version 4.0 (Applied Biosystem, Foster City, STATE, USA). Fragments of 70-500 bp were scored for presence (1) or absence (0) and converted to a binary matrix. Some of the accessions were analyzed twice in order to test reproducibility.

Statistical Analysis
Genetic similarity coefficients among individuals were calculated from the binary data using the Jaccard similarity coefficient sij, Jaccard = n11/(n11 + n01 + n10), where nxy is the number of characteristics that have state x in individual i and state y in individual j. Possible character states are band presence (recorded as 1), band absence (recorded as 0) and missing data (recorded as ?). In order to convert for a dissimilarity coefficient we use D =1sij, Jaccard. A phylogenetic tree was generated from the above matrix using the neighbor-joining (NJ) algorithm [17]. Bootstrap analysis was conducted using 10,000 replicates as applied in FAMD software v.1.25. Analysis of molecular variance (AMOVA), PhiPT value for genetic variability, Mantel test, number of different bands and private bands, percentage of polymorphic loci (%P), and unbiased diversity (uh) were calculated for each population using GenAlEx v.6.4 Software [18]. Candidate loci (AFLP markers) for selection were analyzed by MatSAM v.2 to test their relation with ecological variables [19].

Environmental Variables
Environmental variables were measured during 2008-2010 in sampling sites of both streams. Fluctuations in water elements concentration ( Figure 2) and water characteristics (Figure 3) are represented by the error bars. Micro elements concentration in both streams were NO3-<0.5, Fe2O3 <0.3, PO4-3<1.0 mg/l (averaged all measurements). Sunlight illumination (klux), which was measured along of each stream, has shown differences in illumination, mainly at the Avdat sampling sites. Therefore, Avdat stream might divide into sunny sites (those most exposed to sunlight), and shaded. In the Aqev stream most of the sampling sites were defined as sunny sites. Sunlight values ranged between 50-100 klux in the sunny sites, and 3-10 klux in the shaded, during daytime. The water level varied mainly in the Avdat sampling sites.

Fragment Analysis (AFLP) of Chara Spesies
Four different species were included in the analysis: 19 individuals from Avdat and Aqev streams of C. vulgaris Linnaeus, nine individuals of C. contraria A. Braun ex Kü tzing, two individuals of C. gymnophylla A. Braun from the Dan River and one of Chara sp. from El Qinia pool. A total of 468 AFLPs loci were scored after amplification by six primer combinations in the four species. The number of loci varied among the species, and was higher in C. vulgaris (308) and C. contraria (272) than in the out groups C. gymnophylla (190) and Chara sp. (109). The number of private loci was also significantly higher in C. vulgaris (57) and C. contraria (61), than in C. gymnophilla (27) and Chara sp. (33) (Figure 4). A total of 468 AFLPs loci were scored after amplification by six primer combinations in the four species. The number of loci varied among the species, and was higher in C. vulgaris (308) and C. contraria (272) than in the out groups C. gymnophylla (190) and Chara sp. (109). The number of private loci was also significantly higher in C. vulgaris (57) and C. contraria (61), than in C. gymnophilla (27) and Chara sp. (33) (Figure 4).
Analyses of intra-species diversity of C. vulgaris and C. contraria revealed significant differences between Avdat and Aqev streams. Greater differences were found in Avdat compared to Aqev in: 1) polymorphic loci (% P) ( Figure 5A), and 2) unbiased diversity (uh) ( Figure 5B) parameters in both Chara species. Indeed, the number of private loci found is higher in Avdat compared to Aqev, again, in both species ( Figure 6). No significant differences were found in the overall polymorphic loci (%P) between C. vulgaris (55%) compared with C. contraria (52%).  Interestingly, as found by the AMOVA test, the variance of C. vulgaris among the two streams was relatively small (12.9%), while the variance within streams was significantly higher (87.1%). A similar result was found for C. contraria (14.3% among streams, 85.7% within). When testing C. vulgaris along Avdat's sampling sites, the variance was 40.1 % among sites and 59.9% within them. Because of sample size limitations, differences along sampling sites of C. vulgaris in Aqev and of C. contraria in either of the streams were not analyzed ( Table 1).
The correlation between the stream distance and Chara genetic distance was analyzed using a Mantel test. When considering both Chara species and both streams, the correlation was found to be R = 0.45 (p < 0.0001); higher correlation was found when Aqev C. vulgaris and Aqev C. contraria were considered separately (R = 0.62, 0.66, respectively), and smaller correlations for each of the two species in Avdat (R = 0.4, 0.21 for C. vulgaris and C. contraria, respectively) ( Table 2).  The highest number of candidate loci (AFLP markers) for selection in C. vulgaris was related to sunlight radiation, followed by water level and pH. Three candidate AFLP markers of C. vulgaris were found mostly in relatively low pH stations in both streams (pH variable). Four candidate markers of C. vulgaris were revealed only in the Avdat stream (water level variable). Finally, seven candidate markers of C. vulgaris did not show any stream or stations preference (sunlight intensity variable) (Table 3).

Phylogenetic Tree and Species Similarity (Jaccard Coefficient)
Species separation was found in the NJ tree, based on the Jaccard coefficient. C. gymnophylla and C. vulgaris were found to be the closest species compared to the other pairs. Most of the C. vulgaris and C. contraria individuals were clustered according to streams and sampling sites locations, although several were clustered differently from their sites. Avdat's C. vulgaris individuals were separated into three clusters, which were generally sites-related. The sampling sites differed ecologically in water level and sunlight intensity. A single Aqev's C. vulgaris individual (3-1b) was clustered together with the Avdat's individuals. Aqev's C. vulgaris individuals were mostly separated, indicating a higher genetic distance among them, especially of individual 1-2b. Relatively high genetic distance was found among Avdat C. contraria individuals, compared to Aqev C. contraria individuals (Figure 7). Indeed, the Avdat C. contraria individual (8a-1) was clustered next to the Aqev individuals in the NJ tree, and also according to Inferred ancestry analysis (Structure software). It was clearly shown that C. contraria individual (8a-1) which sampled from Avdat indicated higher probability that ancestor in the last three generations related to Chara population from Aqev stream (probability of Avdat is 0.013; probability of Aqev in the last three generations are respectively 0.187, 0.373, 0.427. Ln Prob of Data = -632.5, Mean value of ln likelihood = -558.9, Variance of ln likelihood = 147). The highest pair-wise genetic similarity (Jaccard coefficient) in Chara was found between C. gymnophylla and C. vulgaris (0.36). Other pair-wise comparisons showed lower similarity, in the range of 0.18-0.20 (Table 4).

Morphology and Taxonomy
The results of the current analysis show that AFLP separated the studied Chara into four species, as shown by the phylogenetic tree ( Figure 7). These species were hitherto delimitated morphologically. C. vulgaris and C. contraria, sampled in the Negev Desert region which is characterize by high sunlight intensity, are morphologically defined as corticated branchlets. In contrast, C. gymnophylla, sampled in the northern Israeli Dan River which receives significantly lower illumination, is characterized as ecorticate branchlets. Our results are supported by Grant and Proctor [5] who have shown that corticated branchlets configuration is heritable and characterize species that grow in daylight and higher sunlight intensity.
A higher genetic similarity (Jaccard coefficient) was found between C. gymnophylla and C. vulgaris (0.36), as expected from relative species [5]. The other pair-wise similarities were smaller, in the range of 0.18-0.20, supporting the morphological observation showing that the species are separated (Table 4). Furthermore, private loci identified in each of the species may be used as markers for species delimitation. The total loci number was higher in C. vulgaris and in C. contraria than in the out groups C. gymnophylla and Chara sp., possibly related to the significant differences in sample size; particularly as C. gymnophylla includes only two individuals and Chara sp. only one, which was mostly used as an out group for experimental control (Figure 4).

Genetic Diversity of C. vulgaris and C. contraria in Avdat and Aqev Streams
Chara vulgaris was found to be more common in the streams compared to C. contraria; some sampling sites contained only one of the species, and other sites contained both. We have not found any evidence for interbreeding between species via "bridger" species, as described by Grant and Proctor [5], but this phenomenon may not be ruled out.
A relatively lower overall polymorphic loci (%P) level was found in C. vulgaris (55%) and C. contraria (52%) species compared to other Chara AFLP studies, which showed significant higher polymorphic loci (%P) level (e.g., C. baltica (99%) and C. hispida (91%) etc. [20] or C. intermedia (83.3%) and C. tomentosa (91.4%) etc. [9]. We suggest that these differences in polymorphic loci (%P) level are consequence of the extreme environmental conditions of the Negev Desert streams, compared to the regions in above mentioned studies. Our findings show that the local ecological conditions, characterizing each sampling sites, in both streams, are the dominant factor behind Chara adaptive genetic divergence, presumably caused by natural selection. The most important factors appear to be sunlight intensity, water level (which is correlated with light intensity), and pH. The greater environmental variety in Chara sampling sites found in Avdat, compared with Aqev stream, resulted in differences of polymorphic loci (%P) level which was significantly higher in Avdat compared with Aqev in both species, respectively (45% vs. 23% in C. vulgaris and 39% vs. 21% in C. contraria) (  (Table 1). These findings indicate that the geographic separation between Avdat and Aqev streams plays a minor role in the intra-species genetic differences found in C. vulgaris and C. contraria, while the role of sampling sites ecology, which has been characterized, is the main cause for their relative higher variance.
Candidate loci (AFLP markers) for selection [19] (Joost et al. 2007) were found to be mostly associated with sunlight intensity, water level and pH among all of the ecology variables considered (Table 3), a finding which further supports that natural selection might play an important role in the adaptive genetic divergence of Chara through local ecological stresses in stream sampling sites. Using a Mantel test, we have found moderate correlation (R = 0.45, P<0.001) between stream distance and Chara genetic distance (Table  3). The phylogenetic tree ( Figure 7) and structure analysis ( Table 2) clearly showed that most Chara individuals in both species were clustered according to their relative location in the stream and according to ecological variables. Relationships between habitat ecology and AFLP markers in Chara were found in several studies. Boegle [1] found that C. intermedia and C. baltica, which showed similar morphology, were established in different salinity regions. Consequently, the salinity differences may be associated with their genetic differences. In addition, a recent study of Boegle [1] suggested that morphological differences among similar genetic species might be a result of biotic and abiotic factors (e.g., differences in sediment nutrients, light viability, or grazing pressure).
Since only a few kilometers separate the streams, gene flow might occur by animals, such as ducks and birds, which are common in the region and are known as oospore transporters [21]. Both structure analysis (Structure 2.3.3) and the phylogenetic tree ( Figure 7) provides some evidence for possible gene flow between the streams. Grant and Proctor [5] suggested that Charophytes are usually closed to new migration, especially when the existing population is established in the area and the incoming oospores might compete with the local ones. We suggest that site ecology is the major limiting factor for oospores establishing in a new site, meaning that gene flow occurs more often in sites with similar ecology. Therefore, the geographic separation between Avdat and Aqev streams is not a substantial factor.
The oospores bank in our study sites is stable, although flush floods often wash the inhabitants of the streams. Bonis [22] stress that higher density of oospores was found in the top layer (2-7.5 cm) of sediment surfaces in a temporarily flooded habitat, which represented annual life cycles. Indeed, higher variation in the density of buried oospores was found in Mediterranean temporary pools, which resulted mainly by oospores germination and their output [23]. One of the main characteristics of the desert climate is the winter flush floods. Ward and Blaustein [24] studied the influence of the flush floods on plants and animals in Aqev temporary pools, and found that the flush floods washed the inhabitants of pools; as a result, species richness decreases sharply, but later increases in pools which retain water.
The sample size in the study was relatively small (31 individuals) and was neither equally distributed between streams nor between species. This is due to limited Chara populations in the streams, and the sampling was the most comprehensive possible under the existing conditions.

Conclusion
The current study of the Israeli Chara species allows us to conclude that environmental factors may play the main role not only in species distribution over climatic differences in Israel, but also in microhabitats, as we revealed in the unique environment of the Negev Desert region. The microgeographical genetic polymorphism divergences identified in our samples as well as other studies such as in "Evolution Canyon" [25][26][27][28] is primarily determined by natural selection leading to adaptive complexes to the microscale ecological local stresses.
Charophyte flora of Israel, presently the most representative for the Eastern Mediterranean, has been occasionally studied since 1970 and more consistently in 2001-2013, revealing a diversity of about 20 species from 45 localities over six ecological regions, including the Arava Valley and Negev Desert. Hypsometrically, the lowermost locality occurs at 347 m b.s.l. in the Dead Sea rift [29]. Ecological characteristics and aquatic communities show that the major influence on charophyte distribution is climate, with sharp North -South gradients of temperature and humidity [29]. Our preliminary results of molecular genetic studies using AFLP analysis show the differences in %P level related to extreme environmental conditions in the Negev Desert compared to the Northern Israel, a typically Mediterranean region. The candidate loci (AFLP markers) for adaptive genetic divergence are those associated with sunlight intensity, water level, and pH, a finding consistent with the climatic differentiation and local environmental stresses playing a decisive role in the contemporaneous charophyte diversity.
In the next round of algological studies, the diverse methods of molecular genetics will be applied to systematics and ecological evaluation of charophyte communities in the Eastern Mediterranean. Figure 7. A phylogenetic consensus tree, built by Neighbor-Joining (NJ) algorithm using the standard Jaccard genetic distance, was tested by bootstraps (10,000 runs). The bootstrap showed 100% support for all the tested splits. Altogether 468 AFLP loci (six primer combinations) were analyzed in 31 Chara individuals