Chromosome-level genome assembly of the synanthropic fly Chrysomya megacephala: insights into oviposition location

Abstract

The oriental latrine fly, Chrysomya megacephala (Diptera: Calliphoridae), is a medically important synanthropic blow fly species characterized by its necrophagy and coprophagy, often observed near carrion and animal feces. Notably, C. megacephala always arrives at carcass earlier than other species. To elucidate the underlying mechanisms behind the host choice in C. megacephala, we present the chromosome-scale genome assembly for this species. The genome size is 816.79 Mb, with a contig N50 of 1.60 Mb. The Hi-C data were anchored to six chromosomes, accounting for 99.93% of the draft assembled genome. Comparative genomic analysis revealed significant expansions in pathways of ligand-gated ion channel activity, passive transmembrane transporter activity, and protein methyltransferase activity, which may be closely associated with host localization and oviposition. After identifying 69 odor-binding proteins (OBPs) in the assembled genome, phylogenetic analysis showed that DmelOBP99b and CmegOBP99b exhibited high homology. Transcriptome analysis demonstrated that the relative expression of CmegOBP99b was consistently the highest during the metamorphosis, and RT-qPCR further confirmed the similar results. Additionally, CmegOBP99b exhibited a strong binding affinity to DMDS (dimethyl disulfide) as determined by molecular docking. To determine the protein expression level of CmegOBP99b in various body parts, we prepared recombinant CmegOBP99b protein and anti-CmegOBP99b polyclonal antibodies. Western blot analysis showed that CmegOBP99b was significantly expressed in the female’s head compared to other parts, which is consistent with RT-qPCR results. Therefore, CmegOBP99b may be the primary odor-binding protein responsible for olfactory recognition and the behavioral coordination of C. megacephala. This study not only provides valuable insights into the molecular mechanisms of oviposition localization in C. megacephala but also facilitates further research into the genetic diversity and phylogeny of the Calliphoridae family.

Introduction

The oriental latrine fly, Chrysomya megacephala (Fabricius, 1794) (Diptera: Calliphoridae), is a medically important synanthropic blow fly species that originated in Australasia and the Pacific and is now distributed worldwide [1]. Adult C. megacephala are attracted to a wide range of substrates, including human food, livestock and human feces, and carrion. Due to its ability to utilize both human feces and carrion as breeding sites, the species’ endophilic behavior contributes to the spread of fecal pathogens, posing a notable public health concern [2]. In addition to reported cases of myiasis caused by C. megacephala [3, 4], strains of Escherichia coli isolated from wild C. megacephala in China has shown 100% resistance to multiple common antibiotics, further emphasizing its potential health threat [5]. Moreover, C. megacephala is considered one of the forensically important flies as it can be used to estimate the postmortem interval (PMI) based on the development of larvae that colonize decomposed corpses [6]. Beyond its medical and forensic relevance, C. megacephala plays a role in modern agriculture as an important pollinator for orchards and vegetables [7]. Additionally, C. megacephala exhibits notable biochemical activities. Its larvae produce native secretions containing antimicrobial peptides that inhibit bacterial growth, such as Staphylococcus aureus [8]. The larvae are also emerging as new sustainable resources for the production of insect proteins, lipids, chitosan, and biofuels [9, 10]. What's more, the larvae of C. megacephala are also effective in biodegrading organic wastes [11, 12].

At present, the morphological characteristics of C. megacephala during non-adult stage are analyzed using both physical and chemical methods. For example, scanning electron microscopy (SEM) provides detailed morphological characters to identify the larvae and puparia [13]. The development age can be estimated using the weathering patterns of hydrocarbons in empty puparia [14]. The cuticular chemical profiles of C. megacephala at different developmental stages are assessed by mid-infrared photoacoustic spectroscopy so as to identify the age of individuals and distinguish different populations [15]. Furthermore, the effects of biological and abiotic factors on the growth and development of C. megacephala have also been extensively studied, such as different constant temperatures [16], flunitrazepam [17], ketamine [18], aluminium phosphide [19, 20], malathion [21], and dietary fat levels [22]. Additionally, a previous study indicated that high concentrations of latex from Parahancornia amapa could alter the post-embryonic development of C. megacephala, suggesting a potential avenue for biopesticide development [23].

Although blow flies are known to be diurnal, the nocturnal oviposition of blow flies is controversial. For example, Chen and Shiao investigated the effects of chronobiology on the reproductive behavior of C. megacephala, and proved that oviposition can occur at night [24]. Additionally, intense competition exists between the larvae of C. megacephala and Chrysomya rufifacies, with significant differences in their oviposition preferences [25]. While C. megacephala females are primarily active during daylight, studies demonstrate that they are more likely to lay eggs in darkness as ambient temperature rises [26]. Females apparently prefer walking to nocturnal oviposition sites [27]. Besides, the occurrence of precocious eggs in C. megacephala has been noted as a topic of concern [6].

C. megacephala is typically one of the earliest colonizers of corpses and the dominant fly species found in maggot masses [28]. Necrophagous flies rely primarily on volatile organic compounds (VOCs) to locate decomposing carcasses [29]. The selection of suitable oviposition sites is closely linked to the ability of gravid females to detect and recognize specific VOCs [30]. Among the common VOCs emitted from carrion and feces, dimethyl disulfide (DMDS) and dimethyl trisulfide (DMTS) are key signals that necrophagous insects use to locate breeding sites [31]. Olfactory proteins play a critical role in the olfactory system, which is commonly used by insects to locate hosts and oviposition sites [32]. The main peripheral olfactory proteins involved in the detection of odorants in insects are odorant-binding proteins (OBPs), chemosensory proteins (CSPs), odorant receptors (ORs), ionotropic receptors (IRs), and sensory neuron membrane proteins (SNMPs). OBPs and CSPs are involved in the initial step of olfaction. The odorant molecules in the environment enter the lymph of antennal sensilla through the pore tubules and interact with the OBPs or combine with CSPs to form a complex [33]. Therefore, olfactory proteins in C. megacephala are essential for colonizing and laying eggs on carcasses. Due to limited availability of genomes, Wang et al. established the developmental expression profiles from egg to adult, identifying two OBPs (Cmeg32081-c4 and Cmeg33593_c0) that may play crucial roles in linking the olfactory system to broader biological processes [34].

