RAD tag sequencing as a source of SNP markers in Cynara cardunculus L
© Scaglione et al; licensee BioMed Central Ltd. 2012
Received: 6 August 2011
Accepted: 3 January 2012
Published: 3 January 2012
Skip to main content
© Scaglione et al; licensee BioMed Central Ltd. 2012
Received: 6 August 2011
Accepted: 3 January 2012
Published: 3 January 2012
The globe artichoke (Cynara cardunculus L. var. scolymus) genome is relatively poorly explored, especially compared to those of the other major Asteraceae crops sunflower and lettuce. No SNP markers are in the public domain. We have combined the recently developed restriction-site associated DNA (RAD) approach with the Illumina DNA sequencing platform to effect the rapid and mass discovery of SNP markers for C. cardunculus.
RAD tags were sequenced from the genomic DNA of three C. cardunculus mapping population parents, generating 9.7 million reads, corresponding to ~1 Gbp of sequence. An assembly based on paired ends produced ~6.0 Mbp of genomic sequence, separated into ~19,000 contigs (mean length 312 bp), of which ~21% were fragments of putative coding sequence. The shared sequences allowed for the discovery of ~34,000 SNPs and nearly 800 indels, equivalent to a SNP frequency of 5.6 per 1,000 nt, and an indel frequency of 0.2 per 1,000 nt. A sample of heterozygous SNP loci was mapped by CAPS assays and this exercise provided validation of our mining criteria. The repetitive fraction of the genome had a high representation of retrotransposon sequence, followed by simple repeats, AT-low complexity regions and mobile DNA elements. The genomic k-mers distribution and CpG rate of C. cardunculus, compared with data derived from three whole genome-sequenced dicots species, provided a further evidence of the random representation of the C. cardunculus genome generated by RAD sampling.
The RAD tag sequencing approach is a cost-effective and rapid method to develop SNP markers in a highly heterozygous species. Our approach permitted to generate a large and robust SNP datasets by the adoption of optimized filtering criteria.
Cynara cardunculus (2n = 2x = 34, haploid genome size ~1.08 Gbp ) an allogamous, highly heterozygous Asteraceae species, includes three taxa: the globe artichoke (var. scolymus), the cultivated cardoon (var. altilis) and their common progenitor the wild cardoon (var. sylvestris) . Globe artichoke contributes significantly to the Mediterranean agricultural economy, and is also cultivated in South America, North Africa, China and USA. Over the past 30 years, a body of evidence has grown that plant-based foods can be effective for the alleviation of several chronic diseases, and globe artichoke in particular has been shown to produce a number of nutraceutically and pharmaceutically active compounds. Extracts from both globe artichoke and cultivated cardoon have exhibited hepatoprotective, anticarcinogenic, antioxidative and antibacterial qualities, and even an inhibition of cholesterol biosynthesis and LDL oxidation [3–6]. Finally, there is increasing interest in developing the species as an energy and oilseed crop [7–10].
Since the first linkage map produced for globe artichoke , a number of other segregating populations have been exploited for genetic mapping, including one generated from a hybrid between a globe artichoke and a cultivated cardoon genotype  and, more recently, one obtained by crossing globe artichoke with wild cardoon . The recent development of a set of gene-based microsatellites  has aided the construction of consensus genetic maps [13, 15, 16]. However, these maps remains insufficiently densely populated for trait mapping and marker assisted selection. Current high throughput sequencing technology, which produces DNA sequence at a rate several orders of magnitude faster than conventional methods, is effective as a platform for SNP (single nucleotide polymorphism) discovery. A particularly efficient protocol, termed "restriction-site associated DNA" (RAD) , in combination with the Illumina Genome Analyzer sequencing device , discovers SNPs by sequencing a large set of restriction fragments [19–21]. Here we report the generation of genomic RAD tags from the three C. cardunculus accessions used as the parents for two of our mapping populations. The RAD tags were used to derive SNP markers some of which were then validated by a Cleaved Amplified Polymorphic Sequence (CAPS) assay. The identified SNPs could be useful to produce denser C. cardunculus genetic maps via high-throughput genotyping technologies. The RAD sequence has also been informative for characterizing the repetitive DNA component of the C. cardunculus genome, in particular allowing some inferences to be made regarding the contribution of DNA methylation in inhibiting its expansion.
Altogether, our data suggest that the RAD procedure, despite its use of GC-rich recognition sites, has produced a random representation of the C. cardunculus genome, and shows that it represents a reliable means of assessing genome complexity.
SNP mining results.
Total SNPs mining (CcRAD1)
"Fully informative" RAD loci (CcRAD2)
Putative testcross markers (CcRAD2)
"Romanesco C3" testcross over "Altilis 41"
"Altilis 41" testcross over "Romanesco C3"
"Romanesco C3" testcross over "Creta 4"
"Creta 4" testcross over "Romanesco C3"
Common intercross markers (CcRAD2)
CAPs markers conversion.
Product size (bp)
Restriction site (bp)
"RomanescoC3" restriction produts
"Altilis 41" restriction products
New LG Alt_22
In crop species where the number of markers available to date is limiting, the use of high throughput sequencing to generate large numbers of genetically informative assays can make a valuable and rapid contribution to linkage mapping, and its major downstream application, marker-assisted selection. RAD tag sequencing based on the Illumina platform has proven to be a highly reliable and cost-effective means of SNP discovery. We were able to identify thousands of putative SNP markers in this way, and the majority of a random sample of 24 was fully validated through conversion to CAPS assays and linkage analysis. Furthermore, the reduction in template complexity generated by the RAD approach greatly facilitates its implementation in mapping-by-sequencing approaches.
A large proportion of the methylation present in DNA occurs in the form of CpG dinucleotides, and there is little evidence for negative selection against these in the many genomes which have been analysed to date [32, 37]. Acquiring genome-wide sequence has given a glimpse of the genome complexity present in C. cardunculus. Even though the RAD tags represent only a sample of the genome as a whole, it was clear that there exists a relationship between the frequency of CpG dinucleotides and the level of sequence repetitiveness, consistent with the known role played by methylation in controlling genome expansion due to transposable element activity [30, 31].
Genomic DNA was extracted from the leaf of the three C. cardunculus accessions, following the protocol described by Lanteri et al. . The three accessions have been used as parents of two F1 populations, made by crossing globe artichoke variety "Romanesco C3" as female with either the cultivated cardoon variety "Altilis 41" or the wild cardoon accession "Creta 4" as male . "Romanesco C3" is a late-maturing variety, which forms large purple-green capitula, each bearing violet coloured florets; "Altilis 41" was selected at the University of Catania  on the basis of its biomass yield potential; its foliage is grey and its florets white. "Creta 4" was collected from a wild population in Crete; it produces a large number of capitula, forms green-violet bracts and violet florets. Each DNA sample was processed into a separate RAD libraries as reported by Baird et al. . Briefly, 300 ng DNA were digested with 20 U of PstI (New England Biolabs, NEB) for 60 min at 37°C in a 50 μl reaction, after which the reactions were heat inactivated by holding at 65°C for 20 min. A 2.5 μL aliquot of 100 nM P1 adaptor (a modified Illumina adapter)  was added to each sample along with 1 μL 10 mM ATP (Promega), 1 μL 10x NEB Buffer4, 1,000 U T4 DNA ligase (Enzymatics, Inc) and 5 μL H2O, and the reaction was incubated at room temperature for 20 min, ending with a heat inactivation step (65°C/20 min). The reactions were then pooled and sheared to an average length of 500 bp using a Bioruptor (Diagenode). The sheared DNA was separated by electrophoresis through a 1.5% agarose gel, and fragments in the 300-800 bp range were isolated using a MinElute Gel Extraction kit (Qiagen). The End-Repair mix (Enzymatics, Inc.) was used to blunten the dsDNA ends, and the samples were re-purified using a MinElute column (Qiagen), following which 15 U Exo-Klenow (Enzymatics, Inc.) were added and the sample incubated at 37°C to generate 3'-adenine overhangs. After subsequent purification, 1 μL 10 μM P2 adapter (a second modified Illumina adapter)  was ligated and the sample purified as above. The concentration of DNA in the eluate was quantified using a Qubit fluorimeter, and a 20 ng aliquot was used for a 100 μL PCR comprising 20 μL Phusion Master Mix (NEB), 5 μL 10 μM P2 and H2O. The 18 cycle PCR amplification regime followed the recommendation of the manufacturer (NEB). After this PCR, the samples were separated by electrophoresis once again through a 1.5% agarose gel, and fragments in the 300-700 bp range were excised from the gel and diluted to 3 ng/μL. The material was analysed on an Illumina Genome Analyzer IIx following the paired ends (2x 54 bp) genomic DNA sequencing protocol suggested by the manufacturer.
The sequences were sorted according to their multiplex identifier tag. A RAD LongRead® contig assembly was generated by a set of algorithms developed at Floragenex Inc. Sequences having more than 5 bases with poor Illumina quality scores (Phred10 or lower) were discarded. Paired reads were collapsed into sequence "clusters" on the basis of single ends (SE) sharing 100% sequence identity. To maximize assembly efficiency, a minimum of 25x and maximum 400x sequence coverage at RAD SE reads were imposed. The variable paired end sequences for each common SE were extracted using the filtered sequence set and compiled for the LongRead® contig construction, using a modified version of the Velvet sequence assembler (v. 1.0.04)  and testing several k-mers in graph construction for each RAD contig. After analysis of the first-pass assembly from each template, "Creta 4" was selected as the reference sequence set. Additional filters were then applied to remove short contigs (< 100 bp in length), low paired end coverage (< 4.0x) or ambiguous contigs (containing N's homopolymers). If more than a single contig (NODE1) was assembled for a given RAD locus, alternative ones were retained in the dataset and labelled accordingly (NODE2, NODE3).
RAD contigs were annotated using Blast2GO software , and were submitted to the NCBI nr protein database where an E-value of 10e-3 or lower were retrieved (20 best hits recorded). Gene names and GIs (gene identifiers) were assigned according to NCBI guidelines, and PIR (Protein Information Resource) identifiers in reference to UniProt, SwissProt, TrEMBL, RefSeq, GenPept and PDB. The annotation was obtained by applying the formula embedded in Blast2GO , setting a threshold score of 55. In the Blast2Go pipeline, GO terms are "transferred" to query sequences only whether a score threshold is reached. This score is calculated basing on both sequence similarity and presence of children node in the directed acyclic graph (DAG). Therefore, in this scenario the first e-value cut-off is used only for the purpose of "collecting" GO-terms, while other more stringent criteria are ruling whether transfer these terms to our sequences. Enzyme codes were retrieved from GO tables and mapped onto KEGG pathways. Transposable elements were detected using RepeatMasker v3.2.9 software http://www.repeatmasker.org, based on the RMBlast algorithm. Default parameters (except for -s flag) were used to search against Viridiplantae repeats.
where CpG represents the observed frequency of CpG dinucleotides and p(C) and p(G) the respective frequencies of each single nucleotide.
MAQ software (v. 0.5.0)  was used to align the paired end reads in the "Creta 4" reference contig set. The alignment threshold was set to a maximum of three nucleotides mismatch between Illumina reads and the reference. Gaps in the alignment of up to 2 nt allowed. Two levels of stringency were applied. In the first (CcRAD1), a comprehensive list of putative SNPs and 1-2 bp indels was populated with a minimum coverage of 6x as threshold prior to uploading to a Microsoft Access relational database; and for the second (CcRAD2), "fully informative" SNPs were defined when a minimum of 1-read allele calling was achieved for each of the three samples. In the latter set, heterozygous SNPs were assessed where the within sample allele frequency ranged from 0.25 to 0.75, together with a minimum coverage of 4x and allele calling for two reads. Sites were assigned as homozygous when the minor allele frequency fell below 0.10.
Candidate SNP markers were categorized as testcross in pair-wise comparisons of genotypes, whether a heterozygous imputation was present for one parent only (testcross) and a homozygous site was predicted for the other. Common intercross markers were defined for loci showing heterozygous states across all the three samples.
A subset of heterozygous SNPs was selected from the "Altilis 41" sequence, and a search carried out for BamHI, EcoRI, EcoRV, NdeI, XbaI, BccI, FokI, XmnI and DraI (6 bp cutters), or TaqI and MseI (4 bp cutters) recognition sites using SNP2CAPS script (v. 0.6) . A predicted fragment size difference of at least 20 bp was imposed to allow detection on standard agarose gels. Locus-specific primers were designed from the BatchPrimer3 web interface , using default parameters but for product size (100-400 bp) and annealing sites (within a 50 bp window at either end of the RAD contig). The resulting assays were applied to a set of 94 F1 segregants from the cross "Romanesco C3" × "Altilis 41" . PCRs were carried out in a 20 μl volume containing 12.5 ng genomic DNA, 1x GoTaq Buffer (Promega), 1.5 mM MgCl2, 0.2 mM dNTPs, 1 U GoTaq (Promega) and 0.5 μM of each primer. The cycling regime was 95°C/5 min, followed by 35 cycles of 95°C/30 s, 55°C/30 s, 72°C/45 s and a final incubation of 72°C/5 min. Amplification was checked by electrophoresis through a 1.5% agarose gel and quantified using a Beckman Coulter spectrophotometer. Restriction reactions (20 μl) comprised 800 ng amplified DNA, 0.3 U restriction enzyme (New England Biolabs), reaction buffer and BSA according to the manufacturers' specifications, incubated for 4 h at 37°C (except for TaqI, where the incubation temperature was 65°C), after which the reactions were heat inactivated (80°C/10 min). The resulting products were electrophoresed through 2% agarose gels.
The CAPS derived genotypic data were incorporated into a pre-existing data set of 273 molecular loci, mainly AFLP and EST-SSRs, already used to generate the cultivated cardoon genetic map [11, 14, 15] including five SNP from genes underlying caffeoylquinic acids synthesis reported by Comino et al.  and Menin et al. ; all maps data are available on request by the authors. Goodness-of-fit between observed and expected segregation ratios was tested by χ2 and only markers fitting or deviating only marginally from expectation (χ2 α = 1 < χ2 ≤ χ2 α = 0.01) were included for mapping. Linkage groups (LGs) were established by JoinMap v4.0 software , on the basis of a LOD threshold of 6.0, using as parameter settings Rec = 0.40, LOD = 1.0, Jump = 5. Map distances were converted to centiMorgans (cM) using the Kosambi mapping function. LGs were drawn and aligned using MapChart v2.1 .
This research was supported by grants from: (i) the National Science Foundation Plant Genome Research Program (No. 0421630), (ii) the Georgia Research Alliance, (iii) the University of Georgia Research Foundation, and (iv) by MIPAAF (Ministero delle Politiche Agricole, Alimentari e Forestali - Italy) through the CYNERGIA ("Costituzione e valutazione dell'adattabilita' di genotipi di Cynara cardunculus per la produzione di biomassa e biodiesel in ambiente mediterraneo") project and CARVARVI ("Valorizzazione di germoplasma di carciofo attraverso la costituzione varietale ed il risanamento da virus") project.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.