Rhus gall aphids (Hemiptera: Aphididae: Eriosomatinae: Fordini) only use Rhus species (Anacardiaceae) as their primary host plants, and each aphid specifically lives on only one or two sister host plants with the same disjunct distribution pattern between Eastern Asia and Eastern North America. We assembled complete mitochondrial genomes and universal single-copy nuclear genes for Rhus gall aphids using a genome skimming method, and the phylogenetic relationships strongly supported the monophyly of the Rhus gall aphids and two genera, Floraphis and Melaphis. However, the relationships among genera are inconsistent across different datasets. The cpDNA analysis on the host plant Rhus species strongly supported the Rhus monophyly and relationship of the Rhus interspecies. Cophylogeny analysis indicated that the origin of Rhus gall aphids was earlier than that of their host plants. However, the divergence time and relationships among some Rhus gall aphid species, particularly those with later divergence times, were consistent with the origin of their corresponding primary host plants. This may suggest that the Rhus gall aphids established the initial association with the ancestors between Rhus and the related groups, or acquired current Rhus hosts through host switching. The divergence time estimations imply that the separation of North America and Eurasia in the Laurasia supercontinent and the disappearance of the Bering Land Bridge, respectively, have played an important role in the divergence of the eastern North American Melaphis and an East Asian lineage. Our current results provide new insights into the coevolution of insects and host plants.

We collected 45 samples representing 14 Rhus gall aphid species and three subspecies and some closely related aphids. We extracted genomic DNA using three to five individuals from each gall sample with the TIANamp Genomic DNA Kit (TIANGEN) following the manufacturer’s instructions. We sent the genomic DNAs of all the samples to the Genomic Sequencing and Analysis Co. (Shanghai, China) for library construction and sequencing using the shotgun genome skimming method on an Illumina HiSeq 4000 platform. Paired-end reads of 2×150 bp were generated with an insert size of 400 bp. After filtering out low-quality and adapter-contaminated reads with Trimmomatic v. 0.35 , high-quality clean reads were used to do de novo  assemblies using the program SPAdes v.3.10.1  with the kmers 21, 33, 55, 77, 99 and 127. We used the mapping method in the software Geneious 10.2.4 to obtain the complete mitochondrial genome from contigs (Kearse et al. 2012). We obtained the complete mitogenome sequence by aligning, comparing and manually checking to the available consensus sequences. We captured single-copy nuclear genes directly from contigs, which was assembled by SPAdes v.3.10.1, and using OrthoDB database of Hemiptera with default parameters. Then we used MAFFT for alignment and obtained the alignment results.