Literature DB >> 31702774

Phylogenomics Identifies an Ancestral Burst of Gene Duplications Predating the Diversification of Aphidomorpha.

Irene Julca1, Marina Marcet-Houben1, Fernando Cruz2, Carlos Vargas-Chavez3, John Spencer Johnston4, Jèssica Gómez-Garrido2, Leonor Frias2, André Corvelo2,5, Damian Loska1, Francisco Cámara1, Marta Gut2,6, Tyler Alioto2,6, Amparo Latorre3,7, Toni Gabaldón1,6,8.   

Abstract

Aphids (Aphidoidea) are a diverse group of hemipteran insects that feed on plant phloem sap. A common finding in studies of aphid genomes is the presence of a large number of duplicated genes. However, when these duplications occurred remains unclear, partly due to the high relatedness of sequenced species. To better understand the origin of aphid duplications we sequenced and assembled the genome of Cinara cedri, an early branching lineage (Lachninae) of the Aphididae family. We performed a phylogenomic comparison of this genome with 20 other sequenced genomes, including the available genomes of five other aphids, along with the transcriptomes of two species belonging to Adelgidae (a closely related clade to the aphids) and Coccoidea. We found that gene duplication has been pervasive throughout the evolution of aphids, including many parallel waves of recent, species-specific duplications. Most notably, we identified a consistent set of very ancestral duplications, originating from a large-scale gene duplication predating the diversification of Aphidomorpha (comprising aphids, phylloxerids, and adelgids). Genes duplicated in this ancestral wave are enriched in functions related to traits shared by Aphidomorpha, such as association with endosymbionts, and adaptation to plant defenses and phloem-sap-based diet. The ancestral nature of this duplication wave (106-227 Ma) and the lack of sufficiently conserved synteny make it difficult to conclude whether it originated from a whole-genome duplication event or, alternatively, from a burst of large-scale segmental duplications. Genome sequencing of other aphid species belonging to different Aphidomorpha and related lineages may clarify these findings.
© The Author(s) 2019. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution.

Entities:  

Keywords:  Aphidomorpha; aphids; gene duplication

Year:  2020        PMID： 31702774      PMCID： PMC7038657          DOI： 10.1093/molbev/msz261

Source DB:  PubMed          Journal:  Mol Biol Evol        ISSN： 0737-4038            Impact factor:   16.240


Introduction

Large-scale gene duplication, including whole-genome duplication (WGD), is a very common phenomenon in eukaryotic genomes. Bursts of gene duplications are considered a major source of evolutionary innovation and have been associated with the increase in biological complexity and adaptive radiations of species (Zhang 2003). In particular, large-scale gene duplications, generally associated with WGDs, have been reported for many eukaryotic lineages including plants (Van de Peer et al. 2017), fungi (Marcet-Houben and Gabaldón 2015), and animals (Taylor et al. 2001). Although large-scale duplication seems less pervasive in animals than in plants, a growing number of studies report such events in animals. Among other lineages, putative WGDs have been described at the base of vertebrates (Ohno 1970; Dehal and Boore 2005; Putnam et al. 2008), and in several lineages of fish (Christoffels et al. 2004; Glasauer and Neuhauss 2014), amphibians (Mable et al. 2011; Session et al. 2016), and arthropods (Jacobson et al. 2013; Kenny et al. 2016; Schwager et al. 2017; Li et al. 2018). Aphids belong to the infraorder Aphidomorpha that includes three families: Aphididae, Adelgidae, and Phylloxeridae (Favret 2013; Nováková et al. 2013; Blackman and Eastop 2000). Aphids and related (Aphidomorpha) species (Becker-Migdisova and Aizenberg 1962) are hemipteran insects that feed on plant sap (Tjallingii 1995). This specialized diet, rich in carbohydrates but poor in nitrogen compounds, has resulted in several adaptations including the establishment of tight relationships with bacterial endosymbionts (Scarborough et al. 2005; Moya et al. 2008; von Dohlen et al. 2017). There are more than 5,000 described aphid species, of which, about 450 have been collected from crop plants, and 100 are considered of significant economic importance (Van Emden and Harrington 2017). Genomes of several aphid species of agricultural interest have been sequenced, including Acyrthosiphon pisum, Myzus persicae, Diuraphis noxia, Aphis glycines, and Sipha flava (International Aphid Genomics Consortium 2010; Nicholson et al. 2015; Mathers et al. 2017; Wenger et al. 2017). However, except for S. flava (subfamily Chaitophorinae), the sequenced aphids belong to a single subfamily, Aphidinae, limiting our understanding of the genomic diversity in this group of insects. Remarkably, most genome analyses in these species have revealed an important number of paralogous sequences and expanded gene families, including amino acid transporters, odorant and gustatory receptor genes, miRNA-specific dicer-1, ago1 genes, and pasha, among others (Smadja et al. 2009; Huerta-Cepas et al. 2010; Jaubert-Possamai et al. 2010; Duncan et al. 2016; Mathers et al. 2017). However, the close relatedness of the sequenced species provides little resolution to the phylogenetic placement of the duplication events, particularly the ancestral ones. Recent studies have focused on assessing patterns of sequence and expression divergence among recently duplicated genes in A. pisum (Fernández et al. 2019) or M. persicae (Mathers et al. 2017). They have also inspected the distribution of old and young A. pisum paralogs along chromosomes, by categorizing the age of genes that are best-reciprocal hits of each other based on the amount of synonymous substitutions (Li et al. 2019). However, we still lack a proper understanding of when the ancestral duplications occurred, and whether they can be linked to phenotypic innovations shared by aphids or related species. To better assess the origin of the paralogous genes of aphids we sequenced the genome of Cinara cedri (Lachninae subfamily, tribe Eulachnini), the first representative genome from an early-branching lineage of the Aphididae family. Cinara species (and most Lachninae) are particular among aphids as they feed on conifers (gymnosperms), whereas all the other genome-sequenced aphids feed on angiosperms. Another clear difference between the Lachninae and the rest of aphids is that two co-obligate endosymbionts (Buchnera aphidicola, Serratia symbiotica) are present in this group, whereas only B. aphidicola is obligate for the rest of aphids (Latorre and Manzano-Marín 2017). We used a phylogeny-based approach (Huerta-Cepas and Gabaldón 2011) to provide the relative timing of aphid duplications in a phylogenetic framework that includes 21 other fully sequenced genomes and two transcriptomes. Our results provide compelling evidence for an ancestral wave of gene duplications, whose origin predates the diversification of all sequenced aphids, adelgids, and phyloxerids, but are subsequent to their divergence from the Coccoidea lineage, ∼106–227 Ma.

Results and Discussion

Genome Sequence of C. cedri

The haploid genome sizes for C. cedri and two other Lachninae species (C. tujafiina and Tuberolagnus salignus, tribes Eulachnini and Tuberolachnini, respectively) were measured using flow cytometry (Johnston et al. 2019) which resulted in estimates of ∼592, 713, and 494 Mb, respectively. For reference, the genome size of A. pisum is 520.8 Mb (International Aphid Genomics Consortium 2010). We used an Illumina pair-end sequencing approach to produce a draft assembly of the C. cedri genome (see Materials and Methods). A rough estimate of the genome size obtained by dividing the total number of 17-mers by the peak 17-mer coverage results in an estimate of 508.6 Mb (supplementary fig. S1, Supplementary Material online), slightly smaller than the flow cytometry estimate. However, the K-mer profiles indicated an appreciable amount of repeated sequences, which makes the assembly from short reads challenging. To obtain a more precise estimate, we used GenomeScope v1.0 (Vurture et al. 2017) and fit the previous K-mer profile to a mixture model. This provided a haploid genome size estimate of 399.76 Mb, which was used to guide our assembly strategy. This analysis also inferred the amount of unique (223.34 Mb) and repetitive (175.4 Mb) content. Separate assemblies, exploring different K-mer sizes, were done with ABySS v1.5.2 (Simpson et al. 2009), and later merged with ASM (Cruz et al. 2016). The continuity of the merged assembly was improved through several rounds of scaffolding, first with ABySS and later with SSPACEv3.0 (Boetzer et al. 2011). Gaps were closed with GapFiller (Boetzer and Pirovano 2012). The length of the final assembly (see Materials and Methods) is 396.03 Mb, and its contig and scaffold N50 are 104,784 bp and 1.23 Mb, respectively. The gene completeness of our assembly is high, as evaluated by BUSCO v3.0.2 (93.9% of 1,658 single-copy, conserved genes in insecta_odb9 data set were present) and CEGMA (100% of 248 eukaryotic core genomes) (Parra et al. 2007; Simão et al. 2015). Notably, 2.5% of the BUSCO genes were duplicated in our assembly. The postassembly K-mer analysis (Mapleson et al. 2017) suggests that these are real paralogs and not the result of assembly artifacts (supplementary fig. S2, Supplementary Material online). The final protein-coding annotation (see Materials and Methods) resulted in 16,996 genes, whose 24,835 transcripts (1.46 transcripts/gene) encode 22,503 unique protein products. Attempts to detect selenoprotein genes with selenoprofiles (Santesmasses et al. 2018) failed, which indicates that the previously described loss of selenoproteins in some aphids (International Aphid Genomics Consortium 2010; Mariotti et al. 2015) is ancient, and had already occurred at the base of the Aphididae lineage. Similarly, the immune repertoire in C. cedri resembles that of other sequenced aphids, which indicates that the reported streamlining of the immune system in aphids (Gerardo et al. 2010) appears at the base of Aphididae lineage (supplementary table S1, Supplementary Material online). Previous aphid genome annotations do not report long-noncoding RNAs (lncRNAs) (International Aphid Genomics Consortium 2010; Nicholson et al. 2015; Mathers et al. 2017). To gain insight on the potential lncRNA content in aphids, we used RNAseq to predict lncRNAs (see Materials and Methods). A total of 13,478 lncRNAs were predicted in the genome of C. cedri. Importantly, 706 lncRNAs are shared between C. cedri and other aphids. Of these 191 appear to form a conserved core within aphids, and some are conserved across insects (fig. 1). Altogether, given its key phylogenetic position, the C. cedri genome provides an important resource to study genome evolution in aphids.
. 1.

Heatmap of Cinara cedri lncRNA conservation among 20 metazoans. The input lncRNA sequences come from C. cedri. The rows represent C. cedri lncRNAs and the columns represent species. Each cell is colorized based on the level of conservation of the lncRNA (green—100% conservation, black—0% conservation [absence of the lncRNA]).

Heatmap of Cinara cedri lncRNA conservation among 20 metazoans. The input lncRNA sequences come from C. cedri. The rows represent C. cedri lncRNAs and the columns represent species. Each cell is colorized based on the level of conservation of the lncRNA (green—100% conservation, black—0% conservation [absence of the lncRNA]).

Aphid Phylomes and Species-Specific Gene Duplications

As our main focus was to assess gene duplication dynamics in aphids, we reconstructed the complete collection of evolutionary gene histories (i.e., the phylome) of C. cedri, A. pisum, M. persicae, D. noxia, Ap. glycines, and S. flava in the context of other sequenced species (supplementary tables S2 and S3, Supplementary Material online, see Materials and Methods). These genes were scanned to infer duplication and speciation events and derive orthology and paralogy relationships among homologous genes per each phylome (Gabaldón 2008). All of the resulting gene trees, alignments, and orthology and paralogy predictions are available for download or browsing at PhylomeDB (PhylomeIDs: C. cedri—701, S. flava—702, Ap. glycines—703, D. noxia—704, M. persicae—705, A. pisum—706) (Huerta-Cepas et al. 2014). To reconstruct the evolutionary relationships among all considered species, we concatenated the protein alignments of 57 gene trees that are present across all considered species (see Materials and Methods). The resulting highly supported topology (fig. 2) was congruent with current views on aphids phylogeny (Nováková et al. 2013; Chen et al. 2016; Rebijith et al. 2017) and places C. cedri as the earliest branching lineage from our set of aphids.
. 2.

Species tree and duplication ratios of the six phylomes. (a) Phylogenetic tree obtained from the concatenation of 57 widespread gene families. In yellow, all the individuals included in this study that belong to the family Aphididae; in green, light blue, and orange, the aphids that belong to the subfamily Aphidinae, Chaitophorinae, and Lachninae, respectively. All omitted bootstrap values are maximal (bootstrap 100%). (b) Zoom out showing the duplication ratios per each phylome: Cinara cedri—orange, Sipha flava—blue, Aphis glycines—purple, Diuraphis noxia—light blue, Myzus persicae—red, Acyrthosiphon pisum—green. The two branches with the higher duplication ratio are marked as A (ancestral to all six aphids) and B (after the divergence of C. cedri and ancestral to the other five aphids). Bars on the right show the percentage of proteins (orange) associated with transposons in each aphid species. Bemisia tabaci is the outgroup.

