DNA methylation differences between stick insect ecotypes

• 2018Methylation_info_spreadsheet_standard.csv: Table with information regarding samples, locations, climatic information, and bisulfite conversion.
• BSseq_pipeline: The series of scripts below were used in the pipeline to process bifulfite reads
  
  • 1.1_parallel_trimmomatic.sh: Runs Trimmomatic to filter raw bisulfite reads
  • 1.2_sampling24k.sh: Samples a number of reads to reduce batch effects on downstream analyses
  • 1.3_bismark.mapping.to.phage.sh: Runs Bismark to map the bisulfite reads to the lambda phage (GenBank J02459)
  • 1.4_bismark.maping.to.tcristinae.sh: Runs Bismark on unmapped reads to the phage to T. cristinae genome (v1.3c2)
  • 1.5_parallel_bismark_methylation_extractor.sh: Runs Bismark function 'bismark_methylation_extractor' to call methylation into cytosine reports tables
  • 1.6_remove.CT.GA.polymorphisms.pl: This script removes the SNPs listed on the file 'variants.raw.CT.GA.bial.noindel.qs20.cov0.whole.genomes.and.radseq.loci' and remove from the cytosine reports
  • variants.raw.CT.GA.bial.noindel.qs20.cov0.whole.genomes.and.radseq.loci: list of C/T and G/A polymorphisms from new and previously published T. cristinae data
  • Bismark_deduplicating_reads.sl: Runs deduplication
• Annotation: This section contains scripts used to annotate the methylation variation at Timema cristinae species level using the 24 individuals
  • 2.1_get.methylation.status.individual.binomial.pl: It calculates the methylation status based on binomial distributions. Run on the cytosine reports on each individual
  • 2.2_retrieve_annotation_augustus.R: Determines the annotation of each methylation position based on the annotation file from Villoutreix et al. 2020.
  • 2.3_retrieve.genes.id.R: Gets the gene id based on the annotation file from Villoutreix et al. 2020, after running 2_retrieve_annotation_augustus.R
  • 2.4_retrieve_annotation_repeatmasker_1.3c2.R: Gets the repeats annotation based on the repeatable elements annotation from Villoutreix et al. 2020.
  • 2.5_retrieve.exon.intron.oreder.on.meth.table.R: Gets the exons and introns in the order (e.g. CDS1, CDS2, etc.)
  • 2.6_level.methylation.exons.introns.R: Estimates methylation levels on different exons and introns.
  • 2.7_enrichment.genomic.features.R: Estimates the ernichment of methylation levels on different genomic features
  • 2.8_GO.enrichment.genes.R: Estimates the enrichment of GO terms in genes that are hypo or hyper methylated.
  • first.batch.compiled.noCT.GA.low5.high60.binomial.annot.12samples.txt: Compilation of the annotation tables. Here, only loci covered by a minimum of 5 reads and maximum of 60 were retained. We also selected loci present at at least 12 samples. This table shows the methylation status at each loci (based on the binomial distribution)
  • first.batch.compiled.12samples.no.intergenic.noCT.GA.low5.high60.annot.txt: Same as above, but here the mean methylation levels were calculated. Intergenic regions were removed to ease the analyses
  • first.low5.high60.noCT.GA.methylation.status.mRNA.GO: Table with gene ids. Formatted to input at script 8_GO.enrichment.genes.R.
• RNA-seq: This section contains scripts used to process RNA-seq data
  • 3.1_cutadapt_filtering.sh: Filters adapters from the data
  • 3.2_trimmomatic_filtering.sh: Runs Trimmomatic
  • 3.3_mapping_array_STAR_relaxed_pe.sh: Runs STAR to map RNA-seq data to T. cristinae reference genome (v1.3c2)
  • 3.4_featureCounts_Tcristinae_genes.sl: Peforms featureCounts function
  • 3.5_plot_expression_methylation.R: Estimates relationship between methylation levels and expression data
• Genome_wide_comparison: This section contains scripts and inputs from the genome-wide analyses
  • 4.1.methylation_genetic_mantel_bayesian.R: Runs mantel tests and Bayesian regressions
  • MethylRaw_CpG_first_run_low2_high60_percentage.no.CT.GA: Table compiling methylation cytosine reports among all 24 samples. Loci with minimum 2 reads and maximum 60 (above 99th quantile) were removed. The methylation levels were calculated as methylated cytosines at a certain locus over the sum of all reads covering it.
  • jtdistance.matrix.txt: Genetic distances between the 24 individuals based on the RAD-seq data
• MACAU: This section contains scripts, inputs and outputs related to MACAU analyses
  • 5.1_methylKit.tiles.R: Runs methylKit and summarizes the data into 1kbp tiles
  • 5.2_tiles.filtering.before.macau.R: Removes tiles that are hypo and hyper metylated following Lea et al. (2016)
  • 5.3_MethylRaw_formatting_first.sh: Formats the table generated by methylKit into MACAU inputs
  • 5.4_macau.first.sh: Runs MACAU
  • 5.5_ibd_methylation_by_host_diff_cutoffs.R: Runs bayesian regressions on the outputs from MACAU
  • 5.6_GO.term.DMR.cutoff.R: Estimates GO enrichments on DMRs from different p-value cutoffs
  • MethylRaw_methylC_first_methyl_tiles_low10: methylated cytosine counts for MACAU
  • MethylRaw_coverage_first_methyl_tiles_low10: coverage input for MACAU
  • 2018Methylation_covariates_batch_first.txt: Covariates input at MACAU, namely: PC1 and PC2 from climate, bisulfite conversion (calculated based on the lambda phage), and sequencing batch
  • 2018Methylation_covariates_batch_first.txt: Host-plant predictor for MACAU
  • 2018Methylation_first_run_relatedness.cXX.txt: Kinship matrix calculated based on RAD-seq. Performed using gemma with default parameters (Zhou et al. 2013).
  • macau.first.output.CpG.tiles.low10.var.adjusted.txt: MACAU output. We disregarded the unassembled scaffolds (i.e. lgNA).
  • tiles.macau.low10.genes.id.order.txt: MACAU table with gene.id formatted to estimate GO enrichment. The annotation was performed similarly to the scripts above

Funding provided by: Royal Society
Crossref Funder Registry ID: http://dx.doi.org/10.13039/501100000288
Award Number: RG140369

Funding provided by: European Research Council
Crossref Funder Registry ID: http://dx.doi.org/10.13039/501100000781
Award Number: R/129639