Draft Aphaenogaster genomes expand our view of ant genome size variation across climate gradients

Sampling and whole-genome sequencing and assembly

Entire colonies of the six Aphaenogaster species were collected by A. Nguyen and C. Penick from field sites in eastern North America (Fig. 2). Ants were identified to species and voucher specimens have been deposited at the Museum of Comparative Zoology, Harvard University. Individuals from each colony were isolated from nest material and debris, weighed, placed in 50 ml Falcon centrifuge tubes, and immediately flash frozen in a −80 °C freezer. Colony weights were: 794 mg (Aphaenogaster ashmeadi), 652 mg (Aphaenogaster floridana), 520 mg (Aphaenogaster fulva), 749 mg (Aphaenogaster picea), 862 mg (Aphaenogaster miamiana), 280 mg (Aphaenogaster rudis 1), and 236 mg (Aphaenogaster rudis 2).

Whole colony DNA was used to have sufficient concentrations for sequencing. DNA was then extracted from each colony using methods developed previously for genomic sequencing of whole colonies of colonial mosquitos (Anopheles spp.) (Neafsey et al., 2010) and sequenced using an Illumina HiSeq 2500 at the Broad Institute (Cambridge, MA, USA). A combination of fragment and jump sequences were used to generate higher quality, long sequence reads.