Species tree and duplication ratios of the six phylomes. (a) Phylogenetic tree obtained from the concatenation of 57 widespread gene families. In yellow, all the individuals included in this study that belong to the family Aphididae; in green, light blue, and orange, the aphids that belong to the subfamily Aphidinae, Chaitophorinae, and Lachninae, respectively. All omitted bootstrap values are maximal (bootstrap 100%). (b) Zoom out showing the duplication ratios per each phylome: Cinara cedri—orange, Sipha flava—blue, Aphis glycines—purple, Diuraphis noxia—light blue, Myzus persicae—red, Acyrthosiphon pisum—green. The two branches with the higher duplication ratio are marked as A (ancestral to all six aphids) and B (after the divergence of C. cedri and ancestral to the other five aphids). Bars on the right show the percentage of proteins (orange) associated with transposons in each aphid species. Bemisia tabaci is the outgroup. We next focused on gene duplications, including large expansions, that occurred specifically in the lineage leading to each aphid. Interestingly, C. cedri, A. pisum, M. persicae, Ap. glycines, and S. flava have similar proportions of proteins that have an in-paralog (resulting from a species-specific duplication): C. cedri—4,670 (28% of the proteome), S. flava—2,832 (21%), Ap. glycines—3,232 (17%), M. persicae—4,097 (22%), A. pisum—5,431 (29%). These events can be assigned to a similar number of inferred specific gene duplication events: C. cedri—1,420, S. flava—899, Ap. glycines—1,153, M. persicae—1,543, A. pisum—1,889. On the contrary, D. noxia only presented a total of 685 proteins (6% of the proteome) with an in-paralog, corresponding to 315 gene duplication events. In all six aphids, the majority of the gene duplication events result in a moderate number of paralogs (2–5 in-paralogs; supplementary fig. S3, Supplementary Material online), and only few represent large gene family expansions (≥10 in-paralogs). The large expansions could be due to the presence of expanded transposable element families (Huerta-Cepas et al. 2010). In the six aphids, an average of 9% of the total number of annotated protein-coding genes are associated with transposons, with A. pisum and M. persicae (Macrosiphini) containing the highest percentages (fig. 2). Moreover, larger expansions in C. cedri and A. pisum (>50 proteins) often include proteins associated with transposons. However, after removing expansions containing at least one paralog annotated with a PFAM domain or a gene ontology (GO) term associated with transposable elements or viruses, the number of duplications remained high (supplementary fig. S3, Supplementary Material online). We performed a functional GO term enrichment analysis of these transposon-free, species-specific paralogs (table 1) for each proteome. DNA and RNA processing terms were enriched among sets of in-paralogs of all species except D. noxia. Moreover, C. cedri in-paralogs were enriched in GO terms associated with olfactory receptor activity, odorant binding, acetyl-CoA transporter activity, and CCR4-NOT. For S. flava, peroxidase activity, methyltransferase activity, beta-glucosidase activity, lipid droplet, CCR4-NOT complex, and response to oxidative stress were enriched. For Ap. glycines, fatty acid synthase activity, SUMO transferase activity, and regulation of JAK-STAT cascade were enriched. For D. noxia, fucose metabolic process and protein glycosylation were enriched. For M. persicae, peroxidase activity was enriched, and for A. pisum, enoyl-reductase, oleoyl-hydrolase, myristoyl-hydrolase, palmitoyl-hydrolase, odorant binding, and response to stress were enriched. These results are consistent with previous results restricted to A. pisum and M. persicae (Huerta-Cepas et al. 2010; International Aphid Genomics Consortium 2010; Mathers et al. 2017).
Table 1.

List of the GO Terms Enriched in the Expanded Protein Families Specific to Cinara cedri, Sipha flava, Aphis glycines, Diuraphis noxia, Myzus persicae, and Acyrthosiphon pisum.

Term Category TermTerm Level Adj. P-ValueTerm Name
Cinara cedri
molecular_functionGO:000122713.56E-18Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367614.52E-14Nucleic acid binding
molecular_functionGO:000367718.29E-05DNA binding
molecular_functionGO:000369014.57E-13Double-stranded DNA binding
molecular_functionGO:000371515.02E-28Obsolete transcription termination factor activity
molecular_functionGO:000372313.25E-07RNA binding
molecular_functionGO:000372512.56E-07Double-stranded RNA binding
molecular_functionGO:000396419.54E-11RNA-directed DNA polymerase activity
molecular_functionGO:000399415.41E-04Aconitate hydratase activity
molecular_functionGO:000419012.57E-05Aspartic-type endopeptidase activity
molecular_functionGO:000435611.36E-04Glutamate-ammonia ligase activity
molecular_functionGO:000449713.44E-06Monooxygenase activity
molecular_functionGO:000452516.85E-06Ribonuclease III activity
molecular_functionGO:000461813.52E-04Phosphoglycerate kinase activity
molecular_functionGO:000498416.99E-06Olfactory receptor activity
molecular_functionGO:000550616.42E-13Iron ion binding
molecular_functionGO:000554911.04E-04Odorant binding
molecular_functionGO:000852113.14E-05Acetyl-CoA transporter activity
molecular_functionGO:001670511.08E-18Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:001802412.57E-05Histone-lysine N-methyltransferase activity
molecular_functionGO:002003713.00E-05Heme binding
molecular_functionGO:003117713.00E-04Phosphopantetheine binding
molecular_functionGO:003149018.29E-05Chromatin DNA binding
molecular_functionGO:004230212.70E-05Structural constituent of cuticle
cellular_componentGO:000078612.76E-10Nucleosome
cellular_componentGO:003001515.41E-04CCR4-NOT core complex
cellular_componentGO:007087718.29E-05Microprocessor complex
biological_processGO:000627811.60E-09RNA-dependent DNA biosynthetic process
biological_processGO:000635311.79E-26DNA-templated transcription, termination
biological_processGO:000680717.05E-05Nitrogen compound metabolic process
biological_processGO:000945218.38E-057-Methylguanosine RNA capping
biological_processGO:001507414.11E-08DNA integration
biological_processGO:001607511.86E-06rRNA catabolic process
Sipha flava
molecular_functionGO:000016612.39E-08Nucleotide binding
molecular_functionGO:000122716.76E-23Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367612.44E-37Nucleic acid binding
molecular_functionGO:000367811.23E-44DNA helicase activity
molecular_functionGO:000369019.51E-18Double-stranded DNA binding
molecular_functionGO:000369612.71E-22Satellite DNA binding
molecular_functionGO:000369711.53E-11Single-stranded DNA binding
molecular_functionGO:000371511.72E-12Obsolete transcription termination factor activity
molecular_functionGO:000372316.15E-37RNA binding
molecular_functionGO:000373011.15E-19mRNA 3′-UTR binding
molecular_functionGO:000388711.15E-05DNA-directed DNA polymerase activity
molecular_functionGO:000453516.78E-04Poly(A)-specific ribonuclease activity
molecular_functionGO:000460111.99E-11Peroxidase activity
molecular_functionGO:000466618.59E-06Prostaglandin-endoperoxide synthase activity
molecular_functionGO:000480113.80E-04Sedoheptulose-7-phosphate: d-glyceraldehyde-3-phosphate glyceronetransferase activity
molecular_functionGO:000520011.10E-10Structural constituent of cytoskeleton
molecular_functionGO:000816811.26E-25Methyltransferase activity
molecular_functionGO:000840817.74E-043′–5′ Exonuclease activity
molecular_functionGO:001052118.42E-07Telomerase inhibitor activity
molecular_functionGO:002003711.05E-07Heme binding
molecular_functionGO:003294718.59E-06Protein complex scaffold activity
molecular_functionGO:004314112.72E-06ATP-dependent 5′–3′ DNA helicase activity
molecular_functionGO:004316915.84E-06Cation binding
cellular_componentGO:000079211.22E-20Heterochromatin
cellular_componentGO:000565713.02E-04Replication fork
cellular_componentGO:000570113.26E-23Polytene chromosome chromocenter
cellular_componentGO:000581118.41E-04Lipid droplet
cellular_componentGO:000585813.51E-04Axonemal dynein complex
cellular_componentGO:000587417.40E-07Microtubule
cellular_componentGO:003001417.64E-05CCR4-NOT complex
cellular_componentGO:003001511.47E-07CCR4-NOT core complex
cellular_componentGO:003052912.48E-14Intracellular ribonucleoprotein complex
biological_processGO:000000213.02E-04Mitochondrial genome maintenance
biological_processGO:000028818.59E-06Nuclear-transcribed mRNA catabolic process, deadenylation-dependent decay
biological_processGO:000028911.41E-04Nuclear-transcribed mRNA poly(A) tail shortening
biological_processGO:000072311.67E-50Telomere maintenance
biological_processGO:000151014.54E-36RNA methylation
biological_processGO:000597511.29E-05Carbohydrate metabolic process
biological_processGO:000626017.46E-04DNA replication
biological_processGO:000628112.55E-24DNA repair
biological_processGO:000635314.95E-12DNA-templated transcription, termination
biological_processGO:000637011.83E-057-Methylguanosine mRNA capping
biological_processGO:000695413.78E-08Inflammatory response
biological_processGO:000697912.18E-10Response to oxidative stress
biological_processGO:000701711.85E-08Microtubule-based process
biological_processGO:000705914.42E-13Chromosome segregation
biological_processGO:000821718.59E-06Regulation of blood pressure
biological_processGO:000945214.36E-457-Methylguanosine RNA capping
biological_processGO:001607013.65E-13RNA metabolic process
biological_processGO:001937118.59E-06Cyclooxygenase pathway
biological_processGO:003026117.60E-17Chromosome condensation
biological_processGO:003150713.26E-23Heterochromatin assembly
biological_processGO:003221118.42E-07Negative regulation of telomere maintenance via telomerase
biological_processGO:003225911.86E-04Methylation
biological_processGO:004480618.42E-07G-quadruplex DNA unwinding
biological_processGO:004572714.44E-20Positive regulation of translation
biological_processGO:005125813.88E-05Protein polymerization
biological_processGO:005197418.42E-07Negative regulation of telomerase activity
biological_processGO:190165715.82E-05Glycosyl compound metabolic process
Aphis glycines
molecular_functionGO:000122712.36E-04Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367613.08E-37Nucleic acid binding
molecular_functionGO:000367719.92E-05DNA binding
molecular_functionGO:000371514.51E-19Obsolete transcription termination factor activity
molecular_functionGO:000388711.12E-04DNA-directed DNA polymerase activity
molecular_functionGO:000431215.19E-08Fatty acid synthase activity
molecular_functionGO:000455318.62E-05Hydrolase activity, hydrolyzing O-glycosyl compounds
molecular_functionGO:000486615.12E-04Endopeptidase inhibitor activity
molecular_functionGO:000486911.04E-05Cysteine-type endopeptidase inhibitor activity
molecular_functionGO:000840815.21E-053′–5′ Exonuclease activity
molecular_functionGO:000842211.12E-04Beta-glucosidase activity
molecular_functionGO:000852115.75E-08Acetyl-CoA transporter activity
molecular_functionGO:001978911.11E-08SUMO transferase activity
molecular_functionGO:003117714.71E-05Phosphopantetheine binding
molecular_functionGO:004302711.55E-05Cysteine-type endopeptidase inhibitor activity involved in apoptotic process
molecular_functionGO:004316911.21E-07Cation binding
molecular_functionGO:004439015.12E-04Ubiquitin-like protein conjugating enzyme binding
molecular_functionGO:006166311.50E-04NEDD8 ligase activity
molecular_functionGO:008972011.55E-05Caspase binding
cellular_componentGO:000565212.75E-04Nuclear lamina
cellular_componentGO:000570514.07E-04Polytene chromosome interband
cellular_componentGO:000587611.91E-04Spindle microtubule
cellular_componentGO:000853711.55E-05Proteasome activator complex
cellular_componentGO:003501212.18E-05Polytene chromosome, telomeric region
cellular_componentGO:007077615.82E-08MOZ/MORF histone acetyltransferase complex
biological_processGO:000151017.60E-13RNA methylation
biological_processGO:000597511.18E-08Carbohydrate metabolic process
biological_processGO:000635311.61E-18DNA-templated transcription, termination
biological_processGO:000728915.12E-04Spermatid nucleus differentiation
biological_processGO:000744613.21E-05Imaginal disc growth
biological_processGO:000945216.31E-187-Methylguanosine RNA capping
biological_processGO:003026119.74E-06Chromosome condensation
biological_processGO:004315411.13E-04Negative regulation of cysteine-type endopeptidase activity involved in apoptotic process
biological_processGO:004642516.74E-04Regulation of JAK-STAT cascade
biological_processGO:004642614.76E-04Negative regulation of JAK-STAT cascade
biological_processGO:007093614.71E-05Protein K48-linked ubiquitination
biological_processGO:009030716.55E-04Mitotic spindle assembly
biological_processGO:009734011.50E-04Inhibition of cysteine-type endopeptidase activity
biological_processGO:190165712.18E-05Glycosyl compound metabolic process
biological_processGO:199000112.82E-07Inhibition of cysteine-type endopeptidase activity involved in apoptotic process
biological_processGO:200127111.50E-04Negative regulation of cysteine-type endopeptidase activity involved in execution phase of apoptosis
Diuraphis noxia
molecular_functionGO:000367816.51E-05DNA helicase activity
molecular_functionGO:000482712.22E-05Proline-tRNA ligase activity
molecular_functionGO:000842411.34E-04Glycoprotein 6-alpha-l-fucosyltransferase activity
molecular_functionGO:004692116.24E-05Alpha-(1->6)-fucosyltransferase activity
biological_processGO:000072315.82E-06Telomere maintenance
biological_processGO:000643312.22E-05Prolyl-tRNA aminoacylation
biological_processGO:003357811.34E-04Protein glycosylation in Golgi
biological_processGO:003607111.34E-04 N-glycan fucosylation
biological_processGO:004636816.24E-05GDP-l-fucose metabolic process
Myzus persicae
molecular_functionGO:000122711.86E-13Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367611.79E-25Nucleic acid binding
molecular_functionGO:000367712.30E-11DNA binding
molecular_functionGO:000369011.17E-09Double-stranded DNA binding
molecular_functionGO:000377715.98E-04Microtubule motor activity
molecular_functionGO:000414912.89E-04Dihydrolipoyllysine-residue succinyltransferase activity
molecular_functionGO:000460111.38E-04Peroxidase activity
molecular_functionGO:000481816.76E-05Glutamate-tRNA ligase activity
molecular_functionGO:000482711.36E-04Proline-tRNA ligase activity
molecular_functionGO:002003715.22E-06Heme binding
biological_processGO:000642412.02E-05Glutamyl-tRNA aminoacylation
biological_processGO:000643311.36E-04Prolyl-tRNA aminoacylation
biological_processGO:000659613.16E-06Polyamine biosynthetic process
biological_processGO:001692512.30E-07Protein sumoylation
Acyrthosiphon pisum
molecular_functionGO:000367614.01E-14Nucleic acid binding
molecular_functionGO:000367812.18E-11DNA helicase activity
molecular_functionGO:000371514.88E-05Obsolete transcription termination factor activity
molecular_functionGO:000417718.98E-13Aminopeptidase activity
molecular_functionGO:000425213.43E-19Serine-type endopeptidase activity
molecular_functionGO:000431317.14E-05[acyl-carrier-protein] S-acetyltransferase activity
molecular_functionGO:000431717.14E-053-Hydroxypalmitoyl-[acyl-carrier-protein] dehydratase activity
molecular_functionGO:000431917.14E-05Enoyl-[acyl-carrier-protein] reductase (NADPH, B-specific) activity
molecular_functionGO:000432017.14E-05Oleoyl-[acyl-carrier-protein] hydrolase activity
molecular_functionGO:000460115.92E-23Peroxidase activity
molecular_functionGO:000474817.45E-05Ribonucleoside-diphosphate reductase activity, thioredoxin disulfide as acceptor
molecular_functionGO:000488812.47E-05Transmembrane signaling receptor activity
molecular_functionGO:000550616.66E-07Iron ion binding
molecular_functionGO:000554914.97E-04Odorant binding
molecular_functionGO:000823419.62E-10Cysteine-type peptidase activity
molecular_functionGO:000823711.20E-10Metallopeptidase activity
molecular_functionGO:001629517.14E-05Myristoyl-[acyl-carrier-protein] hydrolase activity
molecular_functionGO:001629617.14E-05Palmitoyl-[acyl-carrier-protein] hydrolase activity
molecular_functionGO:001670511.03E-06Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:002003711.23E-23Heme binding
molecular_functionGO:003017018.55E-04Pyridoxal phosphate binding
molecular_functionGO:004316915.10E-04Cation binding
cellular_componentGO:000558111.03E-06Collagen trimer
cellular_componentGO:003501217.45E-05Polytene chromosome, telomeric region
cellular_componentGO:004260019.51E-11Chorion
biological_processGO:000072312.18E-11Telomere maintenance
biological_processGO:000597511.46E-06Carbohydrate metabolic process
biological_processGO:000626014.02E-04DNA replication
biological_processGO:000635312.29E-05DNA-templated transcription, termination
biological_processGO:000650811.12E-33Proteolysis
biological_processGO:000685711.00E-04Oligopeptide transport
biological_processGO:000697911.27E-14Response to oxidative stress
biological_processGO:000716611.04E-04Cell surface receptor signaling pathway
biological_processGO:000926316.11E-04Deoxyribonucleotide biosynthetic process
biological_processGO:003519416.11E-04Posttranscriptional gene silencing by RNA
List of the GO Terms Enriched in the Expanded Protein Families Specific to Cinara cedri, Sipha flava, Aphis glycines, Diuraphis noxia, Myzus persicae, and Acyrthosiphon pisum. In order to detect parallel duplications, we searched for orthologs between C. cedri and the other aphids with species-specific duplications. A total of 909 C. cedri genes (26% of the total proteins with in-paralogs) with species-specific duplications have parallel species-specific duplications in at least one of the other aphids (S. flava—364, Ap. glycines—181, D. noxia—57, M. persicae—235, A. pisum—375). Specifically, 694 C. cedri genes share unique parallel duplications with one of the other aphids: S. flava—252, Ap. glycines—91, D. noxia—14, M. persicae—120, A. pisum—217. Interestingly, C. cedri parallel paralogs show enrichments only in four aphids. The parallel duplications shared with S. flava show enrichment for aconitate hydratase activity, l-amino acid transmembrane transporter activity, tricarboxylic acid cycle, aromatase activity, and CCR4-NOT complex. Acyrthosiphon pisum species-specific duplications shared with C. cedri show enrichment for oxidoreductase activity and l-ascorbic acid binding. Cinara cedri duplications shared with Ap. glycines and M. persicae show only five and three enriched terms, respectively (table 2). Interestingly, two proteins show parallel duplications in all the considered aphid species, from which only one has a functional annotation. This protein is associated with UDP-N-acetylglucosamine–peptide N-acetylglucosaminyltransferase 110 kDa subunit-like, which catalyzes the transfer of a single N-acetylglucosamine from UDP-GlcNAc to a serine or threonine residue (O-GlcNAc glycosylation) (Lazarus et al. 2012; Ding et al. 2015). In insects, this type of glycosylation has been shown to be central to a variety of physiological processes, including regulation of the cell cycle, expression of developmental genes, nutrient sensing, response to starvation, insulin signaling, or specification of body size (Vandenborre et al. 2011; Walski et al. 2017). Altogether, these results indicate a high dynamism of aphid gene repertoire and suggest that gene duplication may play a major role in the adaptation of aphid species to their respective environments.
Table 2.

