Basic genome information for Z. cucurbitae (Coquillett), B. dorsalis (Hendel), and C. capitata (Wiedemann) was obtained from the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/).

The genomic data of Z. cucurbitae (Coquillett), B. dorsalis (Hendel), and C. capitata (Wiedemann) were used to identify homologous genes using blastp in combination with OrthoMCL. After obtaining homologous and specific gene families, the number of genes contained in each gene family (the number of gene family copies) was counted in each of the three fruit fly species, and the gene families specific to the three fruit fly species were identified and subjected to functional annotation, KEGG enrichment analysis, and GO enrichment analysis. The gene families corresponding to homologous genes in the three fruit fly species were classified into two gene sets, single-copy homologous genes, and multi-copy homologous genes, and the Ka/Ks values of single-copy homologous genes were calculated, and the single-copy homologous genes were classified into the strong positive select (SP), weak positive select (WP), and negative select (NG) genes according to the Ka/Ks values, and KEGG enrichment analysis and GO enrichment analysis were performed on the SP single-copy homolog genes.

Within an organism, different genes coordinate their biology with each other, and pathway-based analysis helps to further understand the biological functions of genes. KEGG is the main public database on pathways. The pathway significant enrichment analysis uses the KEGG pathway as the unit and applies a hypergeometric test to identify pathways that are significantly enriched in species-specific genes when compared to the entire genome. The p-value of this hypothesis test is calculated as follows: p = 1 − ∑ i = 0 m − 1 ( M i ) ( N − M n − i ) ( N n ) ,

where N is the number of genes with pathway annotation in all unigenes, n is the number of specifically expressed genes in N, M is the number of genes annotated with a specific pathway in all unigenes, and m is the number of specific genes expressed with a specific pathway annotation. After correction for multiple testing, pathways with Q-value ≤ 0.05 were selected as those significantly enriched in specific genes. The q-value and FDR are similar, both being a correction for the p-value. The Q-value of a particular hypothesis test is the minimum of the FDR, and the hypothesis test can be considered significant under these FDR. Enrichment by pathway significance identifies the most important biochemical metabolic pathways and signal transduction pathways involved in specific genes.

Gene Ontology (GO) is an international standardized gene function classification system that provides a dynamically updated set of controlled vocabulary to comprehensively describe the properties of genes and gene products in an organism. GO has a total of three ontologies, which describe the molecular function of a gene, the cellular component, and the biological process involved. Each term corresponds to an attribute.