Although C. megacephala is the most common human-associated fly with medical and forensic importance, the limited availability of genomes has severely hindered the study of its molecular mechanisms. Currently, one or more assemblies of 12 species of the family Calliphoridae are referenced in NCBI [35,36,37,38,39], five of which belong to the same subfamily Chrysomyinae, as C. megacephala. We hereby report a chromosome-level de novo genome assembly of C. megacephala and perform comparative analysis with other dipteran flies with genome available to enrich understanding of the molecular mechanisms involved in the process of oviposition localization in C. megacephala. Moreover, we comprehensively annotated the OBPs family using the high-quality genome assembly, and then a phylogenetic tree was constructed. We identified that CmegOBP99b might be involved in the response to DMDS based on molecular docking analysis, which revealed a strong binding affinity between CmegOBP99b and DMDS. Following this, the expression profiles of CmegOBP99b were further validated at both the mRNA and protein levels. This study provides a valuable resource for understanding the molecular mechanism underlying oviposition localization in C. megacephala, which sheds insight into the biological habits and genetic diversity of blow flies.

Materials and methods

Cultivation of C. megacephala

Adult specimens of C. megacephala were trapped using pork liver bait in Changsha, Hunan Province, China [112º59'E, 28º12'N] (Fig. S1). Species identification was performed by an expert dipterological taxonomist (Lushi Chen) based on descriptive morphological characteristics using a dissecting microscope [40], and then raised in an artificial climate chamber (GPL-250 A, Tianjin Laboratory Instrument Equipment Co. Ltd., Tianjin, China) at 25 ± 1 °C, relative humidity of 70 ± 5%, and a photoperiod regime of 12:12 h light/darkness. Pork liver was used as a medium for larvipositing and larval rearing. In order to reduce genetic variability, mating pairs of adult C. megacephala trapped in the wild were subjected to six generations of highly inbred breeding. The sampled larvae were observed to determine the instar based on the number of clefts in the posterior spiracle. The larval stage includes the first, second, and third instar. The 3rd-instar that just jumped into the sand are defined as being in the wandering stage. The wandering larvae eventually reach metamorphosis, which is referred to as the pupal period. All samples were preserved at –80 °C immediately after processing in liquid nitrogen.

Generation of short read data and genome survey

Genomic DNA was extracted from a single adult female of C. megacephala using the SDS method (Table S1) [41]. Genomic DNA was randomly fragmented using Covaris. The fragmented DNA was selected by Agencourt AMPure XP-Medium kit to an average size of 300–500 bp. The qualified libraries were sequenced on MGISEQ2000 platform. Quality control of the raw data was performed with fastp v0.21.0 [42]. To understand the genomic characteristics of C. megacephala, the 17-mer frequency distribution analysis was performed using jellyfish program before genome assembly to estimate the genome size and heterozygosity with GenomeScope v2.0 [43, 44].

Additionally, we further estimated the genome size of C. megacephala using flow cytometry (BD FacsCalibur, USA) [45]. In order to measure the DNA content of cells, chicken (Gallus gallus domesticus) erythrocytes (1.25 pg) are commonly used as an internal reference. Three adult female C. megacephala specimens were treated and incubated with the DNA fluorochrome propidium iodide (PI). The relative fluorescence of the stained nuclei was then calculated [46]. The cytometer was equipped with an argon ion laser operating at 488 nm. The PI fluorescence was collected using a 620 nm fluorescence-2 (FL2) filter. The results were evaluated based on the average of coefficient of variation (CV) for G1 peaks. The results with a coefficient of variation (CV) of less than 5% were considered as reliable. Histograms were analyzed by Modifit 3.0 software [47].

Genome sequencing and de novo assembly

Genomic DNA was extracted from a single adult female of C. megacephala using the QIAGEN® Genomic kit (QIAGEN, Germany) following the standard procedure (Table S2). SMRTbell target size libraries were constructed using 15 kb preparation solutions according to PacBio’s standard protocol (Pacific Biosciences, CA, USA). Sequencing was performed on a PacBio Sequel II instrument at Haorui Genomics (Xian, China). Raw data was processed through the SMRT Analysis software v10.1.0. Subsequently, the HiFi reads were produced by the CCS subprogram with default parameters. The genome was de novo assembled using an phased string graph with hifiasm v0.15.5 [48]. Completeness of the genome assembly was assessed using BUSCO v5.2.2 [49]. To assess the accuracy of this assembly, the paired end reads were mapped to the assembled genome using BWA v0.7.17 (Burrows-Wheeler Aligner) [50]. We then evaluated the mapping rate and genome coverage of the sequencing reads using SAMtools v1.13 [51]. The base accuracy of the assembly was calculated with bcftools v1.13 [52]. In addition, we applied the GC depth analysis to evaluate whether potential contamination remained during sequencing and the coverage of the assembly, which was implemented in minimap v2.17 [53].