List of GO Terms Enriched in the Parallel Cinara cedri Species-Specific Duplications Uniquely Shared with One of the Other Aphids (Sipha flava, Aphis glycines, Myzus persicae, and Acyrthosiphon pisum).

Term Category TermTerm LevelAdj. P-ValueTerm Name
Sipha flava
molecular_functionGO:000399411.93E-13Aconitate hydratase activity
molecular_functionGO:000449716.19E-11Monooxygenase activity
molecular_functionGO:000453511.84E-04Poly(A)-specific ribonuclease activity
molecular_functionGO:000520011.34E-08Structural constituent of cytoskeleton
molecular_functionGO:000550611.23E-23Iron ion binding
molecular_functionGO:001517913.88E-04 l-amino acid transmembrane transporter activity
molecular_functionGO:001670511.36E-27Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:002003711.46E-22Heme binding
molecular_functionGO:003294713.88E-04Protein complex scaffold activity
molecular_functionGO:005153911.23E-054 Iron, 4 sulfur cluster binding
molecular_functionGO:007033011.14E-04Aromatase activity
cellular_componentGO:000587413.19E-05Microtubule
cellular_componentGO:003001414.67E-05CCR4-NOT complex
cellular_componentGO:003001514.91E-05CCR4-NOT core complex
biological_processGO:000028813.88E-04Nuclear-transcribed mRNA catabolic process, deadenylation-dependent decay
biological_processGO:000028913.41E-05Nuclear-transcribed mRNA poly(A) tail shortening
biological_processGO:000609916.27E-06Tricarboxylic acid cycle
biological_processGO:000640214.17E-04mRNA catabolic process
biological_processGO:000701718.35E-08Microtubule-based process
biological_processGO:001714811.46E-04Negative regulation of translation
biological_processGO:005508519.30E-04Transmembrane transport
Acyrthosiphon pisum
molecular_functionGO:000402017.10E-05Adenylylsulfate kinase activity
molecular_functionGO:000419016.34E-04Aspartic-type endopeptidase activity
molecular_functionGO:000465611.64E-09Procollagen-proline 4-dioxygenase activity
molecular_functionGO:000471911.57E-04Protein-l-isoaspartate (d-aspartate) O-methyltransferase activity
molecular_functionGO:001670213.02E-06Oxidoreductase activity, acting on single donors with incorporation of molecular oxygen, incorporation of two atoms of oxygen
molecular_functionGO:003141811.99E-07 l-ascorbic acid binding
biological_processGO:000010319.52E-05Sulfate assimilation
Aphis glycines
molecular_functionGO:000367611.43E-07Nucleic acid binding
molecular_functionGO:000392418.18E-04GTPase activity
cellular_componentGO:000574113.64E-05Mitochondrial outer membrane
biological_processGO:000805316.78E-13Mitochondrial fusion
biological_processGO:004866211.99E-04Negative regulation of smooth muscle cell proliferation
Myzus persicae
molecular_functionGO:000377713.06E-06Microtubule motor activity
molecular_functionGO:000396413.91E-04RNA-directed DNA polymerase activity
biological_processGO:000701812.38E-04Microtubule-based movement
List of GO Terms Enriched in the Parallel Cinara cedri Species-Specific Duplications Uniquely Shared with One of the Other Aphids (Sipha flava, Aphis glycines, Myzus persicae, and Acyrthosiphon pisum).

High Number of Ancient Gene Duplications Suggests One Ancestral Burst of Large-Scale Genome Duplication

In order to detect waves of ancestral duplications in the evolutionary history of aphids, we used a phylogeny-based phylostratigraphic approach based on a species-overlap algorithm (Huerta-Cepas and Gabaldón 2011) to detect gene duplications and map them onto the species tree (see Materials and Methods). After excluding large expansions (duplications resulting in >5 paralogs), we computed ratios of gene duplications (average number of duplications per gene detected in a given branch of the species tree) for each phylome (fig. 2). Interestingly, in the aphid lineage two branches have high duplication ratios: one present in the ancestral branch of all six aphids (Aphididae family, branch A), probably related to their adaptation to specific diet and life-style, and the other after the divergence of C. cedri and ancestral to the other five aphids (branch B), which could be related to a nonconifer-specific diet. To validate these findings, we analyzed the relative age of the duplications by plotting the ratio of transversions at 4-fold degenerate sites (4DTv) of paralogs mapped at the two branches with high duplication ratios (supplementary fig. S4, Supplementary Material online). We also mapped two speciation events per each phylome by plotting the 4DTv of orthologous gene pairs (see Materials and Methods). Unexpectedly, the distribution of the 4DTv of both waves of duplications detected by phylostratigraphy was fully overlapping within the period of time corresponding to the most ancestral duplication. This suggests that the most recent peak of duplications detected through topological analyses of gene trees may result from more ancestral duplications followed by loss of both paralogs in C. cedri, therefore rendering a topology, that indicates a more recent duplication event. Consistent with this interpretation, 70% of the genes duplicated in branch B do not have an ortholog in C. cedri. From these observations, we conclude that the second apparent duplication peak is actually the result of differential retention of duplicates and our limited sampling of early-branching lineages. Such large levels of differential retention of duplicates have also been observed in other organisms, such as Paramecium (McGrath et al. 2014) and Brassica (Mun et al. 2009). To test for the robustness of the detected ancestral wave of duplications, we applied stronger filters, by considering gene trees containing ancestral aphid duplications with a maximum of five genes per aphid. In all five phylomes, an average of 76% of gene trees passed this filter: C. cedri—11,304 (78%), S. flava—9,781 (75%), Ap. glycines—13,170 (80%), D. noxia—10,379 (85%), M. persicae—11,759 (70%), and A. pisum—12,485 (71%). When the duplication ratios were calculated using this more restricted set of gene trees, only the duplication ratio at the ancestral branch of all six aphids was still apparent (see Materials and Methods; supplementary fig. S5, Supplementary Material online). Taken together, these results suggest that there was one large-scale genome duplication in the evolutionary history of aphids predating the divergence of the Aphididae family, which could be related to adaptive innovations. A functional enrichment analysis of the proteins duplicated in the ancestral branch of the six aphids showed enrichment for annotations related to carbohydrate metabolic process, response to stimulus, olfactory receptor activity, odorant binding, glucuronidation, transmembrane transporter activity, and DNA and RNA processing, among others (table 3).
Table 3.

List of the GO Terms Enriched in the Duplicated Protein Families at the Base of All Six Aphids per Each Phylome.