GO function significant enrichment analysis provides GO function entries that are significantly enriched in specific genes after comparison with reference genes and screens which biological functions are significantly associated with specific genes. This analysis first maps all specific genes to each term of the Gene Ontology database (http://www.geneontology.org/) and counts the number of genes per term to obtain a list of genes with a particular GO function and a count of the number of genes. A hypergeometric test was then applied to identify GO entries that were significantly enriched in differentially expressed genes compared to the whole genomic background, and the p-value of this hypothesis test was calculated as follows: p = 1 − ∑ i = 0 m − 1 ( M i ) ( N − M n − i ) ( N n ) ,

where N is the number of genes with GO annotation in all genes, n is the number of species-specific genes in N, M is the number of genes annotated with a particular GO term in all genes, and m is the number of species-specific genes annotated with a particular GO term. The calculated p-value was corrected by Bonferroni with a threshold value of corrected p-value ≤ 0.05, and the GO term meeting this condition was defined as a GO term significantly enriched in species-specific genes.

In genetics, Ka/Ks or dN/dS indicates the ratio between the rate of non-synonymous substitutions (Ka) and the rate of synonymous substitutions (Ks). This ratio can determine whether there is a selective pressure acting on this protein-coding gene. It is generally accepted that synonymous mutations are not subjected to natural selection, whereas non-synonymous mutations are subjected to the action of natural selection. In evolutionary analysis, it is of interest to know the rate at which synonymous and non-synonymous mutations occur. The following parameters are commonly used: the frequency of synonymous mutations (Ks), the frequency of non-synonymous mutations (Ka), and the ratio of the rate of non-synonymous mutations to the rate of synonymous mutations (Ka/Ks). Single-copy homologous gene protein sequences from different species or subspecies were subjected to multiple sequence alignment by using muscle (http://www.drive5.com/muscle), and then the Ka ratio between each single-copy homologous gene was calculated based on the nucleic acid sequence corresponding to the protein sequence. The Ka, Ks, and Ka/Ks values between homologous genes were calculated using the KaKs_Calculator Toolbox (version 2.0).

The basic information of the genomic data of the three fruit fly species is shown in Table 1.

For the genomes of the three fruit fly species, blastp was used in conjunction with OrthoMCL to identify homologous genes among the three species. The sequences were matched two by two, and pairs with an E-value less than 1e-5 and a coverage of more than 30% were considered homologous genes in the Blast comparison. OrthoMCL was used to assign homologous genes to the same family. As shown in Figure 1, there were 4,656 homologous genes and 955 gene families in the three fruit fly species, 211 and 69 homologous genes and gene families in Z cucurbitae (Coquillett) and B. dorsalis (Hendel), 37,415 and 10,547 homologous genes and gene families in Z cucurbitae (Coquillett) and C. capitata (Wiedemann), and 61 and 24 homologous genes and gene families in B. dorsalis (Hendel) and C. capitata (Wiedemann). Z cucurbitae (Coquillett) had 2,754 specific genes and 1,881 gene families, B. dorsalis (Hendel) had 19701 specific genes and 19701 gene families, and C. capitata (Wiedemann) had 3,160 specific genes and 2,718 gene families.

The unigenes of Z cucurbitae (Coquillett), B. dorsalis (Hendel), and C. capitata (Wiedemann) were annotated to the KEGG database for comparison, and their specific gene KEGG enrichment analysis bubble maps were obtained to identify the pathways that were significantly enriched in specific genes. The highest gene abundance in Z. cucurbitae (Coquillett) was found for the Toll and Imd signaling pathway with about 50 genes, followed by MicroRNAs in cancer and influenza A, all three with Q values less than 0.0005. Among all metabolic pathways in Z. cucurbitae (Coquillett), the largest enrichment factor was neomycin, kanamycin, and gentamicin biosynthesis with an enrichment factor of about 0.7, followed by hematopoietic cell lineage, with an enrichment factor of about 0.46 (Figure 2). The highest gene abundance was found in B. dorsalis (Hendel) with human papillomavirus infection, enriched with 331 genes, followed by the MAPK signaling pathway—fly MAPK and viral carcinogenesis. Among all metabolic pathways in B. dorsalis (Hendel), the largest enrichment factor is the TGF-beta signaling pathway with an enrichment factor of about 1, followed by oocyte meiosis and the Hippo signaling pathway—fly (Figure 3). The highest KEGG enrichment results for C. capitata (Wiedemann) were found for protein digestion and absorption (39 genes), followed by viral carcinogenesis and the Toll and Imd signaling pathway. Among all metabolic pathways in C. capitata (Wiedemann), the largest enrichment factor is circadian rhythm—fly with an enrichment factor of about 0.42, followed by ECM-receptor interaction with an enrichment factor of about 0.34 (Figure 4).

The results of GO enrichment analysis for genes specific to Z. cucurbitae (Coquillett), B. dorsalis (Hendel), and C. capitata (Wiedemann) are shown in Figures 5–7. In the broad category of biological processes, all three species were enriched for more genes in metabolic processes, intracellular processes, and single-tissue processes, and Z. cucurbitae (Coquillett) was enriched most in the broad category of cellular components in the pathways of membranes, cells, and cellular components. B. dorsalis (Hendel) and C. capitata (Wiedemann) were most enriched in the cellular component major category for cellular, cellular component, and organelle pathways, and in the molecular function major category, all three species were more catalytically active and binding enriched for genes.

The genes can be classified according to the Ka/Ks values of single-copy homologous genes. Z. cucurbitae (Coquillett) and B. dorsalis (Hendel) have 508 homologous unigene pairs (Figure 8); and two unigenes have Ka/Ks ratios greater than 1 (Table 2), and the remaining 506 have Ka/Ks ratios less than or equal to 1. Z. cucurbitae (Coquillett) and C. capitata (Wiedemann) have 519 homologous unigene pairs (Figure 9); and two unigenes have Ka/Ks ratios greater than 1 (Table 3). B. dorsalis (Hendel) and C. capitata (Wiedemann) have 517 homologous unigene pairs with Ka/Ks ratios less than 1 (Figure 10). Because most non-synonymous mutations are deleterious and single-copy homologous genes have a low tolerance for deleterious mutations on them due to the lack of copies, most genes have small Ka/Ks values. To identify positively selected genes that are neutralized (masked) by a large number of deleterious mutation sites, genes with Ka/Ks > 1 are defined as SP, genes with Ka/Ks before 0.5 and 1 are WP, and genes with Ka/Ks < 0.1 are NG (purifying selection).

The KEGG enrichment analysis showed that only one gene was enriched among the single-copy homolog genes of Z. cucurbitae (Coquillett) and B. dorsalis (Hendel), and the enriched pathways were fatty acid degradation, fatty acid metabolism, and lipid metabolism (Table 4), and the KEGG enrichment results of the single-copy homolog genes of Z. cucurbitae (Coquillett) and C. capitata (Wiedemann) were similar to those of Z. cucurbitae (Coquillett) and B. dorsalis (Hendel) (Table 5). GO enrichment analysis showed that genes were enriched in the three broad categories of biological processes, molecular functions, and cellular components, and in the broad category of biological processes of Z. cucurbitae (Coquillett) and B. dorsalis (Hendel), genes were enriched in the GO:0044699-single-organism process. Among the molecular functional categories, genes are enriched in GO: 0005488-binding, GO:0003824-catalytic activity, and GO:0005198-structural molecule activity, while the cellular components are enriched in GO:0005623-cell, GO:0044464-cell part, and GO:0043226-organelle (Table 6). In the biological process category of Z. cucurbitae (Coquillett) and C. capitata (Wiedemann), genes are enriched in the GO:0008152-metabolic process and GO:0044699-single-organism process; in the molecular function category, genes are enriched in GO:0005488-binding and GO:0003824-catalytic activity; and the cellular component is enriched in GO:0005623-cell, GO:0044464-cell part, and GO:0032991-macromolecular complex (Table 7). In the biological process category of B. dorsalis (Hendel) and C. capitata (Wiedemann), genes are enriched in the GO:0008152-metabolic process and GO:0044699-single-organism process; in the molecular function category, genes are enriched in GO:0005488-binding and GO:0003824-catalytic activity; and in the cellular component, genes are enriched in GO:0005623-cell and GO:0044464-cell part (Table 8).