The Hi-C library preparation was constructed following a previous protocol [54]. Briefly, fresh head and thorax were fixed with 1% formaldehyde to crosslink DNA–DNA interactions [55]. Cross-linked DNA was digested with the restriction enzyme (Dpn II) [56]. Hi- C libraries were constructed by NEBNext® Ultra™ II DNA library Prep Kit for Illumina (NEB) following manufacturers’ instructions. The libraries were then sequenced for 150 bp paired-end reads on the Illumina HiSeq X Ten platform (San Diego, CA, United States). To construct a chromosomal-level assembly of the genome, the Hi-C raw data was first trimmed by fastp v0.20.0 [42]. After performing quality control, the clean paired-end reads were aligned to the draft assembled genome using the Juicer pipeline v2.3.2 [57]. Hi-C reads were used to anchor the assembled contigs to scaffolds by Bowtie2 in HiC-Pro v2.9.0 [58]. According to the orders and orientations based on the alignment results, these scaffolds were clustered into chromosomes by LACHESIS [59]. The assembled chromosome-level genome was split into bins of 500 kb in order to construct an interaction heatmap for validation and manual correction.

RNA-seq and data analysis

To assist in genome annotation, the transcriptome of C. megacephala in different developmental stages was performed. A total of 15 samples were collected, including the wandering larvae (n = 3, Sample B1), the prepupal stage (n = 3, Sample B2), the middle pupal stage (n = 3, Sample B3), the late pupal stage (n = 3, Sample B4), and the newly emerged adults (n = 3, Sample B5). Three biological replicates were performed for each group. Subsequently, total RNA was extracted using the mirVana miRNA Isolation Kit (Ambion) following the manufacturer’s protocol. DNA contaminants were removed using DNase I (Promega, Madison, WI, USA). RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The libraries were constructed using the TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA) in accordance with the manufacturer’s instructions. Then, these libraries were sequenced on the Illumina sequencing platform (Novaseq 6000), and 150 bp paired-end reads were generated at Major biotech Co Ltd (Shanghai, China). Raw reads were assessed for quality control using Trimmomatic [60], and then mapped to the reference genome using hisat 2.2 [61]. Only aligned reads were further analyzed and annotated using the reference genome. Subsequently, the read counts of each gene were obtained by htseq count [62], and the fragments per kb per million reads (FPKM) value of each gene was calculated using cufflinks [63, 64]. Furthermore, several databases were used to annotate gene functions, including the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Swiss-Prot/TrEMBL/InterPro databases.

Gene prediction and functional annotation

We employed de novo and homology-based approaches to identify the transposable elements (TEs) and tandem repeats in repetitive sequences. RepeatMasker v4.0.772 and RepeatProteinMask were used to align the transposable elements against the known RepBase (v21.0173) [65,66,67]. In the de novo method, a custom repeat library was first constructed using the RepeatModeler v2.0.3 [65]. RepeatMasker v4.0.772 was then used to identify repetitive sequences against these libraries [65]. In addition, tandem repeats were identified using Tandem Repeats Finder v4.0977 [68]. These repetitive sequences were finally combined and condensed to create a non-redundant repeat annotation of the genome. Additionally, non-coding RNAs (ncRNAs) were annotated using BLAST (E-value ≤ 1e − 5) against the Rfam database [69, 70], which includes microRNAs (miRNAs), ribosomal RNAs (rRNAs), small nuclear RNAs (snRNAs), and transfer RNAs (tRNAs). RNAmmer v1.2 was used to annotate the rRNAs and their subunits [71]. tRNAs were identified using tRNAscan-SE v1.3.1 with default parameters [72].

We combined homology searches, de novo prediction, and transcriptome data-based approaches to predict the protein-coding gene structures of C. megacephala. In the homology-based method, protein sequences from five dipteran insects (Chrysomya rufifacies, Phormia regina, Sarcophaga peregrina, Drosophila melanogaster, Lucilia cuprina, and Musca domestica) were used as queries to search in the assembled genome using the GeneWise v2.4.1 [73]. The de novo predictions were performed from the homology-based predictions to train model parameters using the Augustus v3.0 [74], SNAP [75], GlimmerHMM [76], and GeneID v1.4.4 [77]. Meanwhile, transcriptome data was utilized to align against the genome assembly through PASA and TopHat, respectively [78, 79]. Subsequently, we integrated all predicted genes to generate a consensus gene set via EVidenceModeler v1.1.1 [78]. Finally, all gene sets obtained from gene annotation were aligned against the PFAM, NR, NT, Uniprot, KEGG, and GO databases using BLAST v2.2.3179 with an e-value ≤ 1e − 5 [69]. The conserved domains of the predicted proteins were then identified using Batch CD-search [80].

Comparative genomics and phylogenetic analyses

Comparative genomic and evolutionary analyses were performed among C. megacephala and 10 related species, including the species of C. rufifacies, S. peregrina, M. domestica, Lucilia sericata, L. cuprina, P. regina, Calliphora vicina, Sarcophaga bullata, Aldrichina grahami and D. melanogaster. These protein sequences were downloaded from the NCBI database. The extracted protein sequences were then aligned pairwise to identify conserved orthologs using Blastp v2.12.0 with an e-value of 1e-5 [81]. One-to-one orthologs were further identified by OrthoFinder v2.5.4 [82]. Based on the identified orthologous gene sets, a phylogenetic analysis was performed using the shared single-copy genes. Briefly, the orthologous sequences were multiple aligned using MUSCLE [83]. The phylogenetic tree was constructed with 1000 bootstrap replicates using a maximum likelihood method implemented in RAxML v8.2 under the GTRGAMMA substitution model [84, 85]. Afterwards, divergence times between lineages were estimated under the mean substitution rates along each branch using MCMCTREE implemented in PAML v 4.9e [86]. According to the results mentioned above, expansion and contraction of orthologous gene families were further estimated by CAFÉ v4.2 [87].