Term CategoryTermTerm LevelAdj. P-ValueTerm Name
Cinara cedri
molecular_functionGO:000122711.35E-05Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367611.59E-30Nucleic acid binding
molecular_functionGO:000396415.96E-28RNA-directed DNA polymerase activity
molecular_functionGO:000403218.29E-05Alditol:NADP+ 1-oxidoreductase activity
molecular_functionGO:000418513.69E-04Serine-type carboxypeptidase activity
molecular_functionGO:000419711.23E-09Cysteine-type endopeptidase activity
molecular_functionGO:000431614.99E-053-Oxoacyl-[acyl-carrier-protein] reductase (NADPH) activity
molecular_functionGO:000449717.28E-05Monooxygenase activity
molecular_functionGO:000452313.51E-06RNA–DNA hybrid ribonuclease activity
molecular_functionGO:000455317.45E-11Hydrolase activity, hydrolyzing O-glycosyl compounds
molecular_functionGO:000455511.14E-05Alpha, alpha-trehalase activity
molecular_functionGO:000498411.03E-12Olfactory receptor activity
molecular_functionGO:000521515.81E-08Transporter activity
molecular_functionGO:000525416.47E-05Chloride channel activity
molecular_functionGO:000535516.09E-04Glucose transmembrane transporter activity
molecular_functionGO:000550612.92E-09Iron ion binding
molecular_functionGO:000554211.12E-04Folic acid binding
molecular_functionGO:000554913.15E-07Odorant binding
molecular_functionGO:000819413.76E-04UDP-glycosyltransferase activity
molecular_functionGO:000823411.18E-07Cysteine-type peptidase activity
molecular_functionGO:000841711.12E-04Fucosyltransferase activity
molecular_functionGO:000842219.39E-04Beta-glucosidase activity
molecular_functionGO:000851811.12E-04Reduced folate carrier activity
molecular_functionGO:000852113.90E-17Acetyl-CoA transporter activity
molecular_functionGO:001502011.96E-13Glucuronosyltransferase activity
molecular_functionGO:001517111.55E-04Amino acid transmembrane transporter activity
molecular_functionGO:001529511.03E-06Solute:proton symporter activity
molecular_functionGO:001529719.39E-04Antiporter activity
molecular_functionGO:001529919.39E-04Solute:proton antiporter activity
molecular_functionGO:001552813.51E-06Lactose:proton symporter activity
molecular_functionGO:001670511.25E-13Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:001675811.23E-20Transferase activity, transferring hexosyl groups
molecular_functionGO:001678819.89E-04Hydrolase activity, acting on ester bonds
molecular_functionGO:001711013.60E-04Nucleoside-diphosphatase activity
molecular_functionGO:001802414.89E-06Histone-lysine N-methyltransferase activity
molecular_functionGO:001979913.70E-05Tubulin N-acetyltransferase activity
molecular_functionGO:002003712.79E-11Heme binding
molecular_functionGO:002285717.61E-25Transmembrane transporter activity
molecular_functionGO:002289119.33E-06Substrate-specific transmembrane transporter activity
molecular_functionGO:003149013.60E-04Chromatin DNA binding
molecular_functionGO:003519712.40E-04siRNA binding
molecular_functionGO:004316913.94E-06Cation binding
molecular_functionGO:009048211.12E-04Vitamin transmembrane transporter activity
molecular_functionGO:010233619.83E-043-Oxo-arachidoyl-CoA synthase activity
molecular_functionGO:010233719.83E-043-Oxo-cerotoyl-CoA synthase activity
molecular_functionGO:010233819.83E-043-Oxo-lignoceronyl-CoA synthase activity
cellular_componentGO:000853711.12E-04Proteasome activator complex
biological_processGO:000597516.65E-05Carbohydrate metabolic process
biological_processGO:000599111.14E-05Trehalose metabolic process
biological_processGO:000627816.80E-24RNA-dependent DNA biosynthetic process
biological_processGO:000631011.25E-05DNA recombination
biological_processGO:000650812.74E-10Proteolysis
biological_processGO:000662917.28E-05Lipid metabolic process
biological_processGO:000682018.83E-05Anion transport
biological_processGO:000728312.35E-05Spermatogenesis
biological_processGO:000760818.86E-07Sensory perception of smell
biological_processGO:001507411.08E-16DNA integration
biological_processGO:001697319.39E-04Poly(A)+ mRNA export from nucleus
biological_processGO:003016219.39E-04Regulation of proteolysis
biological_processGO:003542813.03E-04Hexose transmembrane transport
biological_processGO:004632313.03E-04Glucose import
biological_processGO:005079012.49E-07Regulation of catalytic activity
biological_processGO:005118011.12E-04Vitamin transport
biological_processGO:005160314.08E-04Proteolysis involved in cellular protein catabolic process
biological_processGO:005269611.12E-04Flavonoid glucuronidation
biological_processGO:005508519.23E-34Transmembrane transport
biological_processGO:005511412.21E-04Oxidation–reduction process
biological_processGO:007192913.70E-05Alpha-tubulin acetylation
biological_processGO:190165713.70E-05Glycosyl compound metabolic process
Sipha flava
molecular_functionGO:000122711.78E-45Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367612.51E-19Nucleic acid binding
molecular_functionGO:000367712.55E-06DNA binding
molecular_functionGO:000367815.62E-15DNA helicase activity
molecular_functionGO:000369011.36E-30Double-stranded DNA binding
molecular_functionGO:000371514.74E-16Obsolete transcription termination factor activity
molecular_functionGO:000396417.50E-07RNA-directed DNA polymerase activity
molecular_functionGO:000431617.50E-073-Oxoacyl-[acyl-carrier-protein] reductase (NADPH) activity
molecular_functionGO:000449713.17E-04Monooxygenase activity
molecular_functionGO:000455311.08E-08Hydrolase activity, hydrolyzing O-glycosyl compounds
molecular_functionGO:000455512.62E-06Alpha, alpha-trehalase activity
molecular_functionGO:000456511.64E-04Beta-galactosidase activity
molecular_functionGO:000498411.45E-14Olfactory receptor activity
molecular_functionGO:000521518.75E-09Transporter activity
molecular_functionGO:000535112.79E-06Sugar:proton symporter activity
molecular_functionGO:000535512.59E-07Glucose transmembrane transporter activity
molecular_functionGO:000550612.01E-05Iron ion binding
molecular_functionGO:000554911.14E-11Odorant binding
molecular_functionGO:000819414.71E-07UDP-glycosyltransferase activity
molecular_functionGO:000852111.51E-10Acetyl-CoA transporter activity
molecular_functionGO:001052117.20E-04Telomerase inhibitor activity
molecular_functionGO:001502014.52E-10Glucuronosyltransferase activity
molecular_functionGO:001661419.74E-05Oxidoreductase activity, acting on CH–OH group of donors
molecular_functionGO:001670519.69E-08Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:001675818.93E-18Transferase activity, transferring hexosyl groups
molecular_functionGO:002003714.91E-12Heme binding
molecular_functionGO:002285713.25E-21Transmembrane transporter activity
molecular_functionGO:002289111.02E-06Substrate-specific transmembrane transporter activity
molecular_functionGO:004230213.45E-07Structural constituent of cuticle
molecular_functionGO:004316914.58E-06Cation binding
cellular_componentGO:004323113.60E-04Intracellular membrane-bounded organelle
biological_processGO:000072311.59E-18Telomere maintenance
biological_processGO:000599112.62E-06Trehalose metabolic process
biological_processGO:000628111.66E-05DNA repair
biological_processGO:000635311.20E-16DNA-templated transcription, termination
biological_processGO:000650816.34E-06Proteolysis
biological_processGO:000685718.42E-05Oligopeptide transport
biological_processGO:000760814.80E-06Sensory perception of smell
biological_processGO:003221117.20E-04Negative regulation of telomere maintenance via telomerase
biological_processGO:003542812.59E-07Hexose transmembrane transport
biological_processGO:004480617.20E-04G-quadruplex DNA unwinding
biological_processGO:004632312.59E-07Glucose import
biological_processGO:005090914.30E-05Sensory perception of taste
biological_processGO:005197417.20E-04Negative regulation of telomerase activity
biological_processGO:005269614.87E-08Flavonoid glucuronidation
biological_processGO:005508516.34E-26Transmembrane transport
Aphis glycines
molecular_functionGO:000122711.41E-09Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367615.84E-43Nucleic acid binding
molecular_functionGO:000369013.51E-05Double-stranded DNA binding
molecular_functionGO:000371513.33E-20Obsolete transcription termination factor activity
molecular_functionGO:000396411.45E-05RNA-directed DNA polymerase activity
molecular_functionGO:000419711.97E-06Cysteine-type endopeptidase activity
molecular_functionGO:000449719.76E-06Monooxygenase activity
molecular_functionGO:000455317.60E-08Hydrolase activity, hydrolyzing O-glycosyl compounds
molecular_functionGO:000498412.42E-20Olfactory receptor activity
molecular_functionGO:000535111.26E-06Sugar:proton symporter activity
molecular_functionGO:000535517.12E-08Glucose transmembrane transporter activity
molecular_functionGO:000550612.39E-07Iron ion binding
molecular_functionGO:000554915.98E-17Odorant binding
molecular_functionGO:000819411.17E-06UDP-glycosyltransferase activity
molecular_functionGO:000823416.93E-04Cysteine-type peptidase activity
molecular_functionGO:000852112.76E-21Acetyl-CoA transporter activity
molecular_functionGO:001502017.06E-13Glucuronosyltransferase activity
molecular_functionGO:001529513.29E-07Solute:proton symporter activity
molecular_functionGO:001670517.31E-10Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:001675814.95E-25Transferase activity, transferring hexosyl groups
molecular_functionGO:001978912.14E-08SUMO transferase activity
molecular_functionGO:001979919.15E-05Tubulin N-acetyltransferase activity
molecular_functionGO:002003712.22E-15Heme binding
molecular_functionGO:002285717.22E-14Transmembrane transporter activity
molecular_functionGO:002289119.76E-06Substrate-specific transmembrane transporter activity
molecular_functionGO:004262611.89E-05ATPase activity, coupled to transmembrane movement of substances
molecular_functionGO:004316911.95E-04Cation binding
molecular_functionGO:005066013.55E-04Flavin adenine dinucleotide binding
molecular_functionGO:009048216.11E-04Vitamin transmembrane transporter activity
cellular_componentGO:000853715.67E-06Proteasome activator complex
cellular_componentGO:003501211.47E-05Polytene chromosome, telomeric region
biological_processGO:000597515.43E-06Carbohydrate metabolic process
biological_processGO:000635313.25E-18DNA-templated transcription, termination
biological_processGO:000650811.08E-05Proteolysis
biological_processGO:000662918.02E-04Lipid metabolic process
biological_processGO:000709511.30E-04Mitotic G2 DNA damage checkpoint
biological_processGO:000744612.35E-06Imaginal disc growth
biological_processGO:000760811.34E-10Sensory perception of smell
biological_processGO:003009719.79E-04Hemopoiesis
biological_processGO:003296816.93E-04Positive regulation of transcription elongation from RNA polymerase II promoter
biological_processGO:003542817.12E-08Hexose transmembrane transport
biological_processGO:004217616.93E-04Regulation of protein catabolic process
biological_processGO:004632317.12E-08Glucose import
biological_processGO:004642511.27E-06Regulation of JAK-STAT cascade
biological_processGO:005079011.27E-06Regulation of catalytic activity
biological_processGO:005090916.62E-04Sensory perception of taste
biological_processGO:005118016.11E-04Vitamin transport
biological_processGO:005269612.03E-07Flavonoid glucuronidation
biological_processGO:005508512.74E-19Transmembrane transport
biological_processGO:007192919.15E-05Alpha-tubulin acetylation
biological_processGO:190165716.11E-04Glycosyl compound metabolic process
Diuraphis noxia
molecular_functionGO:000006415.19E-04 L-ornithine transmembrane transporter activity
molecular_functionGO:000367618.65E-17Nucleic acid binding
molecular_functionGO:000367813.25E-04DNA helicase activity
molecular_functionGO:000371515.79E-09Obsolete transcription termination factor activity
molecular_functionGO:000396418.84E-05RNA-directed DNA polymerase activity
molecular_functionGO:000419714.18E-06Cysteine-type endopeptidase activity
molecular_functionGO:000431613.19E-043-Oxoacyl-[acyl-carrier-protein] reductase (NADPH) activity
molecular_functionGO:000439613.28E-04Hexokinase activity
molecular_functionGO:000449711.51E-10Monooxygenase activity
molecular_functionGO:000455311.67E-07Hydrolase activity, hydrolyzing O-glycosyl compounds
molecular_functionGO:000455513.76E-09Alpha, alpha-trehalase activity
molecular_functionGO:000460112.82E-05Peroxidase activity
molecular_functionGO:000498411.80E-04Olfactory receptor activity
molecular_functionGO:000521513.17E-07Transporter activity
molecular_functionGO:000535112.13E-10Sugar:proton symporter activity
molecular_functionGO:000535511.01E-11Glucose transmembrane transporter activity
molecular_functionGO:000550611.93E-09Iron ion binding
molecular_functionGO:000553613.28E-04Glucose binding
molecular_functionGO:000554915.44E-04Odorant binding
molecular_functionGO:000819411.78E-07UDP-glycosyltransferase activity
molecular_functionGO:000823415.08E-05Cysteine-type peptidase activity
molecular_functionGO:000842413.28E-04Glycoprotein 6-alpha-L-fucosyltransferase activity
molecular_functionGO:000852113.96E-22Acetyl-CoA transporter activity
molecular_functionGO:001502015.59E-13Glucuronosyltransferase activity
molecular_functionGO:001517418.84E-05Basic amino acid transmembrane transporter activity
molecular_functionGO:001518115.19E-04Arginine transmembrane transporter activity
molecular_functionGO:001518915.19E-04 l-lysine transmembrane transporter activity
molecular_functionGO:001529512.27E-05Solute:proton symporter activity
molecular_functionGO:001532618.84E-05Basic amino acid transmembrane transporter activity
molecular_functionGO:001661413.76E-09Oxidoreductase activity, acting on CH–OH group of donors
molecular_functionGO:001670512.96E-16Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:001675811.58E-24Transferase activity, transferring hexosyl groups
molecular_functionGO:001678813.00E-04Hydrolase activity, acting on ester bonds
molecular_functionGO:001687218.84E-05Intramolecular lyase activity
molecular_functionGO:001979912.27E-05Tubulin N-acetyltransferase activity
molecular_functionGO:002003711.52E-22Heme binding
molecular_functionGO:002285713.73E-21Transmembrane transporter activity
molecular_functionGO:002289113.09E-08Substrate-specific transmembrane transporter activity
molecular_functionGO:004262611.91E-07ATPase activity, coupled to transmembrane movement of substances
molecular_functionGO:004316919.84E-06Cation binding
molecular_functionGO:005066012.64E-06Flavin adenine dinucleotide binding
cellular_componentGO:000587417.95E-04Microtubule
cellular_componentGO:001602111.56E-06Integral component of membrane
cellular_componentGO:003146113.19E-04Cullin-RING ubiquitin ligase complex
cellular_componentGO:004323111.28E-05Intracellular membrane-bounded organelle
biological_processGO:000072312.83E-04Telomere maintenance
biological_processGO:000167813.28E-04Cellular glucose homeostasis
biological_processGO:000597519.76E-06Carbohydrate metabolic process
biological_processGO:000599113.76E-09Trehalose metabolic process
biological_processGO:000599313.28E-04Trehalose catabolic process
biological_processGO:000609618.93E-05Glycolytic process
biological_processGO:000635216.52E-05DNA-templated transcription, initiation
biological_processGO:000635315.79E-09DNA-templated transcription, termination
biological_processGO:000650813.17E-08Proteolysis
biological_processGO:000662911.16E-06Lipid metabolic process
biological_processGO:000686518.84E-05Amino acid transport
biological_processGO:000735213.28E-04Zygotic specification of dorsal/ventral axis
biological_processGO:000945213.28E-047-Methylguanosine RNA capping
biological_processGO:003357813.28E-04Protein glycosylation in Golgi
biological_processGO:003542811.01E-11Hexose transmembrane transport
biological_processGO:004632311.01E-11Glucose import
biological_processGO:005079015.25E-04Regulation of catalytic activity
biological_processGO:005090914.68E-04Sensory perception of taste
biological_processGO:005269611.60E-08Flavonoid glucuronidation
biological_processGO:005508516.40E-25Transmembrane transport
biological_processGO:005511411.71E-04Oxidation-reduction process
biological_processGO:007192912.27E-05Alpha-tubulin acetylation
biological_processGO:190335215.19E-04 L-ornithine transmembrane transport
Myzus persicae
molecular_functionGO:000122713.22E-16Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367613.81E-65Nucleic acid binding
molecular_functionGO:000367716.51E-16DNA binding
molecular_functionGO:000369011.01E-08Double-stranded DNA binding
molecular_functionGO:000395015.43E-07NAD+ ADP-ribosyltransferase activity
molecular_functionGO:000418513.40E-04Serine-type carboxypeptidase activity
molecular_functionGO:000419713.56E-09Cysteine-type endopeptidase activity
molecular_functionGO:000425211.33E-05Serine-type endopeptidase activity
molecular_functionGO:000431611.38E-043-Oxoacyl-[acyl-carrier-protein] reductase (NADPH) activity
molecular_functionGO:000439611.87E-05Hexokinase activity
molecular_functionGO:000449712.72E-07Monooxygenase activity
molecular_functionGO:000455311.59E-07Hydrolase activity, hydrolyzing O-glycosyl compounds
molecular_functionGO:000455514.59E-09Alpha, alpha-trehalase activity
molecular_functionGO:000456511.86E-04Beta-galactosidase activity
molecular_functionGO:000460117.64E-06Peroxidase activity
molecular_functionGO:000498411.49E-21Olfactory receptor activity
molecular_functionGO:000521512.40E-04Transporter activity
molecular_functionGO:000535118.69E-06Sugar:proton symporter activity
molecular_functionGO:000535511.92E-07Glucose transmembrane transporter activity
molecular_functionGO:000550617.01E-10Iron ion binding
molecular_functionGO:000553611.87E-05Glucose binding
molecular_functionGO:000554919.69E-19Odorant binding
molecular_functionGO:000819411.02E-05UDP-glycosyltransferase activity
molecular_functionGO:000823413.73E-08Cysteine-type peptidase activity
molecular_functionGO:000823613.71E-04Serine-type peptidase activity
molecular_functionGO:000852112.99E-20Acetyl-CoA transporter activity
molecular_functionGO:001502018.81E-14Glucuronosyltransferase activity
molecular_functionGO:001517411.86E-04Basic amino acid transmembrane transporter activity
molecular_functionGO:001529516.02E-05Solute:proton symporter activity
molecular_functionGO:001532611.86E-04Basic amino acid transmembrane transporter activity
molecular_functionGO:001670512.44E-17Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:001675813.11E-27Transferase activity, transferring hexosyl groups
molecular_functionGO:001676515.72E-05Transferase activity, transferring alkyl or aryl (other than methyl) groups
molecular_functionGO:001687211.59E-07Intramolecular lyase activity
molecular_functionGO:001979911.87E-05Tubulin N-acetyltransferase activity
molecular_functionGO:002003714.37E-20Heme binding
molecular_functionGO:002285711.36E-16Transmembrane transporter activity
molecular_functionGO:002289111.30E-04Substrate-specific transmembrane transporter activity
molecular_functionGO:003117713.23E-04Phosphopantetheine binding
molecular_functionGO:003377715.89E-04Lithocholate 6beta-hydroxylase activity
molecular_functionGO:004230215.89E-04Structural constituent of cuticle
molecular_functionGO:004316911.28E-04Cation binding
cellular_componentGO:003259012.21E-04Dendrite membrane
biological_processGO:000167811.87E-05Cellular glucose homeostasis
biological_processGO:000597512.35E-05Carbohydrate metabolic process
biological_processGO:000599114.59E-09Trehalose metabolic process
biological_processGO:000599311.86E-04Trehalose catabolic process
biological_processGO:000650815.25E-17Proteolysis
biological_processGO:000685715.03E-07Oligopeptide transport
biological_processGO:000686511.86E-04Amino acid transport
biological_processGO:000760816.70E-12Sensory perception of smell
biological_processGO:003542811.92E-07Hexose transmembrane transport
biological_processGO:004632311.92E-07Glucose import
biological_processGO:005079012.38E-08Regulation of catalytic activity
biological_processGO:005089617.80E-04Response to stimulus
biological_processGO:005090916.70E-12Sensory perception of taste
biological_processGO:005160315.47E-05Proteolysis involved in cellular protein catabolic process
biological_processGO:005269611.59E-07Flavonoid glucuronidation
biological_processGO:005508511.33E-13Transmembrane transport
biological_processGO:007192911.87E-05Alpha-tubulin acetylation
Acyrthosiphon pisum
molecular_functionGO:000006415.99E-05 l-ornithine transmembrane transporter activity
molecular_functionGO:000122717.63E-08Transcriptional repressor activity, RNA polymerase II transcription regulatory region sequence-specific binding
molecular_functionGO:000367617.77E-87Nucleic acid binding
molecular_functionGO:000367811.82E-10DNA helicase activity
molecular_functionGO:000371513.16E-08Obsolete transcription termination factor activity
molecular_functionGO:000374319.37E-05Translation initiation factor activity
molecular_functionGO:000395011.71E-07NAD+ ADP-ribosyltransferase activity
molecular_functionGO:000396413.65E-04RNA-directed DNA polymerase activity
molecular_functionGO:000419714.77E-08Cysteine-type endopeptidase activity
molecular_functionGO:000425215.21E-21Serine-type endopeptidase activity
molecular_functionGO:000431613.66E-043-Oxoacyl-[acyl-carrier-protein] reductase (NADPH) activity
molecular_functionGO:000449718.16E-07Monooxygenase activity
molecular_functionGO:000455312.92E-06Hydrolase activity, hydrolyzing O-glycosyl compounds
molecular_functionGO:000455513.35E-05Alpha, alpha-trehalase activity
molecular_functionGO:000456512.51E-05Beta-galactosidase activity
molecular_functionGO:000474812.47E-04Ribonucleoside-diphosphate reductase activity, thioredoxin disulfide as acceptor
molecular_functionGO:000498411.05E-24Olfactory receptor activity
molecular_functionGO:000521513.82E-04Transporter activity
molecular_functionGO:000535113.61E-04Sugar:proton symporter activity
molecular_functionGO:000535519.37E-05Glucose transmembrane transporter activity
molecular_functionGO:000550614.78E-10Iron ion binding
molecular_functionGO:000554911.29E-27Odorant binding
molecular_functionGO:000823416.87E-11Cysteine-type peptidase activity
molecular_functionGO:000840813.61E-043′–5′ exonuclease activity
molecular_functionGO:000841714.30E-04Fucosyltransferase activity
molecular_functionGO:000842416.80E-05Glycoprotein 6-alpha-l-fucosyltransferase activity
molecular_functionGO:000852114.46E-25Acetyl-CoA transporter activity
molecular_functionGO:001502011.56E-15Glucuronosyltransferase activity
molecular_functionGO:001517418.07E-06Basic amino acid transmembrane transporter activity
molecular_functionGO:001518115.99E-05Arginine transmembrane transporter activity
molecular_functionGO:001518915.99E-05 l-lysine transmembrane transporter activity
molecular_functionGO:001529518.16E-07Solute:proton symporter activity
molecular_functionGO:001532618.07E-06Basic amino acid transmembrane transporter activity
molecular_functionGO:001670512.52E-16Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:001675814.79E-27Transferase activity, transferring hexosyl groups
molecular_functionGO:001687215.92E-04Intramolecular lyase activity
molecular_functionGO:001978911.64E-06SUMO transferase activity
molecular_functionGO:001979912.47E-06Tubulin N-acetyltransferase activity
molecular_functionGO:002003711.33E-16Heme binding
molecular_functionGO:002285711.26E-16Transmembrane transporter activity
molecular_functionGO:002289115.49E-05Substrate-specific transmembrane transporter activity
molecular_functionGO:003117713.66E-04Phosphopantetheine binding
molecular_functionGO:003377716.80E-05Lithocholate 6beta-hydroxylase activity
molecular_functionGO:004230214.06E-17Structural constituent of cuticle
molecular_functionGO:004316915.68E-05Cation binding
molecular_functionGO:004692112.03E-04Alpha-(1->6)-fucosyltransferase activity
molecular_functionGO:008001919.50E-07Fatty-acyl-CoA reductase (alcohol-forming) activity
cellular_componentGO:000586813.65E-04Cytoplasmic dynein complex
cellular_componentGO:001628111.57E-04Eukaryotic translation initiation factor 4F complex
cellular_componentGO:003258019.48E-04Golgi cisterna membrane
cellular_componentGO:003438813.65E-04Pwp2p-containing subcomplex of 90S preribosome
cellular_componentGO:004318613.58E-04P granule
biological_processGO:000008219.14E-04G1/S transition of mitotic cell cycle
biological_processGO:000072312.00E-11Telomere maintenance
biological_processGO:000597513.92E-04Carbohydrate metabolic process
biological_processGO:000599113.35E-05Trehalose metabolic process
biological_processGO:000635212.41E-04DNA-templated transcription, initiation
biological_processGO:000635318.24E-09DNA-templated transcription, termination
biological_processGO:000650812.72E-21Proteolysis
biological_processGO:000662919.50E-07Lipid metabolic process
biological_processGO:000685712.68E-11Oligopeptide transport
biological_processGO:000686512.47E-04Amino acid transport
biological_processGO:000709513.77E-05Mitotic G2 DNA damage checkpoint
biological_processGO:000717916.47E-05Transforming growth factor beta-receptor signaling pathway
biological_processGO:000760812.29E-11Sensory perception of smell
biological_processGO:000916613.82E-04Nucleotide catabolic process
biological_processGO:000945212.42E-077-Methylguanosine RNA capping
biological_processGO:000995313.48E-05Dorsal/ventral pattern formation
biological_processGO:001002519.50E-07Wax biosynthetic process
biological_processGO:001062912.41E-04Negative regulation of gene expression
biological_processGO:003357816.80E-05Protein glycosylation in Golgi
biological_processGO:003533619.50E-07Long-chain fatty-acyl-CoA metabolic process
biological_processGO:003542819.37E-05Hexose transmembrane transport
biological_processGO:004570519.37E-05Negative regulation of salivary gland boundary specification
biological_processGO:004632319.37E-05Glucose import
biological_processGO:004636812.03E-04GDP-l-fucose metabolic process
biological_processGO:005079011.45E-08Regulation of catalytic activity
biological_processGO:005090911.01E-05Sensory perception of taste
biological_processGO:005160311.36E-05Proteolysis involved in cellular protein catabolic process
biological_processGO:005269612.47E-04Flavonoid glucuronidation
biological_processGO:005508511.29E-14Transmembrane transport
biological_processGO:007192912.47E-06Alpha-tubulin acetylation
biological_processGO:190335215.99E-05 l-ornithine transmembrane transport
List of the GO Terms Enriched in the Duplicated Protein Families at the Base of All Six Aphids per Each Phylome.

