DNA Methylation in homeologous genes and its correlation with expression for Gossypium petal tissue
Plant Material   Five accessions were used in our study: G. arboreum (2x=2n=26, A2), G. raimondii (2x=2n=26, D5), G. tomentosum (4x=2n=52, AD3), G. hirsutum cv. Maxxa (4x=2n=52, AD1) and a sterile diploid F1-hybrid between A2 and D5 (1x = 1n = 26; F1) (Table 1). The diploid F1-hybrid was created by a hand pollination between reduced gametes of diploids G. arboreum (A2) and G. raimondii (D5), and its somatic cells only contained 13 chromosomes from each extant diploid genome.           Petal tissue was collected from plants growing under controlled greenhouse conditions at the Pohl Conservatory, Iowa State University, USA. Tissue was harvested at time of full petal expansion after dawn but before pollination from three flowers of different plants of each accession (3 biological replicates). Harvested tissue was flash frozen in liquid nitrogen and stored at -80o C until RNA and DNA extraction.     RNA-seq, bisulfite treatment, BS-Seq Libraries  RNA-seq was used to assess gene expression in these petal tissues using Illumina sequencing (Rambani et al., In press). Genome DNA was extracted using Qiagen Plant DNeasy protocols (Qiagen Inc., ). Whole genome bisulfite-sequencing (BS-Seq) libraries were sheared and prepared using the Illumina TruSeq Library Prep kit (Illumina, San Diego, CA) with minor modifications necessitated by the bisulfite treatment. Briefly, five micrograms of genomic DNA from each sample and (a spike-in control of 1-2% unmethylated lambda DNA, (Promega, Madison, WI)) were sheared using the Covaris as per manufacturer instructions. Fragmented DNA was quantified using Picogreen. The cytosine-methylated adapters provided by Illumina were ligated to the blunt ends of fragmented DNA. The ligated DNA was treated with sodium bisulfite using MethylCode™ Bisulfite Conversion Kit (Invitrogen, Grand Island, NY) and a size range (250-350bp) of fragments were selected on a LabChipXT (Caliper Life Sciences, Hopkinton, MA ). Post-size selection samples were enriched by 16 cycles of PCR using 2.5U uracil-insensitive PfuTurboCxHotstart polymerase (Agilent, Santa Clara, CA), 10X PfuTurbo Buffer, dNTPs, and Illumina TruSeq oligo mix (PCR temperatures: 95°C for 2min, then 10 cycles of 98°C for 15 sec, 60°C for 30 sec, 72°C for 4 min, end with 72°C for 10 min). The reaction products were cleaned using Agencourt AMPure XP beads (Beckman Coulter, Inc., Brea, CA). Cleaned PCR products were enriched further using Illumina’s protocol for enrichment of libraries. Libraries were validated on an Agilent Bioanalyzer and quantified using Illumina sequencing primers in a qPCR reaction to sequencing.     DNA Sequencing  The transcriptomes and methylomes of G. arboreum (A2), G. raimondii (D5), G. hirsutum cv. Maxxa, G. tomentosum, and the diploid F1-hybrid were sequenced on an Illumina HiSeq using v.2 chemistry the Huntsman Cancer Institute, University of Utah, Salt Lake City, Utah, USA (Table 1). All reads were trimmed with SICKLE (https://github.com/ najoshi/sickle), using a phred quality threshold of 20.   The unmethylated, trimmed lambda reads were mapped to the lambda phage reference genome sequence using GSNAP (Wu & Nacu, 2010). The number of sequenced cytosines and thymines was tallied at each cytosine position, and the percentage of unconverted cytosines was calculated to estimate the conversion efficiency of the bisulfite treatment.    Mapping and Read Categorization  RNA-seq and bisulfite-treated reads were mapped by GSNAP to the diploid genome reference of G. raimondii [38,39]. To reduce the bias caused by mapping reads from the A- or AT-genomes to the D-genome reference, GSNAP’s SNP-tolerant mapping was used with an index of homoeo-SNPs between various A and D diploid genomes (Page et al., 2013a, In Press). Reads from tetraploids were then categorized into two groups—AT and DT—using PolyCat [40]. The reads from the diploids were also categorized for analyses to obtain unbiased comparisons of gene expression and methylation between the diploid genomes and their respective tetraploid homologous loci.    Methylation Analysis   Mappings of bisulfite-reads were processed by BamBam (REF, website?) to summarize methylation at every cytosine position with at least 5x coverage. This analysis was performed separately for the total of all reads from an accession and separately for the reads categorized as AT and DT. As part of this analysis, the number of converted and unconverted bases mapped to each position was tallied, and a binomial test was performed to determine if the position was methylated (p < 0.05). The results were separated based on the methylation context of the cytosine in question (CG, CHG, or CHH).   The methylation pattern of the “average” gene was estimated by counting the fraction of mapped bases that were cytosines in a sliding window across the gene. The area analyzed included 1 Kbp upstream of the start codon, the exons and introns of the gene normalized to a size of 1 Kbp, and 1 Kbp downstream of the stop codon. The methylated fraction was averaged across all genes to forma a consensus pattern.   Differentially methylated regions were identified by comparing pairs of accessions or genomes in a 10 Kbp sliding window with a step size of 250 bp. A χ2 test of independence was performed for each window (p < 0.05), based on the number of methylated and unmethylated cytosines within the 2 epialleles of the window. Consecutive differentially methylated windows were combined to make full regions. Genes were considered differentially methylated if they overlapped a differentially methylated region.   The phylogenetic relationships between taxa were estimated based on methylation data alone. A distance matrix was constructed by counting the number of methylated vs. unmethylated mismatches between each pair of accessions or genomes. Then a neighbor-joining tree was constructed with PHYLIP on the methylation distances [41]