Subsequently, we calculated the average ratio of nonsynonymous substitution rate (Ka) and synonymous substitution rate (Ks) of protein coding genes. We then conducted a likelihood ratio test (P < 0.05) using Codeml with the branch-site model implemented in the PAML v4.8 package to identify positively selected genes [86]. Finally, functional annotation was performed on specific genes that should significant expansion or contraction. These genes were then selected for KEGG pathway and GO enrichment analysis, along with the positively selected genes. In addition, we searched for species-pairwise synteny blocks using LAST v1293 [88]. We then analyzed and visualized the chromosome synteny between C. megacephala and D. melanogaster based on genome-scale ortholog alignment using MCScanX v1.0 [89].

Real-time quantitative polymerase chain reaction (RT-qPCR) analysis

The specimens were collected in the wandering stage (n = 3, Sample B1), the prepupal stage (n = 3, Sample B2), the middle pupal stage (n = 3, Sample B3), the late pupal stage (n = 3, Sample B4), and the newly emerged female adults (n = 3, Sample B5). Adult females were quickly dissected into head, thorax, abdomen and leg. Total RNA was extracted and measured by the same method as mentioned above. RNA was reverse-transcribed using Goldenstara™ RT6 Cdna Synthesis Mix (TSINGKE, Beijing, China) according to the manufacturer’s protocol. Then, T5 Fast qPCR Mix (SYBR Green I) (TSINGKE, Beijing, China) was applied to RT-qPCR reactions with a 20 μL reaction volume, following the manufacturer’s protocol. The qPCR reactions were performed with 30 s initial denaturation at 95 °C, followed by 40 cycles at 95 °C for 5 s and 62 °C for 30 s, and a dissociation stage with a melting-curve analysis on an ABI7500 real-time PCR system (Applied Biosystems, Foster City, CA, USA).

Candidate differentially expressed (DEGs) genes were obtained from the stable reference genes [90], and transcriptome data of C. megacephala. The primers were designed with Primer3 Plus software. The relative mRNA levels of candidate genes were represented as folds over the expression levels of reference genes. Based on normalization with the reference genes, the 2^−ΔΔCT method was used to determine expression changes [91]. All experiments were performed in three biological replicates.

Preparation of the recombinant candidate protein and high-affinity polyclonal antibodies

In order to explore the molecular mechanism behind C. megacephala’s behavioral response to VOCs, we identified the gene families of olfactory binding proteins. Sequence similarity and phylogenetic analysis were performed for all OBPs. Subsequently, key candidate genes were screened out. Previous studies have found that DMDS is one of the common VOCs and an active antenna compound used by necrophagous insects to locate breeding sites. Therefore, to further explore the function of OBPs related to the recognition of specific volatile substances, we also performed protein–ligand docking with DMDS for OBPs and obtain the candidate proteins by comparing the affinity scores between OBPs and DMDS. Finally, the coding gene of the candidate protein was confirmed.

To express the candidate protein, the protein coding gene was ligated into the NdeI/EcoRI sites of the pET-28b vector with a His tag on the C-terminal end. The full length of OBP cDNA was then amplified. The HE4 fragments and pTT5 vector were double-digested with EcoRI and HindIII, followed by ligation. The ligation mixture was transformed into E. coli DH5a (AtaGenix, Wuhan, China), and positive clones were sequenced. An appropriate amount of plasmid (1–2 μL) was transferred to the prepared BL21 strain and then cultured overnight at 37 °C. Mono-clone was selected and cultured into 3 mL of LB medium containing the appropriate antibiotic. The supernatant of the bacterial culture was collected and labeled as"NPE". After discarding the excess supernatant, 100 μL of 1 × PBS (containing 8 M urea) was added to the precipitate for re-suspension, which was labeled as"DPE". The sample was mixed with 25 μL of 5 × loading buffer and boiled for 10 min before SDS-PAGE detection. To clearly detect the target protein using SDS-PAGE, it was enriched from the culture medium by His-Tag Purification Resin. This resin effectively increased the concentration of the target protein while removing other impurities. The concentration of the purified protein was then determined at 280 nm using Nanodrop 2000, followed by SDS-PAGE analysis. Ultimately, 5.04 mg of recombinant protein was purified, resulting in a final concentration of 1.50 mg/mL.

To generate high-affinity polyclonal antibodies against the candidate protein, two New Zealand white rabbits were immunized with purified antigens of the candidate protein via injection. Pre-immune serum was collected prior to immunization. Two weeks later, the candidate protein was emulsified in Freund's incomplete adjuvant and used for booster injections, administered twice at approximately two-week intervals. Sera from the rabbits were collected after each immunization, and the immunogenicity and serum antibody titer of the candidate protein were detected by ELISA. Eight days after the third immunization, blood samples were drawn from each rabbit to screen for the highest antibody titer. If the titer of the antiserum exceeded 1:32,000, blood was collected and then antiserum was used for antibody purification. Polyclonal antibody purification was performed using antigen-specific affinity purification (CNBr-Sepharose). Finally, the total yield of rabbit polyclonal antibodies was approximately 3.39 mg, with a concentration of 0.87 mg/ml.

Western blot

