Expression and Molecular Evolution of Two DREB1 Genes in Black Poplar (Populus nigra)

Environmental stresses such as low temperature, drought, and high salinity significantly affect plant growth and yield. As selective forces, these adverse factors play essential roles in shaping phenotypic variation in plant populations. Black poplar (Populus nigra) is an economically and ecologically important forest tree species with widely distributed populations and is thus suitable for experiments detecting evolutionary footprints left by stress. Here, we performed expression and evolutionary analysis of two duplicated DREB A1-subgroup (DREB1) genes, PnDREB68 and PnDREB69, encoding transcription factors that are involved in stress responses. The two genes showed partially overlapping but distinct expression patterns in response to stresses. These genes were strongly and rapidly induced by cold stress in leaves, stems, and roots. In leaf tissue, dehydration stress induced the expression of PnDREB68 but not PnDREB69. PnDREB69 displayed more rapid responses and longer expression durations than PnDREB68 under salt and ABA stress, respectively. Based on single nucleotide polymorphism (SNP) analysis, we found significant population genetic differentiation, with a greater F ST value (0.09189) for PnDREB69 than for PnDREB68 (0.07743). Nucleotide diversity analysis revealed a two-fold higher πT for PnDREB68 than for PnDREB69 (0.00563 vs. 0.00243), reflecting strong purifying selection acting on the former. The results suggest that positive selection acted on PnDREB69, as evidenced by neutral testing using Tajima’s D statistic. The distinct selective forces to which each of the genes was subjected may be associated with expression divergence. Linkage disequilibrium (LD) was low for the sequenced region, with a higher level for PnDREB68 than for PnDREB69. Additionally, analysis of the relationship among carbon isotope ratios, SNP classes and gene expression, together with motif and domain analysis, suggested that 14 polymorphisms within the two genes may be candidates for an association study of important traits such as water use efficiency/drought tolerance in black poplar.


Introduction
The growth and development of plants are frequently challenged by environmental stresses such as drought, high salinity, and temperature change [1], resulting in a significant reduction in productivity in both crops and forest trees. These adverse factors represent major selective forces that act on plant phenotypes [2]. Plant species with widely distributed natural populations constantly exhibit phenotypic differentiation in many traits due to adaptation to their specific environments. This intraspecific variation is thought to be largely shaped by the forces of natural selection [3,4]. Thus, understanding the genetic basis underlying adaptive phenotypic variation is an important aim for both evolutionary genetics and molecular breeding in plants [5].
Stress tolerance is a quantitative trait that is often modulated by a group of transcription factor genes. CBF/DREB (C-repeat binding factor/dehydration-responsive element binding) is a gene family belonging to the AP2/ERF (APETALA2/EREBP) superfamily in plants, which encodes a large number of transcription factors [6]. CBF/DREB genes specifically bind to the C-repeat (CRT)/dehydration-responsive element (DRE) motif and regulate the expression of various stress-responsive genes that harbor these elements [7]. The DREB subfamily is composed of six subgroups of genes (A1-A6) [8]. Among these genes, in addition to containing a highly conserved AP2/ERF DNA-binding domain that is ubiquitous to all DREB genes, most A1-subgroup (DREB1) proteins possess conserved sequences including the nuclear localization signal (NLS) sequence PKRAGRTKFRETRHP, the DSAW motif, and the LWSY motif [9]. DREB1 genes, namely CBF1-CBF4, DDF1, and DDF2, were first identified in the model plant Arabidopsis thaliana [10][11][12][13]. Genes belonging to the DREB1 subgroup were subsequently identified and functionally characterized in many annual herbaceous plants, such as rice (Oryza sativa), wheat (Triticum aestivum), barley (Hordeum vulgare), maize (Zea mays), and Atriplex hortensis [7]. Regarding perennial woody plants, studies of DREB1 genes have been only documented in a limited number of species, including Populus spp [14,15], Eucalyptus spp [16,17], and peach (Prunus persica) [18,19].
The majority of studies of DREB genes have focused on the molecular functions of individual gene; an understanding of the mechanisms of functional divergence among DREB genes, and studies aimed at detecting the evolutionary footprints on these genes, are currently lacking. A previous study has proposed that the functional divergence of three duplicated members of the DREB1 gene family (CBF1, CBF2, and CBF3) may have resulted from specific sequence polymorphisms that were produced by contrasting evolutionary forces in Arabidopsis [20]. However, to date, the evolutionary history and molecular genetic analysis of the functional divergence of DREB1 genes have not been reported for any long-lived woody tree species.
Poplars are employed as a model woody species in plant biology studies, particularly studies involving adaptation, secondary growth, and the natural variation of traits [21]. Populus nigra L. (black poplar) belongs to section Aigeiros within the genus Populus (family Salicaceae) [22]. P. nigra is predominantly outbreeding, with a wide distribution across Europe, the Mediterranean coast of Northern Africa, and Western and Central Asia [23]. P. nigra has long been used to obtain interspecific hybrids (e.g., Populus. 6 euramericana) with favorable traits for commercialized planting by crossing this plant with related species, such as Populus deltoides. Natural populations of P. nigra also play important roles in the balance of riparian ecosystems by restoring riverbanks and controlling flooding [24]. Populations of P. nigra are distributed across a wide range of areas, from arid regions in Italy to relatively humid regions in Germany. This ecological feature makes P. nigra a suitable species to investigate the genetic factors that contribute to the natural variation of adaptation-related traits in plants such as drought tolerance. In this study, we combined expression analysis and molecular evolution estimation, including nucleotide diversity, neutrality testing, and linkage disequilibrium, using polymorphism data from two duplicated DREB1 genes, PnDREB68 and PnDREB69, to infer their diverged expression patterns in response to stresses and the selection forces that may have contributed to this divergence. We also examine the potential implications of our results for association studies of drought tolerance/water use efficiency traits.

