Identification of age-related CpG sites from longitudinal avian methylomes

Description

Notes

Methods

Zebra finch blood samples were collected during a long-term experiment where birds were housed in outdoor aviaries (320 × 150 × 210 cm) each containing single sex groups of 18–24 adults . Twenty longitudinal samples were collected from ten known-age individuals (six males and four females) sampled twice during their lifetime at an average interval of 1,470 days between the two sampling points. Jackdaw blood samples were collected from individuals of a free-ranging population breeding in nest-boxes south of Groningen, the Netherlands (53.1708°N, 6.6064°E). We analyzed 22 longitudinal blood samples collected during 2007 – 2021 from 11 known age adults (five males and six females) with an average sampling interval of 2,429 days. Samples were taken between June 2008 and December 2014 and stored in -80^oC in EDTA buffer.

We extracted total-cell DNA using innuPREP DNA Mini Kit (Analytik Jena GmbH+Co) from 3μl (nucleated) red blood cells according to the manufacturer's protocol. Whole genome sequencing was performed by The Hospital for Sick Children (Toronto, Canada) where paired-end Illumina next-generation sequencing (150bp) was carried out on either an Illumina HiSeqX™ (12 zebra finch samples) or an Illumina NovaSeq™ sequencer (eight zebra finch samples and 22 jackdaw samples). Libraries were prepared using the Swift Biosciences Inc. Accel NGS Methyl Seq kit (part no. 30024 and 30096) and the DNA was bisulfite converted using the EZ-96 DNA Methylation-Gold kit from (Zymo Research Inc., part no. D5005) as per the manufacturer's protocol and subsequently subjected to whole-genome amplification.  

Sequences were trimmed using Trim Galore!  v. 0.6.10 (38) in paired-end mode. Visual quality controls of the data were carried out before and after trimming using FastQC v. 011.9 (39) and MultiQC v. 11.14 (40). Because the Swift Biosciences Inc. Accel NGS Methyl Seq kit was used for library preparation, the first ~10 bp showed extreme biases in sequence composition and M-bias, so after checking the M-bias plots, the first 10 bps were further trimmed from each sequence.

Alignments were performed using Bismark v. 0.14.433  using the Bowtie 2 v. 2.4.5 alignment algorithm  for both in silico bisulfite conversion of the reference genomes and alignments. For zebra finch, trimmed reads were aligned against the in silico bisulfite converted zebra finch (Taeniopygia guttata) reference genome (GCA_003957565.4) and the average mapping efficiency was 64.5% (SD: 3.22). Because the jackdaw reference genome did not contain an assembled W chromosome, the jackdaw bisulfite sequencing data was aligned to an in silico bisulfite converted Hawaiian crow genome (Corvus hawaiiensis, GCA_020740725.1, the closest relative species with an assembled W chromosome) the average mapping efficiency was 64.7% (SD: 0.98).

The GFT-formatted annotation files for the zebra finch (GCF_003957565.2) and the Hawaiian crow  (GCF_020740725.1) were converted into BED12 format using the University of California, Santa Cruz (UCSC) utilities gtfToGenePred and genePredToBed (available at https://hgdownload.soe.ucsc.edu/downloads.html#utilities_downloads). AR-CpG sites were annotated using the tool annotateWithGeneParts from the R 4.1.2 package genomation 1.4.1. This tool hierarchically classifies the sites into pre-defined functional regions, i.e., promoter, exon, intron, or intergenic, hereon referred to as annotation categories. The predefined functional regions were based on the annotation information present in the BED12 files accessed with the genomation tool readTranscriptFeatures. Annotations were performed as a hierarchical assignment (promoter > exon). Subsequently, a customized R script was employed to integrate the annotation results of AR-CpG sites with their respective annotation category information.