Proteins were extracted from various body parts, such as the head, head plus antenna, thorax, or abdomen, using RIPA Lysis buffer (ComWin, Beijing, China), along with a 1% cocktail of protease inhibitor (Roche, Basel, Switzerland) and 1 mM PMSF (Beyotime, Shanghai, China). Protein quantification was performed using a BCA protein detection kit (Solarbio, Beijing, China). The supernatants were electrophoresed on a 15% SDS-PAGE gel and then transferred onto polyvinylidene fluoride (PVDF) membranes. The membranes were blocked by 5% BSA for 1 h at room temperature and then incubated overnight at 4 °C with anti-OBP60 pcAb. After washed three times, the membrane was incubated with an HRP-labeled secondary antibody for 2 h at room temperature. Bands were visualized using the enhanced chemiluminescence (ECL) method in a Tanon 5100 system (Tanon, Shanghai, China).

Results

Genome sequencing and assembly

To survey the genome of C. megacephala, 58.74 Gb of raw data were produced, and 48.01 Gb of clean data were retained (Table S3). The analysis of the k-mer estimated the genome size to be 811.75 ~ 842.84 Mb. The heterozygosity ratio was 1.14% ~ 1.37% (Figs. S2-S3 and Table S4). Moreover, the average genome size of three adult samples, as measured by flow cytometry, was 830.00 Mb (Fig. S4 and Table S5). These results indicate a consistent genome size across the samples, reinforcing the reliability of the genome assembly.

We generated 33.81 Gb of raw PacBio data for genome assembly (Table S6). Finally, the de novo genome assembly was 816.79 Mb in size, with a contig N50, the longest contig and contig number of 1.60 Mb, 11.50 Mb and 1,784, respectively (Table 1). The obtained assembly size is in good correlation with k-mer estimations and flow cytometry evaluation. It is the largest of the published genomes in the blowfly species. Based on the mapping of sequence lengths of the assembled genomes, indicating that there is good continuity (Fig. S5). Meanwhile, the results of the assembly's completeness showed that it covered 97.71% of complete BUSCOs and 81.55% of single-copy BUSCOs. Only 1.89% of BUSCOs were missing, indicating that the majority of conserved genes were fully assembled and the genome assembly was highly reliable (Fig. S6 and Table S7). To assess the accuracy of the final genome assembly, the sequencing data from the genome survey was compared, revealing a mapping rate of 97.42% (Tables S8-S9). Meanwhile, the statistical results of PacBio data showed the mapped rate of 99.57% and the coverage rate of 99.96% (Tables S10), indicating that the assembly is well representative of all the raw data. The analysis showed that the average GC content of the genome was 29.20% (Fig. S7). 57.21 Gb of Hi-C raw data was produced, consisting of 381,430,244 reads. After quality control, 56.82 Gb of clean data were obtained (Table S11). Clean reads were aligned to the reference genome, and 90.52% were mapped. The number of contig N50 is 131 (Tables S12-S13). Finally, pseudochromosomes with a total length of 816.19 Mb were accurately anchored into six chromosomes, representing 99.93% of the draft assembled genome (Fig. 1a), which is identical to the karyotype of six chromosomes based on cytological observation in S. peregrina [36, 92] (Fig. 1b and Table S14). Moreover, the six pseudochromosomes in the assembled genome can be aligned against the genome of D. melanogaster (Fig. 1c), further confirming the reliability of our assembly.

Table 1 Genome assembly of 7 common blowfly species

Fig. 1

Chromosome-level de novo genome assembly and comparative genome analysis of C. megacephala with other species. a Genomic landscape from outer to inner: (I) sizes of 6 pseudochromosomes; (II) DNA transposon content; (III) LTR transposon content; (IV) GC content (%); (V) gene density. Densities are calculated in 500 Kb windows. The dorsal image of an adult female is shown in the circle. b Contig contact matrix of the assembled genome. The color bar on the right shows the density of Hi-C interactions from red (high) to white (low), which are the number of contact links at the 500-kb resolution. c Chromosome collinear blocks between C. megacephala and D. melanogaster genome. The best match across the two species is linked by lines with the same color. The pseudochromosome of C. megacephala is labeled “Chr1-6”, and “X, 2L, 2R, 3L, 3R, 4Y” represent the chromosomes of D. melanogaster, respectively. d Venn diagram shows single-copy orthologs and species-specific genes in C. megacephala and other flies. Single: Number of single-copy orthologs; The number represents the sum of exclusive and unclustered genes. Cmeg: Calliphora megacephala, Agra: Aldrichina grahami, Sper: Sarcophaga peregrina, Sbul: Sarcophaga bullata, Preg: Phormia regina, Mdom: Musca domestica, Lser: Lucilia sericata, Lcup: Lucilia cuprina, Dmel: Drosophila melanogaster, Cvic: Calliphora vicina, Cruf: Chrysomya rufifacies

Functional annotation of protein-coding genes

Transcriptome sequencing was performed to assist in functional annotation in the assembled genome. A total of 96.77 Gb of clean data were retained after quality control, including an average of 44.00 Mb of clean reads from each sample (Table S15). 93.00% to 96.12% of the clean reads could be mapped onto the assembled genome (Table S16), indicating a high quality of sequencing. FPKM was calculated and standardized for the analysis of gene expression (Table S17). A total of 16,589 genes were identified, with an average number of 4.11 exons per gene. The average transcript length, average CDS length, and average exon length per gene were 7,624.20 bp, 1,483.97 bp, and 360.81 bp, respectively (Table S18). A total of 16,054 protein-coding genes were annotated with potential functions, accounting for 96.77% of all genes in the assembled genome. We identified 12,062 genes that showed homology to proteins in the PFAM database. A total of 13,950 genes were assigned to the GO database, as well as 15,478 genes in the Uniprot database. Additionally, 11,484 genes were annotated in the NR database (Fig. S8 and Table S19).