Ethics Statement
No specific permissions were required for the described field studies. The Jade Spring Hill Nursery Garden in Beijing is not privately-owned.

Plant Materials
One-year-old P. nigra (clone 'N21') trees were grown in a greenhouse under controlled conditions at the Chinese Academy of Forestry (CAF) for gene expression analysis. To examine molecular evolution based on single nucleotide polymorphisms and phenotypic analysis, a total of 65 unrelated individuals representing the major distribution of P. nigra were planted in a field trial at the Jade Spring Hill Nursery Garden in Beijing. These individuals were initially sourced from natural populations across Europe and Asia, including four geographic regions (WE, Western Europe; CE, Central Europe; SE, Southern Europe; and CA, Central Asia) that cover 15 countries. Detailed information about each sample tree is presented in the Table S1. Young, healthy leaf tissues were sampled from the 65 individuals, quickly frozen in liquid nitrogen, and stored at 270uC until use.

Cloning of cDNAs and Amplification of their Genomic DNA
A genome-wide analysis of the AP2/ERF superfamily of poplar (Populus trichocarpa) revealed six gene members of the DREB1 subgroup [25]. Within this subgroup, two DREB1 genes, DREB68 (XM_002299529) and DREB69 (XM_002298031), are located on chromosome 1 and appear to be tandem duplicates arranged in a tail-to-tail manner. The two genes were not involved in the wholegenome duplication (WGD) event and the intergenic region between them is relatively short; thus, they are suitable for the analysis of expression divergence and molecular evolution. The counterparts in P. nigra were therefore named PnDREB68 and PnDREB69. Total RNA was extracted from leaf tissue using an Ambion Plant RNA Isolation Aid (Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions. Genespecific primers were designed for cDNA cloning via RT-PCR based on the sequences of DREB68 and DREB69 transcripts and the publicly available sequences of the two loci from JGI (http:// genome.jgi-psf.org/Poptr1_1/Poptr1_1.home.html; http://www. phytozome.net/search.php?method = Org_Ptrichocarpa) [26]. Genomic DNA for each individual tree was extracted from leaf tissue using a DNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions. Primer pairs for amplifying genomic DNA in P. nigra were designed according to the genome assembly for P. trichocarpa described above. The webbased program Primer 3 [27] was used to design primers for RT-PCR and genomic DNA amplifications. Expected PCR products were purified and subcloned into the pGEM-T vector (Promega, Fitchburg, WI, USA). To amplify genomic DNA, reactions were carried out in a volume of 50 ml containing 50 ng of template DNA, 0.2 mM of each dNTP, 0.5 mM of each primer, 16 PCR buffer plus MgCl 2 , and 0.5 U LA-Taq polymerase (TaKaRa, Dalian, China). All PCR reactions were run on a PE 9700 Thermal Cycler (Applied Biosystems, Foster City, CA, USA). Detailed information about the primers and PCR programs can be found in Table S2. Using a BigDye Terminator 3.1 Kit (Applied Biosystems), recombinant plasmids from RT-PCR were sequenced, and PCR products produced from genomic DNA templates were purified and subjected to direct sequencing. Sequencing reactions were conducted using an ABI 3730XL Automated DNA Sequencer (Applied Biosystems). Each amplicon was sequenced on both strands to generate a consensus sequence and to minimize PCR and sequencing errors.
Cis-acting elements within the predicted promoter sequences were deduced by searches in the PlantCARE and PLACE databases [28][29][30]. The deduced amino acid sequences of PnDREB68 and PnDREB69 were aligned using Clustal X [31] to perform a comparison analysis of gene structure.

Stress Treatments, qRT-PCR, and Microarray Analysis
Cuttings of P. nigra were cultivated in plastic pots containing an experiment-specific soil mixture in a greenhouse. When the plants were approximately 20 cm tall, they were exposed to cultivation solution containing 300 mM NaCl, 20% polyethylene glycol (PEG6000), or 100 mM abscisic acid (ABA) and grown at a low temperature (4uC) for 0 (normal growth conditions without stress), 2, 8, 24, or 48 h. Each treatment consisted of five replicates. Following the treatments, the leaves, stems, and roots were collected, frozen immediately in liquid nitrogen, and stored until use. Total RNA was extracted from the collected tissues using an Ambion Plant RNA Isolation Aid (Applied Biosystems) according to the manufacturer's instructions. The cDNA was synthesized using a PrimeScript RT Reagent Kit (TaKaRa). Then, qRT-PCR was performed in an ABI 7500 FAST Sequence Detector (Applied Biosystems) with SYBR Green Real-time PCR Master Mix (TaKaRa). Gene-specific primers were designed to amplify 120-130 bp fragments of PnDREB68 and PnDREB69, and parallel PCRs were carried out using a gene-specific primer pair for poplar ACTIN1 (GenBank Accession XM_002298674) as a reference gene. Primer sequences for the real-time PCR assay of PnDREB68, PnDREB69, and ACTIN1 are listed in Table S2. Four technical replicates for each PCR reaction were performed per RNA sample.
To perform comparisons of the orthologs within the genus Populus, expression profiles of DREB68 and DREB69 from P. trichocarpa in a range of samples were analyzed using the publicly available Affymetrix microarray data for poplar from the PopGenExpression server (http://www.barutoronto.ca) [32]. The dataset can also be downloaded from the Gene Expression Omnibus (http://ncbi.nlm.hih.gov/geo) under the accession number GSE13990. Probe sets for P. trichocarpa DREB68 and DREB69 gene models were retrieved using the Probe Match tool in the NetAffx Analysis Center (http://www.affymetrix.com/ analysis/index.affx).

Phylogenetic Analysis
To infer the general evolutionary history of plant DREB1 genes, amino acid sequences of representative members of the A1 group of DREB family (DREB1) proteins from five species (P. trichocarpa [25], Arabidopsis thaliana, Oryza Sativa, Sorghum bicolor, and Vitis vinifera) were randomly selected and retrieved from the wholegenome databases in Phytozome v9.1 at www.phytozome.net (Dataset S1). Thus, the sequences included in this analysis did not specifically consider genes created by WGD or tandem duplication. These sequences, along with those of PnDREB68 and PnDREB69, were subjected to phylogenetic analysis using MEGA 5 [33]. Sequences from members of the A2 group of DREB family (DREB2) proteins in moss (Physcomitrella patens) were also downloaded from Phytozome and used as an outgroup (Dataset S1). The neighbor-joining method was used to build phylogenetic trees with the p-distance amino acid substitution model. Phylogeny testing was performed using 1,000 bootstrap replicates.

Nucleotide Diversity and Population Differentiation
Genomic DNA regions containing PnDREB68 and PnDREB69 were sequenced in 65 unrelated individuals. Sets of alleles for each locus were aligned using Clustal W [31]. Single nucleotide polymorphisms (SNPs) were identified and counted by visually inspecting each sequence alignment for each locus. The corresponding chromatogram files were checked to confirm each polymorphic site and its genotype class. Nucleotide diversity was estimated using DnaSP v5.10 [34] (http://www.ub.es/dnasp/). The levels of nucleotide diversity of the candidate genes were estimated as hw and p, which are based on the number of segregating sites [35] and the average number of nucleotide differences per site between sequences [36], respectively. Insertions and deletions (indels) were excluded from these analyses. Population genetic differentiation was estimated as the fixation index, F ST [37,38] between pairs of populations (geographic regions in this study) using DnaSP. The significance of population differentiation was evaluated using the S nn -test statistic of Hudson [39].

Tests of Neutrality
To test for departures from the standard Wright-Fisher model of neutral evolution, statistics including Tajima's D [40] and Fu and Li's D [41] were calculated for both whole-gene and coding regions, synonymous (syn) sites, nonsynonymous (nonsyn) sites, and silent (synonymous and noncoding) sites using DnaSP. Fu and Li's D is believed to be more sensitive than Tajima's D for detecting evolutionary events such as population expansion and positive selection [41]. Neutral testing was performed using SNP data obtained from each geographic region or whole population.
The significance of these statistics was tested by 1,000 replicates of coalescent simulation based on Pi [42].

Linkage Disequilibrium
Pairwise linkage disequilibrium across the PnDREB68 and PnDREB69 loci was estimated as the squared correlation of allele frequencies, r 2 [43]. The value between pairs of SNPs (minor allele frequency .0.10) was measured using TASSEL software [44], with one-tailed Fisher's exact tests and Bonferroni corrections applied for each locus. The statistical significance of each value was determined using Fisher's exact tests.

Measurement of Carbon Isotope Ratio
Carbon isotope ratio (d 13 C, %), or carbon isotope composition, was proposed as an indicator of water use efficiency, which can reflect drought tolerance in plants [45]. Leaf tissues were collected from the 65 clones in September. This time point was chosen because during September, bud set was observed in all trees, which indicates the end of the growing season. The leaves were dried in an oven at 70uC and ground to a fine powder. The samples then were sent to the MASS Laboratory of Institute of Atom at the Chinese Academy of Agricultural Sciences, and the carbon isotope ratio of each sample ( 13 C/ 12 C sample ) was measured using a MAT-251 Isotope Ratio Mass Spectrometer (Finnigan, Waltham, MA, USA). These ratios were expressed relative to the V-PDB standard (R PDB ), where d 13 C (%) = ( 13 C/ 12 C sample 2 13 C/ 12 C PDB )/( 13 C/ 12 C PDB )61,000. Carbon isotope discrimination values were not shown, as the carbon signature of the trial field was not identified. For each treatment, the ninth leaf from each of three trees was collected and mixed for subsequent analysis. Thus, a total of nine leaves per individual were sampled and three replicated values were thereby obtained. Differences in carbon isotope ratios among 65 individuals, as well as among individuals for each SNP, were evaluated by one-way analysis of variance (ANOVA, a = 0.05). Statistical analyses were performed using Data Processing System (DPS) 3.01 software [46].

Results
Characterization of cDNA and Genomic DNA Sequences at the PnDREB68 and PnDREB69 Loci In P. trichocarpa, the two DREB1 genes examined are annotated in NCBI as DREB68 (XM_002299529) and DREB69 (XM_002298031), and the corresponding gene models in the genome database are annotated as Potri.001G110800 and Potri.001G110700, respectively. Using primers designed with these sequences, the cDNAs with complete coding sequences (CDSs) for the two genes were isolated from P. nigra using RT-PCR (GenBank accession numbers KJ000090 and KJ000091). These two genes were designated as PnDREB68 and PnDREB69 according to the nomenclature for the P. trichocarpa orthologs in NCBI. PnDREB68 and PnDREB69 have 753-bp and 693-bp open reading frames encoding polypeptides of 251 and 231 amino acid residues, respectively. We also isolated 111 bp and 93 bp of the 59UTR and 183 bp and 234 bp of the 39UTR for PnDREB68 and PnDREB69, respectively. These two genes share 57.32% nucleotide sequence and 54.75% amino acid sequence identities. A comparison of PnDREB68 and PnDREB69 with homologs in P. trichocarpa indicated homologous DNA sequence identities of 60.94% and 97.09%, and amino acid sequence homologies of 77.99% and 98.27%, against DREB68 and DREB69, respectively.
We obtained 5,348 bp of the genomic DNA sequence across the PnDREB68 and PnDREB69 loci through four rounds of PCR, each of which produced a fraction of the entire sequenced region ( Figure 1A). This genomic region was re-sequenced in 65 unrelated individuals (GenBank accession numbers KJ019894-KJ020023). The structure of the amplified region is illustrated in Figure 1A, showing a tail-to-tail orientation for transcription of the two genes, which is consistent with that in P. trichocarpa. Within this region, PnDREB69 is intronless, with a 696-bp CDS and a 996-bp 59 upstream sequence, while PnDREB68 possesses one intron (75 bp) that separated the CDS into two fragments (456 bp and 300 bp, respectively), and a 1,395 bp 59 upstream sequence. Regarding the 59-upstream sequences, we assigned the region upstream of the partial 59UTR (obtained from cloned cDNA) as the promoter for PnDREB69 (903 bp) and PnDREB68 (1284 bp), as the transcription start sites (TSS) for the two genes were not experimentally determined. A 1,430-bp intergenic sequence between the translation stop codons of the genes was also amplified ( Figure 1A). Two overlapping PCR primers (RP2 and FP3) (Table S2) covered a 45-bp sequence in the intergenic region ( Figure 1A). We separated the intergenic sequence in the middle of the region comprising the two primers and assigned the upstream and downstream sequences to the PnDREB69 and PnDREB68 loci, respectively, to perform molecular genetic analyses for each gene. The cDNA sequences were further confirmed by analyses of their genomic DNA sequences.
A 903-bp promoter sequence of PnDREB69 and PnDREB68 were compared to infer whether differences in regulatory elements could contribute to the different expression profiles (see below) between the two genes. Various cis-acting elements were detected, in which eleven were found to be shared by PnDREB69 and PnDREB68 with differences in copy number in most cases (Table  S3). Moreover, heat stress-responsive (HSE), drought-inducible (MBS), and auxin-responsive elements were detected in the promoter region of PnDREB68 but not in PnDREB69, while a salicylic acid-responsive (TCA) element was only observed in the PnDREB69 promoter. Alignment of amino acid sequences revealed conserved DREB1 signature sequences shared by both proteins, including a nuclear localization signal (NLS), an AP2 DNA-binding domain, and a DSAW motif ( Figure 1B). A LWSY motif was also identified in PnDREB68, but it was modified to LWSS in PnDREB69. A putative C-terminus activation domain similar to Arabidopsis AtDREB1B/CBF1 could also be observed for PnDREB69 and PnDREB68 ( Figure 1B).

Expression Patterns of PnDREB68 and PnDREB69 under Normal Growth Conditions and under Abiotic Stress
The expression patterns of the duplicated DREB68 and DREB69 genes were determined by quantitative RT-PCR under normal growth conditions, abiotic stress conditions, and ABA treatment. The expression of PnDREB68 and PnDREB69 was rapidly induced by cold stress in leaf and stem tissues (more than two-fold higher than control levels at the 2 h time point), and a slow response to cold stress was observed in roots (with higher relative abundance than the control beginning at the 8 h time point) (Figure 2A, 2B. and 2C). In most cases, PnDREB69 had higher expression levels than PnDREB68 when the plants were exposed to cold stress. When salt stress (NaCl) was applied, PnDREB69 and PnDREB68 displayed elevated expression levels in stems at the 8 and 48 h time point, respectively ( Figure 2B), but the expression levels were reduced in roots ( Figure 2C). Interestingly, in leaf tissue, dehydration (PEG6000) stress significantly induced the expression of PnDREB68 at the early stage of treatment (2 and 8 h) but not the expression of PnDREB69 (Figure 2A), whereas in stem tissue, both genes were slightly induced (with different relative abundances) by dehydration ( Figure 2B). ABA treatment induced the expression of PnDREB69 at all time points in leaves and stems, while PnDREB68 showed elevated expression only at a single time point (24 h in leaves and 48 h in stems) (Figure 2A,2B). Very few transcripts were detected in root tissues under dehydration or ABA treatment conditions ( Figure 2C), which suggested that the expression levels for both genes were reduced. The microarray analysis of orthologs from P. trichocarpa showed that the expression levels of DREB68 were higher than those of DREB69 across most samples, with the exception of xylem tissue. In mature leaf and seedlings grown in continuous darkness, the two genes exhibited inverse trend of transcript accumulations ( Figure S1).

Phylogenetic Relationships and Motif Composition
To investigate the evolutionary relatedness of duplicated PnDREB68/PnDREB69 with the A1 group of DREB (DREB1) genes in other plant species, a NJ phylogenetic tree was constructed based on 26 sequences from P. nigra and seven additional land plant species. The rooted phylogenetic tree (based on DREB1 proteins) shows that PnDREB69 and PnDREB68 have the closest relationships with the homologs in P. trichocarpa (DREB69 and DREB68) and Populus tomentosa (PtCBF2; Figure 3A). The phylogenetic tree forms two well-defined branches. One branch contains only sequences from dicots, while the other branch contains sequences from monocots, suggesting that the diversification of the DREB1 genes occurred before the divergence of dicots and monocots.
Twenty distinct motifs were detected in the 26 represented DREB1 proteins using the MEME program ( Figure 3B). Details of the 20 motifs were listed in Table S4. Common motif composition occurred in most of the closely related members, which suggested similar functional roles of the motifs among these proteins. Motifs 1, 2, and 3 were most commonly present in all tested proteins. Motif 8 was only found in dicots, while motif 19 was uniquely present in monocots. Additionally, motif 10 appeared to be specific to perennial woody plants (poplar and grape). The basal plant species moss had the lowest number of motifs [47], in which motif 17 was not found in other seed plant species.

Nucleotide Diversity and Neutrality Tests
Based on the alignment analysis of the 5,348-bp genomic DNA sequences from 64 unrelated individuals, fifty common SNPs (minor allele frequency .0.10) were detected, including three from coding regions and 47 from noncoding regions. Among these SNPs, 11 were identified within PnDREB69, while 39 (including three located in CDS) were detected from PnDREB68 ( Figure 1A). Two nonsynonymous SNPs were found in the coding region of PnDREB68, causing amino acid changes in the putative activation domain from Val and Ala to Leu and Thr, respectively ( Figure 1A). A total of 18 indels were revealed, in which 14, 2, and 2 were located in the PnDREB68 promoter, PnDREB69 promoter and intergenic region, respectively ( Figure 1A). Of these, eight polymorphic sites (four SNPs and four indels) occurred in several cis-acting elements, including the AT1-motif, Box 4, G-box, and TATA-box (Table 1).
Within the entire sequenced region, an analysis of the SNP data found that the promoters and intergenic region showed the highest level of nucleotide diversity relative to coding regions and intron ( Figure 1A, Figure 4A). When the two loci were considered separately, a higher average nucleotide diversity was observed for PnDREB68, as evidenced by two-fold greater values of h W (0.00777) and p T (0.00563), than those of PnDREB69 (h W = 0.00389, p T = 0.00243) ( Figure 4A, Table 2). The overall nucleotide diversity of samples from each geographic region of PnDREB68 was also found to be higher. Except for the CA samples, larger values of p sil , p syn , and p nonsyn for PnDREB68 from the three remaining geographic regions (WE, CE, and SE) were observed, indicating higher levels of nucleotide diversity at these sites. A greater value for the p nonsyn /p syn ratio was also inferred for PnDREB68 relative to that of PnDREB69 (0.54984 vs. 0.31333). However, compared with PnDREB68, samples from CA at

Population Differentiation
We estimated the genetic differentiation within pairs of P. nigra populations using the fixation index F ST . Both genes showed significant genetic differentiation among populations, with a S nn value of 0.36564 (P = 0.0150) and 0.38653 (P = 0.0030) for PnDREB68 and PnDREB69, respectively. No significant pairwise estimate of F ST was detected between populations. The results show that low to moderate values of F ST were detected for six different population pairs (Western-Southern, Western-Central, Western-Central Asia, Southern-Central, Southern-Central Asia, and Central-Central Asia), ranging from 20.01400 and 0.01953

Linkage Disequilibrium
The extent of linkage disequilibrium (LD) across the PnDREB69 and PnDREB68 loci, including the coding regions, intron, intergenic region, and promoters, were estimated using r 2 as the parameter. A total of 50 common SNPs (frequency .0.10) were included in this analysis, yielding 1,225 pairs of sites for the r 2 calculations. The average level of LD was low, with an r 2 value of 0.140251 for the entire region sequenced. In addition, r 2 values of 0.10641 and 0.21352 were estimated for the PnDREB69 and PnDREB68 loci, respectively. The architecture of LD is illustrated in Figure 5. The LD values varied across different regions, with higher values for 59UTR (average r 2 = 0.3031) and intergenic regions and lower values for exons (average r 2 = 0.0647) and introns (average r 2 = 0.0854). Four (from SNP site 3864 to 4378) and five (from SNP site 4424 to 4688) SNPs in high LD with each other were found to form two distinct haplotype blocks ( Figure 5). Within each haplotype block, LD among the SNPs was high (r 2 . 0.8279).

Identification of Candidate SNPs Affecting the Carbon Isotope Ratio in P. nigra
The expression of PnDREB68 was strongly induced by dehydration (PEG6000) stress in leaves, which was not observed for PnDREB69. PEG6000 treatment is frequently used to mimic drought stress [48]. Therefore, we further measured an important indicator for plant drought tolerance/water use efficiency (WUE), the carbon isotope ratio (d 13 C), to investigate the natural variation of this phenotype and its potential relationship with SNP polymorphisms. The carbon isotope ratio showed abundant phenotypic variation; the values ranged from 225.99 to 230.40% (mean value of 228.49%). ANOVA indicated a significant difference (P,0.05) in d 13 C values among the 65 individuals tested. This result suggests that the carbon isotope ratio trait bears significant natural variation in P. nigra populations and thus is suitable for population genetic analysis. We recorded the genotype (homozygous or heterozygous for each SNP) of each individual tree in our experimental panel (65 individuals), and the carbon isotope ratio of each individual was assigned to each genotyped SNP. For each SNP assayed, phenotypic data from trees with the same genotype were pooled and compared with each other. Significant differences in carbon isotope ratios were found among individuals for four SNPs (M645, M779, M374, and M519) ( Figure 6A). M645 and M519 showed three genotype classes whereas M779 and M374 displayed only two genotype classes. All of these SNPs are located in promoter region of PnDREB68. No correlation was found between carbon isotope ratios and the SNPs located in PnDREB69.
To further examine whether PnDREB68 expression in the different genotype classes of SNPs is correlated with carbon isotope ratio, expression levels of the four SNPs, M645, M779, M374, and M519, were compared using qRT-PCR. In this analysis, leaf tissues from 24 (8 for each of the three genotype classes of M645, M779 or M519, Figure 6A) or 16 (8 for each of the two genotype classes of M374, Figure 6A) individuals treated with PEG6000 for eight hours were used. These individuals were randomly selected from the 65 trees used for SNP identification and d 13 C measurement (Table S1). Among these SNPs examined, significant differences (P,0.01, ANOVA) in the expression levels among the genotype classes were found only in M519 ( Figure 6B). Individuals of AC group showed the highest PnDREB68 transcript levels, and this gene was expressed at the lowest level in the AA group. These data were negatively correlated with the d 13 C values observed for the genotype classes of M519 ( Figure 6A). In these genotype classes, lower (AC group) and greater (AA group) d 13 C corresponded with higher and lower expression of PnDREB68, respectively.

Discussion
Gene duplication is believed to be a source of evolutionary forces leading to functional diversification [49,50]. The consequences of retaining duplicated genes can include the following: (1) retaining the function of its ancestors; (2) loss of gene function (pseudogenization); (3) retaining part of the function of its ancestors (subfunctionalization), and (4) acquiring a new function (neofunctionalization) [51]. Tandem duplication refers to the occurrence of two sequences, one following the other, in the same chromosomal segment. The most remarkable difference in expression patterns revealed in this study is that the expression of PnDREB68 in leaves was dramatically triggered at the early stage of PEG6000 treatment, whereas no significant increase in expression of PnDREB69 was observed under the same stress conditions. Expression of PnDREB69 in leaves and stems was more sensitive to ABA than that of PnDREB68, as evidenced by higher expression of PnDREB69 relative to that of PnDREB68 throughout the duration of treatment. In leaf tissues, PnDREB68 may trigger a series of dehydration/drought responses by regulating the expression of downstream genes, while the ABA signal may be essential for the expression and function of PnDREB69 in both leaves and stems. This phenomenon may suggest functional divergence between the tandemly duplicated PnDREB68 and PnDREB69 genes involving in subfunctionalization and/or neofunctionalization, which need to be experimentally validated. Table 2. Nucleotide diversity and neutrality tests in PnDREB69 and PnDREB68 from P. nigra. The differences in the type and copy number of regulatory elements between the two genes suggested that these differences could be partly responsible for the differential expression of PnDREB69 and PnDREB68. Specifically, the PnDREB68 promoter possesses a drought-inducible (MBS) element that is absent in PnDREB69, which was consistent with the observation that the expression of PnDREB68 was strongly induced by drought stress (PEG6000) in leaf tissues ( Figure 2A). Moreover, the finding of sequence differences in the LWSY(S) motif and the putative activation domain also support the potential functional differences between PnDREB69 and PnDREB68. Expression divergence of orthologs from the close relative, P. trichocarpa, was also revealed by microarray analysis, as reflected by distinct expression patterns across a range of tissues, organs and treatments. The expression of CBF genes (DREB1 members) has been welldocumented in other plant species. In Populus balsamifera subsp. trichocarpa, four CBF genes (CBF1-CBF4) are cold-inducible, with various levels of induction observed among genes, depending on the tissue (leaf or stem) in which they are expressed [14]. Divergent expression was also observed among three CBF genes in wild tomatoes. These genes are all cold-induced, with various transcript accumulation levels observed among accessions, and they are drought-responsive in Solanum peruvianum, although CBF2 is not drought-responsive in Solanum chilense [52]. Genome-wide analysis of the expression pattern of DREB1 genes has also been reported in A. thaliana and rice. The Arabidopsis DREB1 family is composed of six gene members. Among these, CBF1-CBF3 are tandem duplicated and are involved in the cold-response pathway [53], while CBF4 regulates drought adaptation [11], and DDF1-DDF2 function in the salt response and in development [12,13]. In rice (Oryza sativa L.), six out of ten DREB1 genes are expressed in response to low temperatures (4uC or 12uC), and two of these genes are also induced by drought or salt stress [54]. Previous studies on other plant species and our data from P. nigra show that the response to cold stress is a common phenomenon in plant DREB1 genes, whereas other expression behaviors such as responses to drought, salt, and ABA are less common.
The phylogenic relationships and conserved motif analysis revealed diversification and a general evolutionary history of plant DREB1 proteins. Two distinct branches in the phylogenetic tree indicated that DREB1 genes have common ancestors before the split of dicots and monocots ( Figure 3A). MEME analysis revealed an increasing number of motifs from moss to seed plants, and several motifs are limited to specific species (dicots, monocots or woody plants) ( Figure 3B). These data suggests a possible divergence of functions among different clades of plants in which these motifs are involved. Motif 17 in moss is absent in other seed plants ( Figure 3B), which is probably because of a motif-loss event during the evolutionary history of DREB1 proteins.
Evaluations of the level of genetic differentiation can provide information about the extent of population structure, which is crucial for further genetic studies such as association mapping. Significant genetic differentiation in the natural populations was found for both loci (Table 3), where PnDREB69 displayed a higher fixation index (F ST ) value than PnDREB68, indicating a greater level of population differentiation in the former. The degree of genetic differentiation observed in this study was somewhat lower than that previously reported for P. nigra (overall F ST = 0.268 based on AFLP, and 0.081 based on microsatellite analysis) [24], while still confirming the existence of considerable population subdivision in this species. The estimated level of population differentiation in P. nigra is comparable to estimates from other close relative poplar species including Populus tremula (F ST = 0.117 for five loci) [55], P. balsamifera (F ST range, 0.018-0.256 for three loci) [56], and P. tomentosa (F ST range, 0.3552-0.12686 for two loci) [57].
Investigating intraspecific variation of nuclear genes (mainly inferred from the level of nucleotide diversity) in natural populations is essential for understanding the evolutionary forces that act on plant genomes. In this study, we evaluated the nucleotide diversity of two tandemly duplicated DREB1 genes. We detected two-fold higher total nucleotide diversity (p T ) for PnDREB68 than for PnDREB69 ( Table 2), indicating that more relaxed evolutionary constraints occur in PnDREB68 than in PnDREB69. This diversity originated from a much higher value of nucleotide diversity of silent sites (p sil ) for PnDREB68 than for PnDREB69. The higher level of nucleotide diversity of PnDREB68 may be a reflection of strong purifying selection acting on this gene. The levels of nucleotide diversity for the two DREB1 genes are lower than previous estimates for nine loci in P. nigra [58]. A comparison within the genus Populus revealed significantly lower nucleotide diversity in PnDREB loci than in two sucrose synthase genes (PtSUS1/PtSUS2, p T = 0.00924/0.01093) in P. tomentosa [57]. For each gene, the nucleotide diversity of PnDREB69 is similar to that of P. balsamifera (p T = 0.0025) [56] but lower than that of glycosyl-phosphatidylinositol-anchored protein gene (PtCOBL4) in P. tomentosa [59] and P. tremula (p T = 0.0042) [60]. However, PnDREB68 showed comparable levels of diversity relative to PtCOBL4 and to that in P. tremula but was two-fold higher than that in P. balsamifera. This interspecific variation in nucleotide diversity within the same genus can be explained by several factors, including differences in sampling procedures, the genomic regions studied, demographic histories, genetic backgrounds, and muta- tion rates. Variation in nucleotide diversity was also demonstrated in other DREB1 genes, including seven CBF genes in rye (Secale cereale L.; p ranging from 0.0015 to 0.0145) [61] and Arabidopsis thaliana, in which the p value of promoters for CBF genes ranges from 0.00109 to 0.02749, while the nucleotide diversity of transcript units ranges from 0.00264 to 0.00685 [20]. An analysis of orthologs of three CBF genes in wild tomatoes indicated that the average nucleotide diversity ranges from 0.021 to 0.079 in S. chilense and from 0.016 to 0.106 in S. peruvianum [52]. Tajima'D and Fu and Li's D tests are frequently used estimates that determine whether a gene has departed from neutral evolutionary expectations. The generally negative values for Tajima's D indicate an excess of rare variants in the two genes, which may have resulted from a recent population bottleneck (hitchhiking event) or a process such as background selection. Negative and not significant values for Tajima's D and Fu and Li's D were estimated for the entire region, coding sequence, synonymous sites, and nonsynonymous sites of PnDREB68, which is consistent with the results of nucleotide diversity analysis, suggesting that purifying selection was the main force driving the evolution of PnDREB68. By contrast, a much lower level of nucleotide diversity was detected for PnDREB69, with a smaller Figure 5. Pairwise linkage disequilibrium (LD) in sequenced region across the PnDREB69 and PnDREB68 loci. LD was estimated using the squared allelic correlation coefficient (r 2 ) between each pair of SNPs. Boxes below represent the genomic organization of the sequenced region as described in Figure 3. The physical position of each common SNP is labeled on the left and is linked to the corresponding location on the schematic diagram of the genomic region. Each grid in the upper-right matrix represents the level of LD reflected by r 2 for comparison of each SNP pair. doi:10.1371/journal.pone.0098334.g005 nonsynonymous/synonymous nucleotide substitutions (p nonsyn / p syn ) ratio than that observed for PnDREB68. Furthermore, a significant negative value for Tajima's D for this locus was inferred for the coding region and for nonsynonymous sites of PnDREB69, supporting the hypothesis that this gene has undergone positive selection.
The response to dehydration stress in leaf tissue in PnDREB68 was absent from the expression patterns observed for PnDREB69, which may be associated with positive selection acting on PnDREB69. The estimates of LD level for PnDREB68 or PnDREB69, and the entire sequenced region spanning the two loci, are consistent with our previously calculated LD values for nine genes in black poplar, revealing low LD levels in this species [58]. Compared with other species in the genus Populus, the strength of LD of black poplar is slightly higher than that of P. tremula (squared allelic correlation coefficient, r 2 ,0.1; ,200 bp for 77 genes) [60] and the gibberellin 20-oxidase gene (PtGA20Ox) in P. tomentosa (r 2 = 0.1 within a 500-bp region) [62] but much lower than that of P. balsamifera (r 2 = 0.52 for 372 gene fragments) [63]. This interspecific variation in LD levels may be partially attributed to the different gene regions considered, and it also indicates differences in recombination history among species within the genus Populus. Despite its high nucleotide diversity, PnDREB68 exhibited more extended LD than PnDREB69, indicating that more favorable mutations were fixed during the evolutionary history of PnDREB68 than that of PnDREB69. This distinct spectrum of LD extension also explains the expression divergence of the two duplicated DREB1 genes.
DREB family genes have recently been incorporated into association studies of stress-tolerance-related traits in crop species. Two and one nonsynonymous SNPs in ScCbf15 and ScCbf12, respectively, were found to be associated with frost tolerance in rye (Secale cereale L.) [64]. Lata et al. reported an association between a synonymous SNP in the DREB2 subgroup gene SiDREB2 and dehydration tolerance in foxtail millet (Setaria italica L.) [65]. However, molecular evolution and/or association studies of DREB family genes have not previously been documented in any woody tree species. In this study, preliminary analysis of phenotypic data and SNP genotypes revealed four noncoding SNPs that may have affected the carbon isotope ratio in the genetic materials studied. QRT-PCR analysis revealed significant differences in expression levels among the three genotype classes of M519, supporting the contribution of this SNP to the natural variation of carbon isotope ratio of black poplar. Based on the fact that plants can discriminate against 13 C during photosynthesis and the isotopic ratio of 13 C to 12 C (d 13 C) varies in C 3 plants, the important parameter related to drought tolerance, instantaneous water use efficiency (WUE), is found to be positively related to d 13 C [66,67]. This theory has been demonstrated in many plant species [68]. Our results suggest the potential importance of regulatory polymorphisms (located in gene promoters) in contributing to phenotypic variation [69]. Moreover, eight polymorphisms led to nucleotide changes in regulatory elements in PnDREB69 and PnDREB68, including the AT1-motif, Box4, G-box, and TATA-box, which are important for light-regulated gene expression. As cis-acting elements in promoters are crucial for gene expression and changes in their sequences could influence gene regulation, further studies need to be conducted on these polymorphisms to address their potential effects on gene expression and phenotypic variation. Two nonsynonymous SNPs with unknown functions were also found in the putative activation domain of PnDREB68. An activation domain is essential for the activity of a protein, and this has been demonstrated in Arabidopsis CBF1, whose C-terminal 98 amino acids activation domain functions in trans-activation. The 14 polymorphisms described above, therefore, are candidates for further association studies of important traits in P. nigra, such as water use efficiency/drought tolerance.
In summary, we presented expression and molecular evolution analysis of two duplicated DREB1 genes in black poplar, P. nigra. We detected expression divergence; the most remarkable difference was that the expression of PnDREB68 varied in response to dehydration stress, which was not observed in PnDREB69. This result may suggest that functional divergence have occurred between the two genes. PnDREB68 showed much higher nucleotide diversity than PnDREB69, suggesting that strong purifying selection has acted on the former gene. Evidence of positive selection on PnDREB69 was inferred by neutrality tests, which revealed a significantly negative value of Tajima's D for its coding region. The different evolutionary forces acting on PnDREB69 and PnDREB68 may explain the expression divergence observed between these genes. Additionally, the two genes also displayed different degrees of population genetic differentiation and linkage disequilibrium. Finally, a set of polymorphisms in PnDREB68 and PnDREB69 may be candidates for association analysis of important traits, such as drought tolerance. Figure S1 Expression of P. trichocarpa DREB69 and DREB68 genes across a range of tissues, organs, and treatments. The patterns of relative transcript accumulation of the two genes were determined by microarray analysis. Red indicates higher levels and blue indicates lower levels of transcript accumulation. Each column represents the average of biological triplicates. ML, mature leaf; YL, young leaf; RT, root; CL, seedlings grown in continuous light; DL, seedlings grown in continuous darkness and then transferred to light for 3 h; CD, seedlings grown in continuous darkness; XY, differentiating xylem; FC, female catkins; and MC, male catkins. (TIF)   Dataset S1 Amino acid sequences of DREB1 and outgroup proteins used in phylogenic analysis. (DOC)