The Ancestral Wave of Duplications Predates the Divergence of Aphids and Adelgids

Given the long branch subtending Aphididae, and to provide a narrower placement of the ancestral duplication wave, we expanded our taxonomic sampling by including the transcriptomes of two additional hemipteran insects from the suborder Sternorrhincha of taxonomic importance for our group of study: the adelgid (Adelgidae) Adelges tsugae (accession number: PRJNA242203) and the scale insect (Coccoidea) Paratachardina pseudolobata (Christodoulides et al. 2017). Most phylogenies show that the Adelgidae family is a sister group of the Phylloxeridae family (Heie and Wegierek 2009; Vilcinskas 2016). Thus, with the inclusion of Ad. tsugae we can obtain a general image of the Aphidomorpha lineage. With this increased species set, we reconstructed an expanded C. cedri phylome and species tree (fig. 3). The duplication analysis on the expanded data set initially resulted in two waves of duplications (fig. 3 one (0.18 duplications/gene) still specific to the Aphididae lineage, and another (0.36 duplications/gene) at the base of the Aphidomorpha lineage (Aphididae and Adelgidae families). However, it has been previously observed that, due to their incompleteness, transcriptomic data sets make difficult the correct placement of duplications (Jiménez-Guri et al. 2013). To account for this, we repeated the analysis considering only gene trees that included the two species with transcriptomic data sets. In this stringent set, the duplication ratio in the ancestral branch of the Aphidomorpha was still high (0.29), whereas the one in the branch subtending Aphididae disappeared (0.07). To confirm the presence of a single ancestral peak we again analyzed ratios of 4DTv for the pairs of paralogs mapped at the ancestral branch of the Aphidomorpha lineage, and for the orthologous pairs found between C. cedri and Ad. tsugae, and P. pseudolobata (see Materials and Methods). The distribution of the 4DTv values of paralogs and orthologs is used to estimate the relative order of duplication and speciation events, respectively. Consistent with our phylogenomic analyses described above, the ancestral peak of duplications is placed before the divergence of C. cedri and Ad. tsugae, and after the divergence of C. cedri and P. pseudolobata. From these results, we conclude that the large-scale gene duplication observed at the long branch subtending Aphididae in the full-genome data set occurred before the divergence of the Aphidomorpha group, and after the separation of this lineage from P. pseudolobata (Coccoidea). A dating analysis (see Materials and Methods) situates this duplication wave over a putative long temporal period, 106–227 Ma. As these times are molecular estimates, additional analysis should be necessary to place a more accurate time scale for this duplication event.
. 3.

Species tree, duplication ratio, and 4DTv of Cinara cedri. (a) Phylogenetic tree which included the two species with transcriptomes (Adelges tsugae and Paratachardina pseudolobata). Taxonomic groups are indicated by different colors. The bottom line represents the divergence time in Ma. Bars in the nodes indicate the uncertainty around mean age estimates based on 95% credibility intervals. All omitted bootstrap values are maximal (bootstrap 100%). Red dots mark the calibration points used to estimate the divergence times. (b) Zoom out of the Sternorrhyncha group. Duplication ratios are indicated in each branch for each set of gene trees: set 1 and set 2. The green star marks the position of the large-scale duplication event. (c) 4DTv of paralogous genes of C. cedri in the branch where the large-scale duplication event is marked in figure 3. 4DTv of orthologous pairs between C. cedri and Adelges tsugae, and P. pseudolobata are shown with different colors.

Species tree, duplication ratio, and 4DTv of Cinara cedri. (a) Phylogenetic tree which included the two species with transcriptomes (Adelges tsugae and Paratachardina pseudolobata). Taxonomic groups are indicated by different colors. The bottom line represents the divergence time in Ma. Bars in the nodes indicate the uncertainty around mean age estimates based on 95% credibility intervals. All omitted bootstrap values are maximal (bootstrap 100%). Red dots mark the calibration points used to estimate the divergence times. (b) Zoom out of the Sternorrhyncha group. Duplication ratios are indicated in each branch for each set of gene trees: set 1 and set 2. The green star marks the position of the large-scale duplication event. (c) 4DTv of paralogous genes of C. cedri in the branch where the large-scale duplication event is marked in figure 3. 4DTv of orthologous pairs between C. cedri and Adelges tsugae, and P. pseudolobata are shown with different colors. A functional analysis of the C. cedri proteins duplicated at the base of the Aphidomorpha lineage was largely consistent with the analysis of the data set that included only complete genomes (see above, table 3). A general analysis of both results (tables 3 and 4) shows many GO terms enriched with key functions for aphid, phylloxerid, and adelgid biology. These insects base their diets strictly on phloem sap, which requires very specific adaptations (Douglas 2006). In this regard, genes duplicated ancestrally are enriched in carbohydrate metabolism and metabolite transporters, which may be related to the need for efficient exploitation of phloem sap, which is rich in sugar but poor in other essential nutrients (Douglas 2006). Essential amino acids in aphids are provided by microbial symbionts (Baumann 2005). In this context, ancestral duplications are also enriched in amino acid transporters, which may allow reallocation of these essential nutrients to enhance the amino acid supply (Hansen and Moran 2011). Other important adaptations of phloem-sap feeding are the adaptation to plant secondary metabolites. Glutathione S-transferases play an important role in the detoxification of many substances including allelochemicals from plants (Francis et al. 2005). Genes associated with glutathione S-transferases are also duplicated at the ancestral branch of Aphidomorpha (e.g., see supplementary fig. S6a, Supplementary Material online). Similarly, genes duplicated ancestrally are enriched in functions related to UDP-glycosyltransferases, which are a major class of drug-metabolizing enzymes and play an important role in the detoxification of a large number of xenobiotics (Bock 2016). In aphids, UDP-glycosyltransferase may confer tolerance to thiamethoxam (Pan et al. 2015). Other ancestral duplications of genes are involved in wax biosynthesis, which may be related to maintaining water balance, and preventing desiccation (Chung and Carroll 2015). Enrichment in fatty acyl reductases (e.g., see supplementary fig. S6b, Supplementary Material online) may be related to not only wax biosynthesis but also components of insect cuticular hydrocarbons and pheromones (Tupec et al. 2019). Moreover, ancestral duplications of genes involved in smell and sugar taste perception may have facilitated detection of suitable plant host or development of alarm pheromones (Zhang et al. 2017). This ancestral wave of duplication is also enriched in functions related to growth and molting. One such example is the genes associated with ecdysis triggering hormone receptor (supplementary fig. S6c, Supplementary Material online), which are crucial for the activation of the ecdysis sequence (Roller et al. 2010). Other functional classes enriched among ancestral paralogs are associated with DNA and RNA processing, and may be fundamental to the maintenance of the genomic and phenotypic plasticity observed in aphids (Mathers et al. 2017).
Table 4.

List of GO Terms Enriched in the Duplicated Protein Families at the Base of the Aphidomorpha Group.

Term CategoryTerm Term Level Adj. P-ValueTerm Name
molecular_functionGO:000006412.38E-04 l-ornithine transmembrane transporter activity
molecular_functionGO:000373012.38E-04mRNA 3′-UTR binding
molecular_functionGO:000403219.45E-06Alditol:NADP+ 1-oxidoreductase activity
molecular_functionGO:000403519.08E-04Alkaline phosphatase activity
molecular_functionGO:000418513.29E-04Serine-type carboxypeptidase activity
molecular_functionGO:000419715.74E-16Cysteine-type endopeptidase activity
molecular_functionGO:000449718.22E-08Monooxygenase activity
molecular_functionGO:000455311.02E-15Hydrolase activity, hydrolyzing O-glycosyl compounds
molecular_functionGO:000455515.24E-08Alpha, alpha-trehalase activity
molecular_functionGO:000521512.31E-16Transporter activity
molecular_functionGO:000525413.45E-06Chloride channel activity
molecular_functionGO:000535112.46E-06Sugar:proton symporter activity
molecular_functionGO:000535515.52E-07Glucose transmembrane transporter activity
molecular_functionGO:000550618.41E-12Iron ion binding
molecular_functionGO:000554214.33E-07Folic acid binding
molecular_functionGO:000819412.43E-05UDP-glycosyltransferase activity
molecular_functionGO:000823419.57E-12Cysteine-type peptidase activity
molecular_functionGO:000823711.46E-06Metallopeptidase activity
molecular_functionGO:000826016.12E-043-Oxoacid CoA-transferase activity
molecular_functionGO:000842211.50E-04Beta-glucosidase activity
molecular_functionGO:000851814.33E-07Reduced folate carrier activity
molecular_functionGO:000852119.33E-25Acetyl-CoA transporter activity
molecular_functionGO:001502012.42E-15Glucuronosyltransferase activity
molecular_functionGO:001511611.50E-04Sulfate transmembrane transporter activity
molecular_functionGO:001517113.58E-05Amino acid transmembrane transporter activity
molecular_functionGO:001517413.81E-05Basic amino acid transmembrane transporter activity
molecular_functionGO:001518112.38E-04Arginine transmembrane transporter activity
molecular_functionGO:001518912.38E-04 l-lysine transmembrane transporter activity
molecular_functionGO:001529511.61E-09Solute:proton symporter activity
molecular_functionGO:001529712.33E-06Antiporter activity
molecular_functionGO:001532613.81E-05Basic amino acid transmembrane transporter activity
molecular_functionGO:001552819.87E-09Lactose:proton symporter activity
molecular_functionGO:001629816.12E-04Lipase activity
molecular_functionGO:001649116.42E-04Oxidoreductase activity
molecular_functionGO:001662012.71E-05Oxidoreductase activity, acting on the aldehyde or oxo group of donors, NAD or NADP as acceptor
molecular_functionGO:001670511.95E-14Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
molecular_functionGO:001671711.21E-06Oxidoreductase activity, acting on paired donors, with oxidation of a pair of donors resulting in the reduction of molecular oxygen to two molecules of water
molecular_functionGO:001674716.42E-04Transferase activity, transferring acyl groups other than amino-acyl groups
molecular_functionGO:001675815.44E-27Transferase activity, transferring hexosyl groups
molecular_functionGO:001678715.09E-05Hydrolase activity
molecular_functionGO:001711017.12E-06Nucleoside-diphosphatase activity
molecular_functionGO:001979912.63E-07Tubulin N-acetyltransferase activity
molecular_functionGO:002003718.50E-19Heme binding
molecular_functionGO:002285711.30E-38Transmembrane transporter activity
molecular_functionGO:002289113.70E-10Substrate-specific transmembrane transporter activity
molecular_functionGO:004262615.04E-07ATPase activity, coupled to transmembrane movement of substances
molecular_functionGO:004316914.04E-10Cation binding
molecular_functionGO:005066015.48E-04Flavin adenine dinucleotide binding
molecular_functionGO:005268918.97E-05Carboxylic ester hydrolase activity
molecular_functionGO:008001918.02E-07Fatty-acyl-CoA reductase (alcohol-forming) activity
molecular_functionGO:009048214.33E-07Vitamin transmembrane transporter activity
molecular_functionGO:010233612.97E-043-Oxo-arachidoyl-CoA synthase activity
molecular_functionGO:010233712.97E-043-Oxo-cerotoyl-CoA synthase activity
molecular_functionGO:010233812.97E-043-Oxo-lignoceronyl-CoA synthase activity
cellular_componentGO:000588718.49E-11Integral component of plasma membrane
cellular_componentGO:000853711.27E-06Proteasome activator complex
cellular_componentGO:001602112.85E-43Integral component of membrane
cellular_componentGO:004318611.97E-04P granule
cellular_componentGO:004323113.81E-05Intracellular membrane-bounded organelle
biological_processGO:000151018.02E-07RNA methylation
biological_processGO:000597514.78E-13Carbohydrate metabolic process
biological_processGO:000599115.24E-08Trehalose metabolic process
biological_processGO:000650811.84E-22Proteolysis
biological_processGO:000662916.45E-11Lipid metabolic process
biological_processGO:000663318.18E-04Fatty acid biosynthetic process
biological_processGO:000685719.96E-05Oligopeptide transport
biological_processGO:000686513.81E-05Amino acid transport
biological_processGO:000728314.73E-05Spermatogenesis
biological_processGO:000815212.80E-08Metabolic process
biological_processGO:000916613.12E-07Nucleotide catabolic process
biological_processGO:000945212.51E-077-Methylguanosine RNA capping
biological_processGO:001002511.73E-07Wax biosynthetic process
biological_processGO:001697319.45E-06Poly(A)+ mRNA export from nucleus
biological_processGO:003533611.73E-07Long-chain fatty-acyl-CoA metabolic process
biological_processGO:003542812.63E-07Hexose transmembrane transport
biological_processGO:004275916.18E-05Long-chain fatty acid biosynthetic process
biological_processGO:004632312.63E-07Glucose import
biological_processGO:005079013.19E-11Regulation of catalytic activity
biological_processGO:005118014.33E-07Vitamin transport
biological_processGO:005160311.40E-08Proteolysis involved in cellular protein catabolic process
biological_processGO:005269617.61E-06Flavonoid glucuronidation
biological_processGO:005508514.51E-56Transmembrane transport
biological_processGO:005511411.35E-15Oxidation-reduction process
biological_processGO:007050712.09E-04Regulation of microtubule cytoskeleton organization
biological_processGO:007192912.63E-07Alpha-tubulin acetylation
biological_processGO:190165719.80E-06Glycosyl compound metabolic process
biological_processGO:190335212.38E-04 l-ornithine transmembrane transport
List of GO Terms Enriched in the Duplicated Protein Families at the Base of the Aphidomorpha Group.

Use of a Chromosome-Level Assembly of A. pisum

Recently, a chromosome-level assembly for the A. pisum genome became available (Li et al. 2019). The presence of high numbers of paralogs complicates the process of genome assembly. Depending on the nature of the data and the assembly algorithms and parameters used, recently duplicated paralogs can be (partially) collapsed into a single sequence or, conversely, divergent alleles of the same loci can be separated into distinct sequences (Gabaldón and Alioto 2016). These issues had always been a concern when assessing the high levels of duplications in A. pisum and other aphid species (International Aphid Genomics Consortium 2010). To confirm our previous results with a more contiguous version of the assembly, we repeated the A. pisum phylome (A. pisum2, PhylomeID 707), this time using the annotation of the newly released, chromosome-level assembly (Li et al. 2019). The results of the analysis of this A. pisum2 phylome are in agreement with the results of the A. pisum1 phylome (PhylomeID 706). The duplication ratio at the ancestral branch of all aphids is still high (0.69). Also, the number of specific gene duplication events (1,825) is similar to that in the phylome 706 (1,889), and the percentage of proteins that have an in-paralog is the same (29%). However, the availability of the AL4 assembly allowed us to analyze the paralogs in the context of chromosomes. Interestingly, a high percentage of the duplicated proteins are present in the assembled chromosomes: For the aphid-ancestral wave of duplication, 82% of duplicated proteins are present in the chromosomes, and for the species-specific A. pisum wave, 76%. If we analyze the number of pairs of paralogs (expansions will form more than one pair) present in the same chromosome, a higher percentage was observed in the species-specific duplications (chromosome X—47%, A1—33%, A2—42%, A3—34%) with respect to the aphid ancestral duplications (chromosome X—34%, A1—28%, A2—19%, A3—12%). In both cases chromosome X has the highest percentage, whereas A3 has the lowest. Moreover, from the total number of proteins duplicated in both waves, chromosome X has the highest percentage, followed by chromosome A1 (fig. 4). These results are in agreement with recent studies (Li et al. 2019).
. 4.

Acyrthosiphon pisum paralogs mapped onto the chromosomes. (a) Percentage of proteins with paralogs at the two waves of duplications (age1: A. pisum species-specific, age6: ancestral to all aphids) per each chromosome. (b) Links between pairs of paralogs at age1. (c) Links between pairs of paralogs at age6. The chromosomes are shown in different colors: X—red, A1—blue, A2—gray, A3—green. Each chromosome was plotted independently for each age.

Acyrthosiphon pisum paralogs mapped onto the chromosomes. (a) Percentage of proteins with paralogs at the two waves of duplications (age1: A. pisum species-specific, age6: ancestral to all aphids) per each chromosome. (b) Links between pairs of paralogs at age1. (c) Links between pairs of paralogs at age6. The chromosomes are shown in different colors: X—red, A1—blue, A2—gray, A3—green. Each chromosome was plotted independently for each age. In order to analyze the distribution of paralogs along the chromosomes, we plotted all the pairs per chromosome and per wave of duplication (fig. 4). As previously noticed (Li et al. 2019), the paralogs are distributed throughout all the chromosomes. However, some blocks of paralogous pairs can also be observed. An evident block is shared between chromosomes X and A2 in both waves of duplications (fig. 4). Finally, we searched for footprints of these large-scale duplication events in the relative gene order of paralogs (i.e., synteny conservation), which rendered no significant result. A comparison of C. cedri against itself using Symap did not reveal any conserved region (see supplementary fig. S7, Supplementary Material online). The analysis was repeated by comparing C. cedri with A. pisum and Bemisia tabaci. We observed that some conserved blocks between C. cedri and A. pisum. However, when we compare C. cedri with Be. tabaci, a more distant relative, the gene order conservation disappears (see supplementary fig. S8, Supplementary Material online). Similarly a comparison of the chromosome-level assembly of A. pisum to Be. tabaci revealed no apparent conserved synteny block. We repeated the analysis using i-ADHore (Proost et al. 2012) and found few collinear segments (9) (see Materials and Methods). We repeated both analyses for the chromosome-level assembly of A. pisum with similar results. We also compared gene order between C. cedri and M. persicae and between C. cedri and the more distantly related Be. tabaci. The comparison between the first two showed a moderate gene order conservation, though the fragmentation of both genomes makes it difficult to assess whether the conservation is real or just an artifact. i-ADHore detected only 43 conserved segments between the two species, much less than the 509 segments found between more closely related A. pisum and M. persicae, indicating a quick degradation of synteny. Comparison between C. cedri and the less related Be. tabaci shows that all signs of gene order conservation have been lost between these two genomes, which is confirmed by i-ADHore. This result was also observed between the chromosome-level A. pisum and Be. tabaci, which indicates that the loss of gene order conservation is unlikely to be due to the fragmentation of the two genomes. The patterns observed in the 4DTv analysis (see above) indicate that the duplication event likely occurred soon after the divergence between aphids and P. pseudolobata. Therefore, if in fact there was a WGD, it is likely that further rearrangements and additional duplications have blurred the syntenic conservation between duplicated genes. These results may be influenced by the fragmentation of the genomes, as only the genome of A. pisum is at chromosome level. Yet, from our observations we do not see many differences when comparing C. cedri to A. pisum and M. persicae, leading us to believe that although we may be missing syntenic blocks and we could not provide an exact number of such blocks, the observed trend would remain similar if we had chromosome-level assemblies. Altogether, our results point to the presence of one major wave of ancestral duplications in the aphid lineage, predating the diversification of Aphidomorpha. This ancestral wave of duplications occurred in addition to other lineage-specific duplications, and to many recent species-specific duplications, highlighting a high genomic plasticity in aphids.

Conclusions

Recently, many WGDs have been described in insects (Li et al. 2018). Moreover, several independent studies have shown a burst of gene duplication in different species of aphids (Huerta-Cepas et al. 2010; Mathers et al. 2017; Li et al. 2018), but the origin of these duplications has been thus far unclear. Here we present the genome sequence of an early-diverging aphid (C. cedri) and its comparison with 20 complete animal genomes, including five sequenced aphid genomes, and the transcriptomes of two other phylogenetically important species. Taken together, our phylogenomic results provide compelling evidence for the existence of a large-scale gene duplication event predating the divergence of aphids, adelgids, and presumably phylloxerids (i.e., of Aphidomorpha). Genes duplicated at this large-scale event are enriched in functions that are relevant to aphids, adelgids, and phylloxerids, which share traits such as a phloem-sap-based diet or tight association with endosymbionts. The availability of the genome sequence of C. cedri, belonging to the subfamily Lachninae, an early-branching lineage within Aphididae, which feeds on gymnosperm, helps situate other major genomic events, such as the loss of selenoproteins and the streamlining of the immune repertoire closer to the base of the aphid lineage. Additional genome sequences, particularly of species of the Phylloxeridae family would help to confirm these findings and provide additional evolutionary insights. Our results underscore the use of phylogenomic approaches to study ancient duplication events (Marcet-Houben and Gabaldón 2015; Julca et al. 2018), and we do not discard the possibility of at least one ancestral WGD in the evolutionary history of aphids, even in the absence of syntenic conservation. Although synteny has been the traditional approach to uncover ancient WGDs (Ohno 1970; Wolfe 2015), it relies on a signal that is blurred by subsequent genomic rearrangements and may not apply equally well to different lineages. It has been suggested that the holocentric nature of aphid chromosomes may promote rapid reshuffling of gene order (Blackman 1980; Mandrioli et al. 2016; Li et al. 2019), and this is consistent with our findings. In our case, the lack of syntenic conservation between Be. tabaci and the aphids indicates that—at the level of resolution of our methods and the current level of completion of the compared genomes—a large number of rearrangements have occurred. This, coupled with large amounts of lineage-specific gene expansions, has likely degraded the initial syntenic relationships originated at the putative ancestral wave of gene duplications.

Materials and Methods

Heterozygosity Analysis and DNA Extraction

A C. cedri population collected from a single cedar tree (Cedrus libani) in the spring of 2011 in Lliria (Valencia, Spain, 39°38′47,0826″N, 0°37′ 51,3006″O) was introduced and maintained on a cedar tree at the facilities of the Institut Cavanilles de Biodiversitat I Biologia Evolutiva at the University of Valencia (Spain). For the DNA extraction, all the individuals were collected during March 2012 from a colony from the tree, and because they are parthenogenetic, it is likely that all come from the same female. To analyze heterozygosity in the population, the cytochrome oxidase I (COI) mitochondrial gene of 45 individual insects, collected between July 2011 and March 2012, was amplified by polymerase chain reaction (PCR) using the primers LCO1490 and HCO21980 (Folmer et al. 1994) and sequenced. All analyzed sequences were identical. Approximately 50–60 apteral adult females were gently homogenized and used for DNA extraction, using the “Purification of DNA from insects using the DNeasy Blood & Tissue Kit” (Quiagen) according to the manufacturer’s instructions.

Flow Cytometric Genome Size Estimates

Flow cytometric genome size was estimated for nuclei isolated from whole insect tissue stained with propidium iodide and scored by flow cytometry for relative red fluorescence of the 2C peaks of the sample and a co-prepared standard (see Johnston et al. 2019). In brief, the anterior one-third of each sample was placed into 1 ml of Galbraith buffer in a 2-ml Kontes Dounce tissue homogenizer along with the head of a lab-reared Drosophila virlis standard (1C = 328 Mb). Nuclei were released from the sample and standard by 15 strokes of the “A” Dounce pestle. The resultant solution was filtered through a 40-µm nylon mesh, stained with 25 mg/ml propidium iodide and held for 2 h in the dark and cold to allow comparable levels of dye saturation in the sample and standard. The average channel number of the 2C peaks of the sample and standard were scored and the 1C (gametic) DNA calculated as the ratio of the 2C peaks of the sample and standard times 328 Mb. At least 500 nuclei were scored under each peak with a maximum coefficient of variation (CV) of 3.0 for each. The mean and standard error were based on estimates for individuals and standards in separate co-preparations.

Genome Sequencing and Assembly

The short-insert paired-end libraries were prepared with the NO-PCR protocol. TruSeqDNA Sample Preparation Kit v2 (Illumina Inc.) and the KAPA Library Preparation Kit (Kapa Biosystems) were used. In short, 2.0 µm of sheared genomic DNA was end-repaired, adenylated, and ligated to Illumina-specific indexed paired-end adaptors. The DNA was size selected with AMPure XP beads (Agencourt, Beckman Coulter) in order to reach the fragment size of 220–550 bp. The final libraries were quantified by Library Quantification Kit (Kapa Biosystems). The library was sequenced using the TruSeq SBS Kit v3-HS (Illumina Inc.) in paired-end mode, 2 × 101 bp, in one sequencing lane of a HiSeq2000 (Illumina Inc.) according to standard Illumina operation procedures with a yield of >30 Gb of raw data. Primary data analysis, image analysis, base calling, and quality scoring of the run were processed using the manufacturer’s software Real Time Analysis (1.13.48) and followed by the generation of FASTQ sequence files by CASAVA. In addition, two Mate Pair (MP) libraries with insert sizes of 3,000 (MP3000) and 5,000 bp (MP5000) were constructed according to a modified Illumina protocol. In brief, after genomic DNA fragmentation, circularization of the DNA was performed in the presence of a biotinylated 454 double-stranded linker. Thereafter, the standard Illumina matepair preparation method was followed. The libraries were sequenced on the Illumina HiSeq2000 platform in paired-end mode, which outputs 101-bp reads (2 × 101 bp). The yield was at least 11 Gb for both MP libraries. Postprocessing of sequence reads involved trimming of the linker sequence (TCGTATAACTTCGTATAATGTATGCTATACGAAGTTATTACG and reverse complement) using cutadapt v2.5 (Martin 2011) with -e 0.05 and -O 10 options. Only pairs for which at least one mate was trimmed (i.e., contained the linker and was thus a true mate pair (MP) and not paired-end (PE) contamination) were kept for scaffolding. Then, we used gem-mapper (Marco-Sola et al. 2012) to detect reads matching contaminants. Contaminated reads were filtered before the assembly. The genome size and complexity were estimated using Jellyfish v1.1 (Marçais and Kingsford 2011) and GenomeScope v1.0 (Vurture et al. 2017), resulting in estimates of 508.6 and 399.76 Mb, respectively. The latter estimate was used to guide our assembly strategy. The genome was assembled as follows. SGA preqc analysis (Simpson 2014) was used to determine optimal K-mer length (k = 65) for de Bruijn graph construction. Bracketing around this optimum, multiple assemblies with shorter and longer K-mers were made using ABySS v1.5.2 (Simpson et al. 2009) merged with ASM, an OLC-like assembly-merging software to obtain contigs (Frias and Ribeca, 2016; Cruz et al. 2016). The merging parameters were: –anchor 125 –anchor-spacing 10 –min-anchor 25 –coverage 2 –divergence 0.03 –anchorsxchunk 50000000 –repeat-resolution-depth 0 –path-expansion-depth 0 –consensus-type majority. The sequencing libraries were mapped with gem-mapper (Marco-Sola et al. 2012), and these mappings were used for scaffolding the merged contigs with ABySS using parameters -n 4 -s 200 -N 8 -S 130-2000 -k 67 -l 36 -q 10. This intermediate assembly was refined by decontamination, consistency check (Cruz et al. 2016), and then discarding scaffolds shorter than 4 kb already contained in longer scaffolds (i.e., unique mappings with 0% mismatches detected with gem-mapper). This refined assembly was re-scaffolded with SSPACEv3.0 (Boetzer et al. 2011; Cruz et al. 2016), and gaps closed with GapFiller (Boetzer and Pirovano 2012). Afterwards, the PE library (PE400) was mapped against the target assembly (gap-filled and decontaminated) with BWA mem v0.7.7 (Li and Durbin 2009) and then performed variant calling with samtools v0.1.19 mpileup (Li et al. 2009). Only single-nucleotide substitutions or indels having at least ten reads supporting the alternative allele were used to produce an alternative reference. The alternative FASTA sequence was obtained using GATK v3.5 FastaAlternateReferenceMaker (McKenna et al. 2010). Finally, the assembly was named internally cinced3, and the gene completeness was evaluated with CEGMA v2.4 (Parra et al. 2007), which searches for 248 core-eukaryotic genes, and BUSCO v3.0.2 (Simão et al. 2015) using the insecta_odb9 that includes 42 species and 1,658 genes.

RNA Extraction and Transcriptome Sequencing

RNA samples were prepared from 400 adults of C. cedri females from the clonal population mentioned above. Samples were obtained from aphid heads (absence of endosymbionts) and dissected bacteriocytes (presence of the two endosymbionts, B. aphidicola and Se. symbiotica). Total RNA extraction was performed using the “TRI Reagent Solution” Kit (Ambion), following the manufacturer’s instructions. In addition, two SOLID libraries were prepared from the same tissues.

Genome Annotation

A combination of the Program to Assemble Spliced Alignments (PASA v2.0.2) and Evidence Modeler (EVM v1.1.1) (Haas et al. 2008) was used to obtain consensus coding sequence models using three main sources of evidence: aligned transcripts, aligned proteins, and gene predictions. Finally, noncoding RNAs were annotated employing CMsearch (Cui et al. 2016), tRNAscan-SE (Lowe and Eddy 1997), and lncRNAs were obtained from the PASA-assemblies without protein-coding gene annotations that were longer than 200 bp. The Piper-R NF pipeline (https://github.com/cbcrg/piper-nf) was used to detect lncRNA conservation between C. cedri and 20 other metazoans (supplementary table S3, Supplementary Material online). Then, a heatmap was plotted including all the lncRNAs that are present in at least one of the other species using the R package gplots (Warnes et al. 2016).

Identification of Genes of the Immune System of C. cedri

A database of the genes involved in the immune system of well-studied insects was generated by updating and expanding the database from Insect Innate Immunity Database (Brucker et al. 2012). The source references for the jewel wasp (Nasonia vitripennis [Werren et al. 2010]), the honey-bee (Apis mellifera [Evans et al. 2006]), the fruit fly (Drosophila melanogaster [De Gregorio et al. 2001]), the African malaria mosquito (Anopheles gambiae [Christophides et al. 2002]), and the pea aphid (A. pisum [Gerardo et al. 2010]) were revised, and the genes they described were retrieved. To increase completeness, the immune system genes from the red flour beetle (Tribolium castaneum [Zou et al. 2007]), the diamondback moth (Plutella xylostella [Xia et al. 2015]), the tobacco hornworm (Manduca sexta [Cao et al. 2015]), the head louse (Pediculus humanus [Kang et al. 2015]), the Florida carpenter ant (Camponotus floridanus [Gupta et al. 2015]), the Asian citrus psyllid (Diaphorina citri [Arp et al. 2016]), and the silkworm (Bombyx mori [Tanaka et al. 2008]) were also added. EggNOG-mapper (Huerta-Cepas et al. 2017) was used to identify orthologs among the selected species using the Arthropoda (artNOG) data set followed by manual curation where discrepancies were observed.

Aphid Phylomes Reconstruction

The phylomes of C. cedri, A. pisum, M. persicae, D. noxia, Ap. glycines, and S. flava were reconstructed using the PhylomeDB pipeline (Huerta-Cepas et al. 2011). For A. pisum, two phylomes were reconstructed using the two proteomes available (Acyr 2.0 [International Aphid Genomics Consortium 2010] and AL4 [Li et al. 2019]). In brief, for each protein-coding gene in each aphid genome we searched for homologs (Smith–Waterman BLAST search, e-value cutoff <1e-05, minimum contiguous overlap over the query sequence cutoff 50%) in a protein database containing the proteomes of 21 species for C. cedri and a subset of 14 species for the other aphids (supplementary table S2, Supplementary Material online). The most similar 150 homologs were aligned using three different programs (MUSCLE [Edgar 2004], MAFFT [Katoh et al. 2005], and KALIGN [Lassmann and Sonnhammer 2005]) in both forward and reverse orientations. These six alignments were combined using M-COFFEE (Wallace et al. 2006) and trimmed with trimAl v.1.3 (Capella-Gutiérrez et al. 2009) using a consistency cutoff of 0.16667 and a gap threshold of 0.1. Phylogenetic trees were built using maximum likelihood approach as implemented in PhyML v3.0 (Guindon et al. 2010) using the best fitting model among seven different ones (JTT, LG, WAG, Blosum62, MtREV, VT, and Dayhoff). The two models best fitting the data were determined based on likelihoods of an initial neighbor joining tree topology and using the AIC criterion. We used four rate categories and inferred fraction of invariant positions and rate parameters from the data. All alignments and trees are available for browsing or download at PhylomeDB with the PhylomeID: C. cedri—701, S. flava—702, Ap. glycines—703, D. noxia—704, M. persicae—705, A. pisum 1—706, and A. pisum2—707 (Huerta-Cepas et al. 2014).

Prediction of Gene Duplications, and Orthology and Paralogy Relationships

Orthology and paralogy relationships were predicted based on phylogenetic evidence from each aphid phylome. We used ETE v3.0 (Huerta-Cepas et al. 2016) to infer duplication and speciation relationships using the species overlap method (Gabaldón 2008) and a topology-based phylostratigraphic method to date duplication events (Huerta-Cepas and Gabaldón 2011). In brief the species-overlap algorithm identifies internal nodes as duplications if the two daughter clades show any overlap in species, and the phylostratigraphic method assigns a relative age to that duplication as the last common ancestor of all the taxa contained in the two daughter clades. Species-specific duplications (expansions) were computed as duplication that map only to the seed species of each phylome (C. cedri, A. pisum, M. persicae, D. noxia, and Ap. glycines). In order to reduce redundancy, expansions that overlap in more than 50% of their sequences were fused together using a UPGMA clustering. Duplication ratios were calculated by dividing the number of duplications mapped to a given node in the species tree by all the gene trees that contain that node. In all the cases, duplication frequencies that include expansions larger than five in each phylome were excluded. Due to ancestral expansions in other species that can affect the duplication ratio, an additional filter was applied by removing all the gene trees that contained aphid duplications with more than five sequences in any of the aphid species included (C. cedri, A. pisum, M. persicae, D. noxia, Ap. glycines, and S. flava). Duplication ratios were calculated again using this new subset of gene trees. All orthology and paralogy relationships are available through PhylomeDB (Huerta-Cepas et al. 2014).

Incorporation of Transcriptomic Data in the Phylome

In order to increase the taxonomic sampling in the Sternorrhyncha group we decided to include two species, which have their transcriptome available: Ad. tsugae (accession number: PRJNA242203) and P. pseudolobata (Christodoulides et al. 2017). Adelges tsugae belongs to the family Adelgidae, a clade inside the Aphidomorpha. Paratachardina pseudolobata is a scale insect that belongs to the superfamily Coccoidea. The transcriptome of P. pseudolobata was obtained from the whole body of female samples, and the assembly is mostly complete according to BUSCO (89% of highly conserved arthropod sequences were present as single-copy or duplicated transcripts in the assembly) (Christodoulides et al. 2017). The transcriptome of Ad. tsugae was also obtained from the whole body of female samples. Although these data are not published, the gene completeness was evaluated with BUSCO 3.1.0 (Simão et al. 2015) and we found that 79% of the conserved genes of the insecta_odb9 data set were present as single-copy or duplicated transcript in this assembly. The Transcriptome Shotgun Assembly (TSA) file of both species was downloaded, and the prediction of proteins was obtained by selecting the longest open reading frame for each transcript (>100 amino acids). The incorporation of the transcriptomic data was done using the following pipeline. First, a BLASTP was performed from the seed protein against a database that contained the two transcriptomes. Then, the results were filtered based on three thresholds: e-value < 1e-05. Overlap between query and hit had to be at least 0.3, and sequence identity > 40.0%. Proteins that passed these filters were incorporated into the raw alignment of the phylome using MAFFT v7.222 (–add and –reorder options) (Katoh et al. 2005). Then, gene trees were reconstructed using this new alignment following the same procedure as described above. Finally, these gene trees were filtered in order to remove unreliably placed transcriptome sequences (set1) and filtered again to keep only trees that contained both species, Ad. tsugae and P. pseudolobata (set2).

GO Term Enrichment

GO terms were assigned to the five aphid proteomes using Interproscan v5.34-73.0 (Jones et al. 2014) and the annotation of orthologs from the PhylomeDB database (Huerta-Cepas et al. 2014). GO term enrichment analysis was performed in each phylome using an in-house adaptation of FatiGO (Al-Shahrour et al. 2007) by comparing annotations of the duplicated proteins specific to the aphid (and at the branches subtending the wave of duplications) against all the other proteins encoded in the aphid genome.

Species Tree Reconstruction

All the gene trees that contained one-to-one orthologs per each of the 21 species included in the C. cedri phylome were kept (57 gene trees). The alignments used for the reconstruction of these gene trees were concatenated, and a species tree was reconstructed using the amino acid substitution model LG implemented in RAxML v8.1.17 (Stamatakis 2014) and 1,000 bootstrap replicates. In addition, a species tree including the two transcriptome data (Ad. tsugae and P. pseudolobata) was also obtained by concatenating, from the previous set, the gene trees that included at least one of the new species (52 gene trees).

Transversion Rate at 4-Fold Degenerate Sites

The distribution of the 4DTv was used to estimate the relative age of speciation events and duplications. Per each phylome, the 4DTv of pairs of paralogs at the branches subtending the high duplication rates (branches marked as A and B, fig. 2) were calculated. Also, for the C. cedri phylome we included the 4DTv of orthologs of C. cedri with S. flava, and Be. tabaci. For the other five phylomes (S. flava, Ap. glycines, D. noxia, M. persicae, A. pisum1) the 4DTv of orthologs for each aphid with Be. tabaci, and C. cedri were also included. After the two transcriptome data sets were included, the 4DTv was calculated for pairs of paralogous genes of C. cedri at the base of the Aphidomorpha group, and pairs of orthologous genes between C. cedri and P. pseudolobata, and Ad. tsugae.

Divergence Times

To place a time scale on the maximum likelihood phylogeny, we used the Bayesian-relaxed molecular clock approach as implemented in PhyloBayes v4.1c (Lartillot et al. 2013). An uncorrelated relaxed clock model was applied, and six constraints specified to the most recent common ancestor, including fossil specimens and secondary calibrations: Aphidomorpha (135 Ma, fossil evidence [Heie 1987; Havill et al. 2007]), Aphididae (80–100 Ma, previous molecular date estimates [Von Dohlen 2000] and fossil remains [Heie 1987, 1999; Heie and Wegierek 2011]), Aphidinae (70 Ma, fossil record [Heie 1987; Hong 2002]), Hexapoda (425 Ma, fossil evidence [Grimaldi and Engel 2005]), Holometabola (300 Ma, fossil evidence [Labandeira and Phillips 1996]), and Acari (410 Ma, fossil record [Hirst 1923; Dubinin 1962; Dunlop and Selden 2009]). These calibration constraints were used with soft bounds (Yang and Rannala 2006) under a birth–death prior, and a prior on the root of the tree (a mean of 560 Ma and a standard deviation of 100 Ma) (Liu et al. 2014). Two independent MCMC chains were run for 100,000 cycles, sampling posterior rates and dates every ten cycles. The initial 25% were discarded as burn-in. Posterior estimates of divergence dates and associated 95% credibility intervals were then computed from the remaining samples of each chain.

Physical Mapping of A. pisum Paralogs into the Chromosomes

Paralogs obtained from the A. pisum2 phylome (PhylomeID 707) were mapped into the four chromosomes (X, A1, A2, and A3). Links between the pairs of paralogs were visualized using Circos v0.69-6 (Krzywinski et al. 2009).

Synteny

SyMap v4.2 (Soderlund et al. 2011) was used to search for gene order conservation within the C. cedri genome, between C. cedri and its ancestor Be. tabaci, between C. cedri and A. pisum and between A. pisum and Be. tabaci. Only scaffolds of 0.5 Mb or longer were used for the comparison. Default parameters were used. i-ADHore (Proost et al. 2012) was also used to search for collinear segments within genomes and between different genomes. It was run for the individual genomes of C. cedri and the chromosome-level assembly of A. pisum. Then it was run between C. cedri and M. persicae and C. cedri and Be. tabaci, and finally, repeated between A. pisum and the same three species. i-ADHore was run using gg2 as the alignment method, with a tandem_gap of 75, gap_size of 30, cluster_gap of 35, q_value = 0.75, prob_cutoff = 0.01, and anchor_points = 3.

Data Availability

The assembly and annotation are also hosted at http://denovo.cnag.cat/ccedri where a jbrowse genome browser with most of the data can also be accessed. Click here for additional data file.
  105 in total

1.  Scaffolding pre-assembled contigs using SSPACE.

Authors:  Marten Boetzer; Christiaan V Henkel; Hans J Jansen; Derek Butler; Walter Pirovano
Journal:  Bioinformatics       Date:  2010-12-12       Impact factor: 6.937

2.  Evolution of host specialization in the Adelgidae (Insecta: Hemiptera) inferred from molecular phylogenetics.

Authors:  Nathan P Havill; Robert G Foottit; Carol D von Dohlen
Journal:  Mol Phylogenet Evol       Date:  2006-11-17       Impact factor: 4.286

3.  Glutathione S-transferases in the adaptation to plant secondary metabolites in the Myzus persicae aphid.

Authors:  Frédéric Francis; Nicolas Vanhaelen; Eric Haubruge
Journal:  Arch Insect Biochem Physiol       Date:  2005-03       Impact factor: 1.698

4.  Expansion of the miRNA pathway in the hemipteran insect Acyrthosiphon pisum.

Authors:  Stéphanie Jaubert-Possamai; Claude Rispe; Sylvie Tanguy; Karl Gordon; Thomas Walsh; Owain Edwards; Denis Tagu
Journal:  Mol Biol Evol       Date:  2010-02-23       Impact factor: 16.240

Review 5.  The UDP-glycosyltransferase (UGT) superfamily expressed in humans, insects and plants: Animal-plant arms-race and co-evolution.

Authors:  Karl Walter Bock
Journal:  Biochem Pharmacol       Date:  2015-10-08       Impact factor: 5.858

6.  GenomeScope: fast reference-free genome profiling from short reads.

Authors:  Gregory W Vurture; Fritz J Sedlazeck; Maria Nattestad; Charles J Underwood; Han Fang; James Gurtowski; Michael C Schatz
Journal:  Bioinformatics       Date:  2017-07-15       Impact factor: 6.937

7.  Comparative transcriptomics of early dipteran development.

Authors:  Eva Jiménez-Guri; Jaime Huerta-Cepas; Luca Cozzuto; Karl R Wotton; Hui Kang; Heinz Himmelbauer; Guglielmo Roma; Toni Gabaldón; Johannes Jaeger
Journal:  BMC Genomics       Date:  2013-02-24       Impact factor: 3.969

8.  The draft genome of whitefly Bemisia tabaci MEAM1, a global crop pest, provides novel insights into virus transmission, host adaptation, and insecticide resistance.

Authors:  Wenbo Chen; Daniel K Hasegawa; Navneet Kaur; Adi Kliot; Patricia Valle Pinheiro; Junbo Luan; Marcus C Stensmyr; Yi Zheng; Wenli Liu; Honghe Sun; Yimin Xu; Yuan Luo; Angela Kruse; Xiaowei Yang; Svetlana Kontsedalov; Galina Lebedev; Tonja W Fisher; David R Nelson; Wayne B Hunter; Judith K Brown; Georg Jander; Michelle Cilia; Angela E Douglas; Murad Ghanim; Alvin M Simmons; William M Wintermantel; Kai-Shu Ling; Zhangjun Fei
Journal:  BMC Biol       Date:  2016-12-14       Impact factor: 7.431

9.  Fast and accurate short read alignment with Burrows-Wheeler transform.

Authors:  Heng Li; Richard Durbin
Journal:  Bioinformatics       Date:  2009-05-18       Impact factor: 6.937

10.  Insect Innate Immunity Database (IIID): an annotation tool for identifying immune genes in insect genomes.

Authors:  Robert M Brucker; Lisa J Funkhouser; Shefali Setia; Rini Pauly; Seth R Bordenstein
Journal:  PLoS One       Date:  2012-09-12       Impact factor: 3.240
View more
  6 in total

1.  PhylomeDB V5: an expanding repository for genome-wide catalogues of annotated gene phylogenies.

Authors:  Diego Fuentes; Manuel Molina; Uciel Chorostecki; Salvador Capella-Gutiérrez; Marina Marcet-Houben; Toni Gabaldón
Journal:  Nucleic Acids Res       Date:  2022-01-07       Impact factor: 16.971

2.  Chromosome-level genome assembly for the horned-gall aphid provides insights into interactions between gall-making insect and its host plant.

Authors:  Hong-Yuan Wei; Yu-Xuan Ye; Hai-Jian Huang; Ming-Shun Chen; Zi-Xiang Yang; Xiao-Ming Chen; Chuan-Xi Zhang
Journal:  Ecol Evol       Date:  2022-04-21       Impact factor: 3.167

3.  Detoxification Gene Families at the Genome-Wide Level of Rhus Gall Aphid Schlechtendalia chinensis.

Authors:  Hongli He; M James C Crabbe; Zhumei Ren
Journal:  Genes (Basel)       Date:  2022-09-10       Impact factor: 4.141

4.  Improved Genome Assembly and Annotation of the Soybean Aphid (Aphis glycines Matsumura).

Authors:  Thomas C Mathers
Journal:  G3 (Bethesda)       Date:  2020-03-05       Impact factor: 3.154

5.  Genome Sequence of the Banana Aphid, Pentalonia nigronervosa Coquerel (Hemiptera: Aphididae) and Its Symbionts.

Authors:  Thomas C Mathers; Sam T Mugford; Saskia A Hogenhout; Leena Tripathi
Journal:  G3 (Bethesda)       Date:  2020-12-03       Impact factor: 3.154

6.  Chromosome-Scale Genome Assemblies of Aphids Reveal Extensively Rearranged Autosomes and Long-Term Conservation of the X Chromosome.

Authors:  Thomas C Mathers; Roland H M Wouters; Sam T Mugford; David Swarbreck; Cock van Oosterhout; Saskia A Hogenhout
Journal:  Mol Biol Evol       Date:  2021-03-09       Impact factor: 16.240
  6 in total

北京卡尤迪生物科技股份有限公司 © 2022-2023.