The results of the de novo and homology-based predictions showed that 534.24 Mb of repetitive sequences were identified, covering 65.41% of the assembled genome (Table S20). DNA transposons (189.90 Mb) represented the most abundant transposable elements (TEs), accounting for 23.25% of the genome (Table S21). Furthermore, 77 miRNAs, 990 rRNAs, 101 snRNAs, and 1,519 tRNAs were identified in the assembled genome (Table S22).

Gene family identification and phylogeny analysis

A total of 10,282 gene families were identified in the assembled genome, covering 15,764 genes. Among these, 146 gene families were unique to C. megacephala. Besides, 730 unclustered genes were identified (Table S23). We then identified 1,598 single-copy orthologs to construct phylogenetic trees, including 1,197 specific genes for C. megacephala (Fig. 1d and Table S24), whilst the significant variation in species-specific genes further indicates a high degree of specificity between species. Phylogenetic analysis indicated that ten fly species were clustered together into a large branch and strongly supported by bootstrap (ML bootstrap percentage, BP = 100). Additonally, C. megacephala and C. rufifacies were clustered more closely than the other species. As an outgroup taxon, D. melanogaster (Diptera: Drosophilidae) was clearly separated (Fig. S9). Moreover, among the Calliphoridae family, the estimated divergence time between them was 23.2 million years ago (Mya) (95% HPD: 21.3–25.7 Mya) within the Late Paleogene epoch. The diversification of C. megacephala and C. rufifacies occurred at 8.90 Mya (95% HPD: 6.30–11.70 Mya) (Fig. S10).

Gene family expansion and contraction

This approach revealed 3,324 contracted gene families in C. megacephala, and 2,295 gene families had a higher degree of expansion (Fig. 2 and Table S25), which exhibited significant expansion and contraction compared to other closely related species. The corresponding genes were identified from these gene families, which were used for enrichment analysis of KEGG and GO, respectively. The expanded genes were primarily enriched in GO pathways, notably including the structural constituent of the eye lens, ligand-gated channel activity, passive transmembrane transporter activity, defense response to Gram-positive bacteria, antibacterial humoral response, response to bacteria, and sensory perception of taste. The KEGG pathways mainly included the drug metabolism—cytochrome P450, biosynthesis of secondary metabolites, glycine, serine, threonine metabolism, and glutathione metabolism, etc. (Tables S26-S29). Besides, a total of 322 positively selected genes were identified in C. megacephala. Subsequently, KEGG and GO enrichment analyses were performed (Tables S30-S31). The pathways mainly included the DNA-binding transcription factor activity, RNA polymerase II-specific, GTPase binding, Hippo signaling pathway – fly, etc.

Fig. 2

Phylogenetic tree of C. megacephala and other species (A. grahami, S. peregrina, S. bullata, P. regina, M. domestica, L. sericata, L. cuprina, C. vicina, C. rufifacies, D. melanogaster) based on the analysis of contracted and expanded gene families. The number near each branch indicates the number of significantly expanded (green) and contracted (red) gene families for each clade. The black numbers show the divergence times. Blue: the number of gene families unchanged in size across all branching gene families. Orange: the number of gene families changed in any branch gene family

Insecticide resistance genes

In order to gain a better understanding of insecticide resistance, the insecticide resistance genes of C. megacephala were identified in the assembled genome. Three gene families are involved in the detoxification of xenobiotics in insects, primarily including the cytochromes P450 monooxygenases (CYPs) and glutathione S-transferases (GSTs) [96]. The largest group is the cytochrome P450 s, in which a total of 121 genes were identified. This represents a significant expansion of P450 s compared to D. melanogaster (86) [39], but a contraction compared to Musca domestica (146) [38]. Additionally, a total of 42 GST genes were predicted from the assembled genome, which is larger than in other species, indicating that this family has significantly expanded. Moreover, insect Cys-loop ligand-gated ion channels (cysLGICs) play an important role in the new targets of insecticides, whose members are essential in mediating chemical synaptic transmission [97, 98]. For the CysLGICs superfamily of C. megacephala, a total of 24 genes were identified, which is similar to other insects, such as M. domestica (23) and D. melanogaster (23).

Identification of OBP gene family and phylogenetic analysis

The results showed that a total of 69 OBPs were identified in the assembled genome of C. megacephala, compared to 41 OBPs in S. peregrina. In addition, the number of OBPs in D. melanogaster and M. domestica was 52 and 87, respectively. As expected, C. megacephala shared the roughly conserved members of these families with D. melanogaster (Fig. 3 and Table S32). Cluster analysis of OBP family members in four species revealed that DmelOBPs, MdomOBPs, and SperOBPs were highly similar to those of CmegOBPs (Fig. S11 and Table S33). Most OBPs clustered together, indicating the presence of species-specific duplicates as well as species-specific single-copy genes during evolution. Transcriptome analysis further indicated that the expression level of CmegOBP99b was the highest among all OBPs during metamorphosis from larva to adult (Fig. 4a). Additionally, by comparing its molecular docking affinity with DMDS (Fig. S12), we conclude that CmegOBP99b may be the DMDS-binding OBP of C. megacephala.

Fig. 3

Identification and phylogenetic analysis of odorant-binding proteins (OBPs) family. Phylogenetic tree of OBPs family, which was constructed using IQTREE with Maximum Likelihood method. C. megacephala (Cmeg) and D. melanogaster (Dmel) are marked in red and green, respectively. Bootstrap values greater than 70% (1000 replications) were displayed

Fig. 4