Raw sequences were processed to remove chimeric and contaminant sequences, screened for contaminants by BLAST searches (using blastn) to identify sequences with likely matches to non-target species (primarily Wolbachia and Mycoplasma), and assembled using ALLPATHS-LG (version r48559) (Gnerre et al., 2011). Additional assembly processing using PILON (version 1.13) (Walker et al., 2014) was applied to reduce base-call errors and gaps in coverage. On average, across all seven genomes, PILON reduced coverage gaps by 3.1% or 3.9 Mb. GAEMR (http://www.broadinstitute.org/software/gaemr/) software produced summary statistics of the final assembled genomes. Once assembled, repeat regions in the Aphaenogaster genomes were detected and masked using Repeatmasker (version 4.0.5 Institute for Systems Biology). Genome completeness was assessed for all ant genome assemblies using BUSCO version 3.0.2 (Waterhouse et al., 2018) with the Hymenoptera BUSCO lineage (ODB9 from February 13, 2017) with the honey bee starting genome and HMMER version 3.1b1, BLAST+ version 2.2.31 and AUGUSTUS version 3.0.3 dependencies. For each genome we ran BUSCO on the final contaminant removed, repeat-masked Aphaenogaster assemblies and the most recent published assembly for the NCBI published genomes.

Analysis of genomes along climate gradients

After masking repeat regions, we applied MASH distance (Ondov et al., 2016) to measure pairwise dissimilarity of genomic sequences. The MASH method extends a data compression and dimensionality-reduction algorithm to generate estimates of sequence similarity with low computational overhead. Briefly, the pairs of genomic sequences were pre-processed into sets of k-mers of size 21 with the size of the non-redundant hashes retained set to 1,000. These settings have been demonstrated to provide good representation of genomic similarity with minimal computational costs (Ondov et al., 2016). These sets were then used to estimate the Jaccard similarity coefficient (the ratio of shared k-mers to total k-mers) of subsampled k-mer pairs of genomes. This unbiased estimate of the Jaccard similarity (J) was then used to calculate the dissimilarity of the two genomes (D) as D = 1 − J. All Jaccard similarity estimates had p-values less than 10⁻¹⁴, which is below the recommended 10⁻³ probability of observing values of J due to chance.

We used multivariate correlation analyses to examine biogeographic patterns of ant genomes. Mantel tests of multivariate correlation of distance matrices were used to examine correlations among ant genomes and climate variables. Specifically, we used directional (H_◦: Mantel r ≤ 0) partial mantel tests, which calculate the correlation between two distance matrices while controlling for the covariance with other matrices (Goslee & Urban, 2007). First, we examined the correlations between genomic similarity (MASH distance), whole-genome size similarity (Euclidean distance of assembly size in total base pairs) and climate variables (also using Euclidean distance). Via partial Mantel tests, we were able to isolate the correlation between genome size and climate by controlling for spatial autocorrelation and potential phylogenetic patterns by including geodesic and MASH distances as terms.

We obtained previously sequenced ant whole genome and climate data from publicly available databases. Whole genome sequences for ants were obtained from the NCBI Genome database (accessed August 2018, see Table S1). Climatic variables for each sampling location was obtained from the WorldClim database (version 2.0) at a 2.5 arc minute spatial resolution from the years 1970 to 2000 (Fick & Hijmans, 2017). Although used in the previous analyses of ant genomes, two species, (Wasmannia auropunctata and Monomorium pharaonis), which did not have published location information, were excluded from biogeographic analyses.

Using a permutational multivariate analysis of variance (PerMANOVA) procedure, we parsed the individual variables that were correlated with both genome size and MASH similarity. PerMANOVA is a flexible multivariate analog of ANOVA that permits the use of a wider set of similarity metrics to be used for the response matrix (Anderson, 2001), such as the MASH distance. We ran a total of 10,000 permutations of the original distance matrices for each statistical permutation procedure. We chose a subset of all possible climate variables available via WorldClim for this analysis. A visual inspection of the sampled climate variable correlations indicated that the primary climate variables, mean annual temperature (MAT), minimum temperature of the coldest month (Tmin), maximum temperature of the hottest month (Tmax), annual precipitation (PA) and precipitation seasonality (PS), represented the majority of climate variation (Fig. 3). Based on this, we only included these variables, along with latitude and longitude coordinates, as factors in the PerMANOVAs.

It is important to note that we are using assembly length as an indicator of genome size. As genome size estimates are generally used to set assembly size targets for whole genome sequencing efforts (Hare & Johnston, 2011), we expect there to be a high degree of correlation between assembly size and genome size. Also, as a test of the potential relationship between assembly size and true genome size, we examined the correlation between the average assembly sizes of ant genera that overlapped with flow cytometry estimates of those published in Tsutsui et al. (2008). We applied leave-one-out outlier detection (Cook & Weisberg, 1982) to least squares regression analysis examining the predictive potential of assembly length and flow cytometry size estimates at the genus level and found one strong outlier (Fig. S1), Dinoponera quadriceps, whose assembly size was less than half the flow cytometry estimate (259.7 vs. 554.7 Mb). After removing this species from the regression analysis we found a strong, significant relationship between assembly length and genome size (R² = 0.76, F_1,5 = 15.99, p-value = 0.010), which supports the use of assembly length as an indicator of genome size. Based on this, we did not include D. quadriceps in any of the analyses examining patterns of assembly length.

Data, computation and statistics

The raw and assembled genome sequences are currently archived at Harvard Forest (Petersham, MA, USA) and in NCBI’s genome database (Genome accessions NJRK00000000–NJRQ00000000 and BioSample accessions SAMN06894590–SAMN06894596). Genomic distance (MASH) computations were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. All analyses were conducted in R Core Team (2017). Analytical scripts for the project have been versioned and archived (DOI: 10.5281/zenodo.1341982) and are available online at https://zenodo.org/record/1341982. We used the vegan (Oksanen et al., 2016) and ecodist (Goslee & Urban, 2007) packages in R for multivariate analyses.

Genome quality and composition

DNA extractions yielded substantial amounts of high quality DNA with concentrations and quality scores ranging from 3.45–5.39 to 4.05–4.27 ngμL⁻¹, respectively. All genome assemblies displayed good coverage, with an average of 70% across all genomes (Table 3). Across all species, the length of the shortest contig at 50% of the genome (i.e., N50) was 18,864 bases; average assembly GC content was 38.18%; and average genome size was 370.45 Mb. Using GAEMR’s BLAST feature to conduct a search of the contigs against the NCBI’s nucleotide sequence database, we discovered that 38.98% and 22.04% of the top hits were “ant” and Aphaenogaster, respectively. Using the Hymenoptera BUSCO set, the Aphaenogaster genomes displayed a high degree of completeness. Between 96.1% and 97.6% with an average of 97.0% ± 0.2 (S.E.) of Hymenoptera BUSCOs were completely represented. The Aphaenogaster genomes were among the more complete genomes, relative to the other sequenced ant genomes from NCBI, which ranged from 88.2% to 98.5% complete Hymenoptera BUSCOs with a an average of 96.3% ± 0.5 (S.E.) (Fig. 4A). Also, the sizes of the Aphaenogaster genomes were within the range of other ant genomes based on size from both flow cytometry (Tsutsui et al., 2008) and the previously sequenced ant genomes available in NCBI (Fig. 4B).

Using the MASH genomic distances, we observed patterns of genomic similarity that were in line with expectations of ant relatedness. Sequences formed groups that corresponded with subfamilies (Fig. 5). Aphaenogaster clustered with other genera from the Myrmicinae and, in general, subfamily level clustering tended to follow previously observed patterns of subfamily relatedness (Bolton, 2006; Moreau et al., 2006; Ward, 2014). The Aphaenogaster sequences formed a single cluster containing only Aphaenogaster species and displayed intra-generic levels of genomic variance comparable to other genera (e.g., Trachymyrmex spp.). The separation of the two Aphaenogaster rudis species was initially surprising, as these two samples were collected at the same site (Duke Forest, Durham, NC, USA) and were identified as the same species based on their morphological characteristics (Ellison, 2012; DeMarco & Cognato, 2016). However, two recent studies of targeted gene regions have demonstrated the polyphyletic nature of Aphaenogaster rudis. One study of the evolution of the subfamily Myrmicinae observed that the genus as a whole could be split into at least four different lineages (Ward et al., 2015). Another, more detailed study of the genus in North America found that multiple individuals of Aphaenogaster rudis separated out into distinct groupings, each with other species, specifically, individuals of Aphaenogaster rudis from North Carolina (USA) were observed to form distinct clusters with individuals of Aphaenogaster carolinensis, Aphaenogaster miamiana, Aphaenogaster lamellidens and Aphaenogaster texana (DeMarco & Cognato, 2016).

Biogeographic patterns of ant genomes

After controlling for both spatial autocorrelation and potential phylogenetic patterns, we found a marginally significant, positive correlation between ant assembly size similarity and climate similarity (Mantel R = 0.12, p-value = 0.068). Although assembly size similarity and MASH genome similarity were not significantly correlated (p-value = 0.126), we included MASH as a covariate in addition to geodesic distance because previous research indicated that genome size is associated with phylogenetic relatedness (Alfsnes, Leinaas & Hessen, 2017). We found that different spatial and climatic variables were associated with the size similarity of ant genomes. Longitude but not latitude was a significant predictor of assembly size similarity (Table 4). Temperature of the coldest (Tmin) and hottest (Tmax) month and total PA, were significant predictors, but neither MAT nor PS were significant predictors of assembly size similarity. Overall, Tmin and PA had the strongest relationships with assembly size similarity with both exhibiting generally negative relationships with size (Fig. 6). When the newly sequenced Aphaenogaster genomes were excluded from the analysis, only PA was a significant predictor of assembly size similarity (R² = 0.31, p-value = 0.020 and see Table S2).