DNA extraction and Wolbachia infection rate
Tuberculatus aphids were collected from regions throughout Japan and around Mt. Kariwangsan of South Korea (Table S2 and Figures 1 and 2). A species was regarded as ant-attended if aphids offered honeydew directly from their anus to attending ants. Eleven ant-attended and 12 non-attended species were obtained (Table 1 and Table S1). Because it was difficult to identify the three ant-attended aphid species (T. fulviabdominalis, T. indicus, and T. pilosulus) and the seven non-attended aphid species (T. higuchii A- and B-types, T. kashiwae A- and B-types, T. yokoyamai, T. sp. D, and T. sp. F), those species were identified through the sequence (Table S1). Sampling was conducted on viviparous females, which appears from April to September. Since Tuberculatus aphid species parthenogenetically produce nymphs in summer, several nymph individuals on a leaf are a high possibility of a clone. Therefore, aphids were collected from several leaves in a tree, to avoid collecting clonal aphids. They were collected in 99.5% ethanol and stored at -20 ℃. Before DNA extraction, the collected aphids were dissected to check for the presence of parasitoid wasps. Total DNA was extracted from each dissected aphid (whole body) with the Wizard genomic DNA purification kit (Promega, Tokyo, Japan). Since the 16S rRNA gene is highly conserved in a wide variety of microorganisms, it was used for polymerase chain reaction (PCR) amplification to determine the presence or absence of Wolbachia. In the small-scale experiment, using a gene map of the 16S rRNA locus of Wolbachia (Simões et al 2011), seven pairs of primers were selected and tested for each of 23 Tuberculatus species, in which two to three individuals per species were tested (Table 3). One pair of primers, 16SWolbF (16S-3f) (Casiraghi et al 2001) and WspecR (16S-2r) (Werren and Windsor 2000), was identified as the most appropriate for assessing the 23 species (Table 3) because it was able to amplify Wolbachia at the maximum number of species (seven species) of the 23 species. After the small-scale experiment, a full-scale experiment using the pair of primers was conducted on all collected samples (Table 1). To check whether DNA extraction was successful, the barcoding region (in mitochondrion) of primer pairs, LCO1490 and HCO2198, was also used (Table 2). Because more than 90% of individuals of T. macrotuberculatus in the Ishikari site harboured Wolbachia (Yao 2019), one individual of the species from the site was used for a positive control sample for Wolbachia detection. PCR was performed in 10 µL volume which included 2 µL of 5×KAPATaq Extra buffer (Nippon Genetics, Tokyo, Japan), 1 µL 25 mM MgCl2, 0.3 µL dNTP mixture (10 mM of each), 0.5 µL of 10 µM of each primer, 1 µL template DNA, and 0.05 µL KAPATaq Extra DNA polymerase (5 units/ µL). Reaction cycle parameters were: 94°C for 1 min; 40 cycles of 94°C for 20 sec, 50°C for 20 sec, and 68°C for 1 min, followed by a final extension of 68°C for 1 min. When PCR products had faint bands, the samples were rechecked by PCR in 20µL reaction volume. If the bands were false, nothing was amplified in 20µL reaction volume. The PCR product was checked using 1.5% agarose gel electrophoresis with ethidium bromide stain illuminated by UV light. The Wolbachia infection rate of each species was defined as the percentage of individuals amplified with the Wolbachia-specific primer out of all individuals amplified with the barcoding region primer. The correlation between the Wolbachia infection ratio in each collection site and geographical distance was tested by Mantel test (Mantel 1967) using the package vegan in R (R Development Core Team 2008). Mantel test was applied to the species that was collected from more than a single site and showed mean infection rates of less than 100%.
Phylogenetic trees for Tuberculatus aphids
A phylogenetic tree of the 23 Tuberculatus aphid species was constructed from the nucleotide sequences of a mitochondrion gene of a partial of cytochrome oxidase subunit I (COI) (940bp) from DDBJ (DNA Data Bank of Japan) (Table 1). Besides the COI gene, a partial of the nuclear gene of 18S rRNA (approx. 670bp) was amplified and used to construct phylogenetic trees. For reading the sequences of 18S rRNA gene, PCR was performed in 20µL reaction volume with a pair of primers (Ns1 and Ns2a; Table 2), the same reagents, and reaction cycles, as mentioned in the previous section, but changed annealing temperature to 47°C. PCR products were purified and sent to sequencing service (using Sanger sequencing) (Eurofin, Japan). The sequence data of the 18S rRNA gene (515bp) were deposited in the DDBJ and accession numbers were listed in Table 1. In addition to sequences of COI and 18S rRNA genes, a combined sequence of COI and 18S rRNA genes (1,455bp) was used for the construction of phylogenetic trees. The appropriateness of the combined sequence was checked by a homogeneity test implemented in PAUP* (P > 0.05). Judging from the phylogenetic position of Calaphidini aphid species (Lee et al 2017, 2021), specimens of two aphid species, Tinocallis zelkowae and a species of Takecallis, were included as outgroups. Neighbor-joining (NJ), most parsimonious (MP), and maximum likelihood (ML) analyses were performed using PAUP* 4.0b10 PPC and 4.0a 169 (Swofford 2002). For MP analysis, all characteristics were equally weighted. MP trees were searched heuristically, with 1,000 random addition replications using tree bisection-reconnection (TBR) branch swapping. To assess confidence in clades, bootstrap tests were performed for NJ and MP trees using a full heuristic search. Replicates with TBR branch swapping were 10,000 and 1,000 for NJ and MP, respectively. For the ML tree, parameters were chosen based on the Akaike Information Criterion, as implemented in Modeltest ver 3.7 (Posada and Crandall 1998). The GTR + I + G model was selected for the sequence of COI and the combined sequence of COI and 18S rRNA genes. The Trn +I model was selected for the sequence of the 18S rRNA gene. ML trees were searched heuristically with TBR branch swapping. For the bootstrap test on ML, 1,000 replicates were performed using fast stepwise addition as a starting option. For Bayesian analysis, MrBayes ver. 3.2.0 (Ronquist and Huelsenbeck 2003) was performed with the GTR + I + G model selected by MrModeltest ver. 2 (Nylander 2004). The sample frequency and burn-in fraction were set to 100 and 25%, respectively. The numbers of generations were 20,000 for the COI gene, 1,220,000 for 18S rRNA gene, and 100,000 for the combined sequence of COI and the 18S rRNA genes.
Phylogenetic Pairwise Comparisons
The extent of ant association was categorized as either 0 (non-attendance) or 1 (facultative and obligate ant-attendance). Wolbachia infection rates in some species had a score of zero, and so were log-transformed. To examine the correlation between ant association (binary data) and Wolbachia infection rates in Tuberculatus species (continuous data), phylogenetically independent contrasts were calculated using pairwise comparisons (Read and Nee 1995; Maddison 2000) implemented in Mesquite ver. 3.61 (Maddison and Maddison 2019). First, taxa were paired as sets of terminal taxa (a set of comparison) with contrasting characteristics (ant-attendance vs. non-attendance). The pair of terminal taxa had to be phylogenetically separate from each other. Then, taxa pairs were counted up to maximal pairings that were equal to the number of parsimony-counted steps in the characteristic treated as unordered in a dichotomous tree. Finally, all maximal pairings were placed on a dichotomous tree, using an algorithm that performed two traversals through the tree. Wilcoxon matched-pairs signed-ranks test was applied on each pair to test whether the second characteristic (Wolbachia infection rate) had consistently higher or lower values to the first characteristic (ant-attendance and non-attendance).
Wolbachia supergroups
For Wolbachia that were detected in aphids (Table S2), the PCR products were sequenced with the same primers (16S-3f and 16S-2r) (Table 2). PCR products were purified with FastGene Gel/PCR Extraction Kit (Nippon Genetics, Tokyo, Japan). The cycle sequencing reaction was performed with a 5 µL volume consisted of 2 µL of Quick Start Mix (Beckman Coulter, Tokyo, Japan), 0.5 µL of 10 µM forward or reverse primers, and 2.5 µL of 10 ng/µL template DNA. The reaction cycle was 40 cycles of 94°C for 20 sec, 50°C for 20 sec, and 60°C for 1 min. DNA sequencing was analyzed using the CEQ2000XL DNA Analysis System (Beckman Coulter, Tokyo, Japan). The length of sequences that were successfully read through the samples were from about 500bp to 900bp. Multiple sequence alignment including the sequences of 16 Wolbachia supergroups (A, B, C, D, E, F, H, I, J, K, L, M, N, O, Q, S) that were cited by Bing et al (2014) (A to O), Ren et al (2020) (O found in aphids), Glowska et al (2015) (Q), and Lefoulon et al (2020) (S) (Table 4) was processed with the Clustal W (Thompson et al 1994) on the DDBJ. Supergroup P was not included in multiple sequence alignment because it was insufficient sequence length for the lower region of the gene. After multiple sequence alignment, the length of sequences was 471bp. Because of the insufficient length of sequences, it was expected that an unreliable model would be selected by the DNA base substitution model test. This study aims to determine what kind of Wolbachia supergroup is found in Tuberculatus aphids so that the unweighted pair group method with arithmetic mean (UPGMA) emphasized clustering formation was applied to construct Wolbachia phylogenetic tree. The bootstrap test was performed using the UPGMA method and 1000 replicates with TBR branch swapping.