Expression of recombinant CmegOBP99b and development of polyclonal antibodies against CmegOBP99b. a Expression profiles of OBPs family. The color scale on the right shows the FPKM values for each developmental stage from red (high) to white (low), standardized by min–max normalization. B1, B2, B3, B4 and B5 represent the wandering stag, the prepupal stage, the middle pupal stage, the late pupal stage and the newly emerged adults, respectively. b Developmental expression patterns of CmegOBP99b in different developmental stage of C. megacephala. c Tissue-specific expression patterns of CmegOBP99b in adult females. d The protein was mainly expressed in the inclusion body since the expression plasmid contained CmegOBP99b natural signal peptide. MW: protein marker; 1: inclusion body protein; 2: Quality control diagram of amplified and purified protein (SDS-PAGE). e Western blot assay showed the results of CmegOBP99b rabbit polyclonal antibodies, 1: CmegOBP99b rabbit polyclonal antigen protein (16.81 KDa). f WB analysis of CmegOBP99b polyclonal antibodies in different tissues. 1–3: abdomen, 4–6: thorax, 7–9: head plus antenna, 10–12: head

Validation of transcriptome data and expression of recombinant CmegOBP99b

To validate the transcriptome result, CmegOBP99b was selected for RT-qPCR analysis. The primers used for RT-qPCR were shown (Table S34). These results are consistent with the RT-qPCR results (Fig. 4b). The expression level of CmegOBP99b was highest in the prepupal stage, then down-regulated after entering the pupal stage, and finally up-regulated after emerging into adults. Moreover, the tissue-specific expression pattern of CmegOBP99b in females was performed (Fig. 4c). CmegOBP99b was significantly expressed in the female’s head compared to other parts.

Subsequently, CmegOBP99b was amplified and cloned into pET-28b to generate the expression plasmid, which was transfected into BL21 strain and cultured. Afterwards, purified CmegOBP99b proteins were collected. The SDS-PAGE analysis showed that the recombinant CmegOBP99b protein was mainly expressed in the inclusion body since the expression plasmid contained CmegOBP99b natural signal peptide. The molecular weight of the recombinant CmegOBP99b was consistent with the predicted value (16.81 KD) (Fig. 4d). In order to produce more efficient polyclonal antibodies, after three consecutive immunization, ELISA assay showed that the serum titer was greater than 64,000 after the third immunization (Table S35). The purity of CmegOBP99b polyclonal antibodies was identified by western blot analysis (Fig. 4e and Fig. S13). To determine the expression of the CmegOBP99b protein in different body parts of C. megacephala, western blot analysis indicated that it was highest in the head and head plus antenna, and lowest in the thorax and abdomen (Fig. 4d and Fig. S14).

Discussion

It is well known that the primary task of necrophagous flies colonizing on the decomposed corpses is to lay eggs in order to reproduce. The process of flies laying eggs on carrion is driven by a series of intricate physiological mechanisms, jointly regulated by external chemical cues, neural signal transmission and integration [99]. During the decomposition process, carrion releases a series of volatile organic compounds (such as ammonia, amines, sulfides, etc.), which are keenly detected by the fly’s sensory organs—particularly the olfactory receptors located on its antennae and mouthparts [100]. Afterwards, the detected chemical signals are transmitted via the olfactory receptors to the fly’s central nervous system, where the information is integrated, subsequently activating neural circuits related to reproductive behavior and prompting the fly to decide to seek an appropriate site for oviposition [101]. Therefore, olfaction plays a crucial role in the detection and interpretation of various semiochemicals or odorants in the environment, guiding animals toward essential behaviors such as locating food, avoiding predators, mating, and nurturing offspring [102]. A complex and sensitive olfactory system has been developed during the long-term evolution. Like many animals, insects rely on their sense of smell to navigate through the environment towards food sources and mates [103].

Insect feeding behavior involves a wide range of activities, such as initial activation, orientation, recognition, and feeding [104]. For example, necrophagous flies are known to colonize and breed on decomposing organic matter, with olfactory cues playing a pivotal role in host choice selection [33, 105, 106]. Phylogenetic analyses have shown that the adaptive evolution in insects has led to diverse families of chemical receptors, with the ability to detect smells and tastes depending on a wide array of chemoreceptors and related proteins. These proteins are encoded by at least four major gene families, among which OBPs play a central role [107], and which are small, globular, secreted proteins that deliver hydrophobic odorants to the receptors on the membranes of sensory neurons within the sensory organs, primarily located on the antennae [108]. Given the importance of olfaction in insect behavior, which rely on olfactory cues for detecting cadaveric volatiles, it is critical to improve genomic resources for C. megacephala. A high-quality genome assembly is essential to gain a better understanding of the molecular mechanisms underlying olfactory response in C. megacephala.

We hereby report the first chromosome-level draft genome of C. megacephala. Until now, the genomes of carcass-associated flies have been rarely published, mainly including C. rufifacies, S. peregrina, M. domestica, L. sericata, L. cuprina, P. regina, C. vicina, S. bullata, and A. grahami [35,36,37,38, 95, 109]. Firstly, for the estimation of C. megacephala's genome size, the results indicate that the k-mer depth distribution estimation was highly consistent with that of flow cytometry. The de novo genome assembly was 816.79 Mb in size, with a contig N50, the longest contig and contig number of 1.60 Mb, 11.50 Mb and 1,784, respectively. The Hi-C sequence data was used to assist in the draft genome assembly. A total of 1,742 contigs were anchored to six chromosomes, and the final chromosome-level draft genome assembly was 816.19 Mb with a scaffold N50 of 152.02 Mb. Furthermore, 16,054 protein-coding genes were annotated in the assembled genome. Comparative genomic analysis was then performed with other published dipteran flies to enhance our understanding of adaptive evolution. Functional enrichment analyses of expanded genes, positively selected genes, and shared genes indicated that, in addition to the pathways (ligand-gated ion channel activity and ion gated channel activity) that have significantly expanded, the pathways of the passive transmembrane transporter activity, protein methyltransferase activity, and DNA-binding transcription factor activity were also significantly enriched. Although this study shows that they are key components of cadaver localization, further studies are needed to determine the role of specific genes associated with these pathways, and to analyze their specific functions in reaching cadavers earlier than other flies. Moreover, enrichment analysis revealed a significant expansion of xenobiotics metabolism, specifically cytochrome P450 and drug metabolism-cytochrome P450 in the C. megacephala lineage. These pathways appear to be associated with insecticide resistance and the ability to thrive in decomposed corpses.

