A new genome sequence resource for five invasive fruit flies of agricultural concern: Ceratitis capitata, C. quilicii, C. rosa, Zeugodacus cucurbitae and Bactrocera zonata (Diptera, Tephritidae)

De novo genome assembly

An inbred lab colony of each of the following tephritid species was established in an artificial setting and larvae were collected for subsequent sequencing: Ceratitis capitata, C. quilicii, C. rosa, Zeugodacus cucurbitae and Bactrocera zonata. Inbred specimens of C. quilicii, C. capitata and C. rosa were produced at Citrus Research International in Mbombela and were originally sourced from wild flies collected in Ermelo (-26.516021, 29.996168), Burgershall (-25.112083, 31.087778) and Mbombela (-25.452258, 30.970778), Mpumalanga Province, South Africa respectively in 2020 (C. rosa) and 2021 (C. capitata and C. quilicii). Species identity was confirmed by Marc De Meyer (C. quilicii) and Aruna Manrakhan (C. capitata and C. rosa). Inbred lines for Z. cucurbitae and B. zonata were already present at the facilities of CIRAD, Réunion for more than 150 generations and could thus be used for our purposes. Pupae of all species supplied for sequencing originate from a parent x F1 backcross to increase homozygosity. The sequencing and assembly process can be described by three consecutive steps: generation of PacBio CCS reads and primary assembly with Hifiasm, generation of Hi-C (specifically, Dovetail™ Omni-C™ reads) coupled with secondary assembly using HiRise and lastly, generation of an RNAseq library for ab initio genome annotation. Only the assemblies of C. capitata, C. quilicii and Z. cucurbitae comprised the HiRise scaffolding and annotation steps.

De novo PacBio assembly and filtering

A de novo assembly was constructed using ±38.8 Gb of PacBio CCS reads resulting in a coverage of around 70x of the tephritid genome ( Table 1). The obtained PacBio reads were used as input to Hifiasm v0.15.4-r347 (Cheng et al. 2021) with default parameters. Blast results of the Hifiasm output assembly against the nucleotide BLAST database (https://blast.ncbi.nlm.nih.gov/) were used as input for blobtools v1.1.1 (Laetsch and Blaxter 2017) and scaffolds identified as possible contamination were removed from the assembly. Finally, purge_dups3 v1.2.5 (Guan et al. 2020) was used to purge haplotigs and contig overlaps. The final assembly was checked for its completeness using BUSCO using the diptera_odb10 dataset (Manni et al. 2021).

Chromosome conformation capture and HiRise scaffolding

To construct a Dovetail™ Omni-C™ library, chromatin was fixed in place with formaldehyde in the nucleus and then extracted. Fixed chromatin was digested with DNAse I, chromatin ends were repaired and ligated to a biotinylated bridge adapter followed by proximity ligation of adapter containing ends. After proximity ligation, crosslinks were reversed and the DNA purified. Purified DNA was treated to remove biotin that was not internal to ligated fragments. Sequencing libraries were generated using NEBNext Ultra enzymes and Illumina-compatible adapters. Biotin-containing fragments were isolated using streptavidin beads before PCR enrichment of each library. The library was sequenced on an Illumina HiSeqX platform to produce approximately 30x sequence coverage.

The input de novo assembly and Dovetail™ Omni-C™ library reads (MQ > 50) were used as input data for HiRise, a software pipeline designed specifically for using proximity ligation data to scaffold genome assemblies (Putnam et al. 2016). Dovetail™ Omni-C™ library sequences were aligned to the draft input assembly using bwa (https://github.com/lh3/bwa). The separations of Dovetail™ Omni-C™ read pairs mapped within draft scaffolds were analyzed by HiRise to produce a likelihood model for genomic distance between read pairs, and the model was used to identify and break putative misjoins, to score prospective joins, and make joins above a threshold.

Ab initio genome annotation

Firstly, repeat families in the three tephritid genome assemblies (C. capitata, C. quilicii and Z. cucurbitae) were identified de novo and classified using the software package RepeatModeler2 (Flynn et al. 2020, the original version of RepeatModeler is free and available at https://github.com/Dfam-consortium/RepeatModeler/blob/master/RepeatModeler). The custom repeat library obtained from RepeatModeler2 was used to discover, identify and mask the repeats in the assembly using RepeatMasker (Version 4.1.0, available at https://github.com/rmhubley/RepeatMasker). Secondly, coding sequences from Bactrocera dorsalis, Ceratitis capitata and Drosophila melanogaster available on GenBank were used to train the ab initio model in AUGUSTUS (version 2.5.5) by performing six rounds of optimization. Likewise, the same coding sequences were used to train an independent ab initio gene model using SNAP (Korf 2004). Furthermore, RNAseq reads were mapped onto the genome using the STAR aligner software (Dobin et al. 2013). MAKER (Campbell et al. 2014), SNAP and AUGUSTUS (with intron-exon boundary hints provided from RNAseq) were then used to predict genes in the repeat-masked reference genome. To help guide the prediction process, SwissProt peptide sequences from the UniProt database (https://www.uniprot.org/) were downloaded and used in conjunction with the protein sequences from the aforementioned species to generate peptide evidence in the Maker pipeline (Campbell et al. 2014). Only genes that were predicted by both SNAP and AUGUSTUS were retained in the final gene sets. To help assess the quality of the gene prediction, AED scores were generated for each of the predicted genes as part of the MAKER pipeline. Genes were further characterised for their putative function by performing a BLAST (Ye et al. 2006) search of the peptide sequences against the UniProt database. tRNA were predicted using the software tRNAscan-SE (Lowe & Chan 2016, available at: https://lowelab.ucsc.edu/tRNAscan-SE/).

Phylogenetic tree reconstruction

We inferred orthogroups using OrthoFinder v2.5.5. (Emms & Kelly 2019) for the three fruit fly species with an annotated genome assembly in this study (C. capitata, C. quilicii and Z. cucurbitae). In addition, we downloaded protein sequence data for Drosophila melanogaster Meigen (GCA_000001215.4), Anopheles darlingi Root (GCA_000211455.3), Musca domestica Linnaeus (GCF_030504385.1), Rhagoletis pomonella (Walsh) (GCF_013731165.1) and Bactrocera tryoni (Froggatt) (GCF_016617805.1). Sequences were aligned using Diamond and gene trees were inferred using fasttree. The STAG algorithm combined with the STRIDE rooting methods, implemented in OrthoFinder, was then used to infer a species tree with realistic branch lengths from the full set of gene trees (Emms & Kelly 2017). A time-calibrated tree was constructed by transforming the species tree rendered by Orthofinder into a ultrametric tree and calibrating it based on the split between A. darlingi and the rest of the taxa (240.8 MYA) as inferred from TIMETREE5 (timetree.org).