Striking similarities in olfaction are observed across a wide range of insect species, suggesting the possibility of an optimal solution for detecting and distinguishing odors [110]. Therefore, studying the orthologous OBPs in D. melanogaster could provide valuable insights into the biological roles of candidate OBPs [32]. Comparative genomic analysis of the OBPs gene families in four species revealed a highly dynamic evolutionary pattern, characterized by significant gene gains and losses. The number of OBPs varies widely across these species, reflecting their diversity. Additionally, insect feeding behavior appears to be closely linked to species-specific OBP duplications [111]. Interestingly, in this study, most OBP members exhibit clear orthologous relationships not only in the Drosophila genomes but also across three other species, indicating that the OBP gene family is highly conserved during evolution. This observation suggests the presence of functional constraints that maintain the clusters [112]. Only a few OBPs were not clustered, suggesting the occurrence of species-specific gene duplications and single-copy genes. It is well established that OBPs capable of binding the same volatile substances should cluster into the same orthogroup, but other OBPs are also clustered together. Based on these results, it is speculated that OBPs may recognize odorants in a combinatorial manner. However, the number of OBPs involved in olfactory recognition can also affect the sensitivity of the olfactory sense to an odorous substance, thus altering the preferences of insects [113]. Additionally, although we concluded that CmegOBP99b might play a critical role in the recognizing DMDS based on the molecular docking affinity and could enhance the sensitivity of C. megacephala in detecting cadaveric VOCs, it is necessary to use affinity analysis, such as microcalorimetry can further confirm the binding affinities between CmegOBP99b and DMDS.

Wang et al. demonstrated that the DmelOBP99 group is associated with the recognition of olfactory responses to benzaldehyde and acetophenone [114]. DmelOBP99b is sensitive to and may affect the nutritional and reproductive status of both females and males [115]. Swarup et al. found that the effects of suppressing OBP expression on behavioral responses to odorants were sexually dimorphic. For example, the disruption of DmelOBP99b expression altered behavioral responses to 2-ethylpyrazine in females and to acetophenone, benzaldehyde, citral, and d-carvone in both sexes. Additionally, the suppression of multiple OBPs often affects responses to specific odorants [116]. Swarup et al. further showed that RNAi-mediated inhibition of DmelOBP99b led to a significant increase (47%) in the intake of berberine and papaverine. Additionally, the inhibition of DmelOBP99b gene expression affects feeding behavior differently in males and females This finding is consistent with the proposed roles of OBPs in transporting taste substances to bitter receptors and isolating them to restrict their access to these receptors [117]. These results suggest that DmelOBP99b is essential in mediating olfactory behavioral responses. In this study, the evolutionary analysis revealed that DmelOBP99b and CmegOBP99b exhibited a high degree of homology, as they clustered together in the same branch of the evolutionary tree. Based on the identification of the specific attractants for C. megacephala, CmegOBP99b, which exhibits a strong binding affinity for DMDS, was selected for further investigation. Moreover, olfactory proteins have been demonstrated to influence various aspects of insect biology, including nutrient uptake, lifespan, and behavioral changes during developmental stages [118, 119]. The transcriptional regulation of OBPs affects the olfactory preference of insects, as evidenced by the differential expression of OBPs in various body parts and developmental stages [120]. In this study, the relative expression of CmegOBP99b is consistently highest during the metamorphosis, but the discrepancy between qPCR analysis and transcriptome data may be due to which these were not conducted on the same batch of samples. Therefore, it still requires further exploration. Moreover, to determine the expression of CmegOBP99b protein in various body parts of female adults, RT-qPCR analysis indicated that CmegOBP99b was significantly expressed in the female’s head compared to other parts. Furthermore, western blot analysis also revealed that the highest expression of CmegOBP99b protein was observed in the head. These results indicate that a higher expression level of OBPs can improve the sensitivity of C. megacephala in recognizing specific cadaveric VOCs.

Conclusions

This study highlights the crucial role of olfaction in the survival and behavior of necrophagous insects, emphasizing the complexity and adaptability of the olfactory system developed over long-term evolution. By providing the first chromosome-level draft genome of C. megacephala, we not only gained a deeper understanding of the molecular mechanisms underlying this species but also provided essential data for comparative genomic studies. The investigation into CmegOBP99b particularly underscores its significance in recognizing VOCs associated with decomposition, indicating that OBPs play a vital regulatory role in olfactory recognition, host localization, and mating behavior in insects. Future research should focus on exploring the expression patterns and functional variations of OBPs and related genes under different environmental conditions, elucidating their roles in ecological adaptation, behavioral shaping, and responses to external environmental changes. Additionally, with advancements in genomics and molecular biology techniques, coupled with bioinformatics analyses, we will gain broader insights into the mechanisms of insect olfaction. These studies not only hold substantial significance for insect ecology and evolutionary biology but also provide a theoretical foundation for developing novel biological control strategies, particularly in managing pest populations associated with decomposition.