Patterns of genomic and allochronic strain divergence in the fall armyworm, Spodoptera frugiperda (J.E. Smith)

Abstract Speciation is the process through which reproductive isolation develops between distinct populations. Because this process takes time, speciation studies often necessarily examine populations within a species that are at various stages of divergence. The fall armyworm, Spodoptera frugiperda (J.E. Smith), is comprised of two strains (R = Rice & C = Corn) that serve as a novel system to explore population divergence in sympatry. Here, we use ddRADSeq data to show that fall armyworm strains in the field are largely genetically distinct, but some interstrain hybridization occurs. Although we detected F1 hybrids of both R‐ and C‐strain maternal origin, only hybrids with R‐strain mtDNA were found to contribute to subsequent generations, possibly indicating a unidirectional barrier to gene flow. Although these strains have been previously defined as “host plant‐associated,” we recovered an equal proportion of R‐ and C‐strain moths in fields dominated by C‐strain host plants. As an alternative to host‐associated divergence, we tested the hypothesis that differences in nightly activity patterns could account for reproductive isolation by genotyping temporally collected moths. Our data indicates that strains exhibit a significant shift in the timing of their nightly activities in the field. This divergence in phenology creates a prezygotic reproductive barrier that likely maintains the genetic isolation between strains. Thus, we conclude that it may be ecologically inaccurate to refer to the C‐ and R‐ strain as “host‐associated” and they should more appropriately be considered “allochronic strains.”

disruptive selection and population-specific directional selection can also act on sympatric populations resulting in divergence. For example, the availability of a novel host plant (Feder et al., 1988;Filchak et al., 2000;Rice, 1984), sexual selection resulting in assortative mating (Turner & Burrows, 1995), and divergent selection acting on phenology (Fukami et al., 2003;Santos et al., 2007Santos et al., , 2011 are all mechanisms that have been shown to lead to population divergence in sympatry. Because speciation takes many generations and is nearly impossible to study in real-time for multicellular organisms, researchers rely on populations that are undergoing various stages of speciation to gain insights into this process. The fall armyworm, Spodoptera frugiperda (J.E. Smith) (Lepidoptera: Noctuidae), is a moth species native to the Western Hemisphere. In the United States, this insect only overwinters in south Texas and south Florida, and these locations serve as the source for fall armyworms reinvading the northern US and Canada each year (Sparks, 1979). No evidence of a return migration has been found for this species. Thus, at the end of the season, individuals occurring north of these overwintering sites die and their genetic variation is lost (Nagoshi & Meagher, 2008). Importantly, this species is comprised of two morphologically identical but genetically distinct sympatric strains that have been previously described as "hostassociated". These strains were originally discovered by Pashley (1986) and named for the crops on which they were found, corn and rice (Pashley, 1986). Although this species is highly polyphagous feeding on up to 353 host plants from 76 plant families (Montezano et al., 2018), the larvae of the corn-strain, or C-strain, are considered to be primarily associated with corn and sorghum, whereas the rice-strain, or R-strain, are more commonly associated with pasture grasses, turf grass, and rice.
Although the two fall armyworm strains can be consistently identified using genetic markers, hybridization has been reported in both the lab and field, with field studies suggesting that 16% of moths collected show inconsistencies between multiple strain-specific genetic markers (Prowell et al., 2004). In these field assessments, the majority of putative hybrids have the maternally inherited mitochondrial markers from the R-strain, indicating a directional mating bias (Nagoshi, 2010;Prowell et al., 2004). In the lab, single pair matings have been conducted in both directions with evidence of reduced fertility among interstrain hybrids (Dumas et al., 2015;Kost et al., 2016). Still, these genetic markers have not become homogenous across the two strains, suggesting these strains remain genetically distinct despite occasional gene flow.
Despite being described as "host-associated", the host ranges of the two fall armyworm strains largely overlap, with evidence of asymmetric host use between strains (Groot et al., 2010. So, although the C-strain is more commonly found in association with corn and sorghum, the R-strain can also be found feeding on these hosts. It is uncommon, however, for the C-strain to be found feeding on smaller grasses such as pasture grasses and turf (Machado et al., 2008;Nagoshi, 2010). Several studies have used behavioral assays to assess whether the strains show a strong preference or fitness benefit when fed on different host plants, but no consistent differences have been found Pashley, 1988b;Pashley et al., 1995). This frequent overlap in habitat use and lack of strong host association suggests that differences in host plant use are unlikely to be the only factor maintaining genetic differentiation between these strains. Given the limited empirical evidence that host plant differences are involved in divergence, several studies have suggested that it may be more appropriate to refer to these strains as incipient species or genetic forms rather than host strains (Dumas et al., 2015;Juárez et al., 2014;Kergoat et al., 2012).
As an alternative to host-association, allochrony in nightly mating has been hypothesized as a mechanism underlying strain divergence. Pashley et al. (1992) found that strains exhibited differences in the timing of their mating activities, with the C-strain becoming active early in the scotophase (i.e., night) and the R-strain becoming active much later in the scotophase (Pashley et al., 1992). This temporal difference has been observed in subsequent laboratory mating assays Schöfl et al., 2009), and linked to heritable polymorphisms in the circadian rhythm modulator gene, vrille . If this temporal divergence also occurs in the field, the fall armyworm could be an excellent study system for assessing incipient allochronic divergence and speciation in action.
In this study, we used double digest restriction site-associated DNA sequencing data to better elucidate the patterns of divergence and gene flow between S. frugiperda strains collected from five locations across the central US. We then used temporal collection data to test the hypothesis that C-and R-strain fall armyworm moths in the field exhibit significant temporal differences in their nightly activity periods. Our study provides new insights into the patterns of genomic divergence and reproductive isolation of fall armyworm strains in the United States.

| Central US insect collection for sequencing
Spodoptera frugiperda moths were collected using universal moth traps baited with Scentry PSU 2-component lures (Scentry Biologicals) and containing Hercon Vaportape. Strains do not significantly differ in their attraction to pheromone lures, and therefore no bias in strain sampling is expected (Unbehend et al., 2014). Each trap was placed in or around corn and sorghum fields at five locations across the central US. Multiple sampling times were selected throughout the year that roughly corresponded to seasons when moths were present at each location (Table 1)

| DNA extraction
Prior to DNA extraction, the thorax was isolated from each specimen and surface sterilized in 95% ethanol. Tissues were tapped dry, placed individually in 2-ml Eppendorf tubes, and then frozen in liquid nitrogen. Sterilized plastic pestles were used to macerate the frozen thorax tissue. The Qiagen Gentra Puregene Tissue Kit was used to extract DNA following the manufacturer's protocol. The concentration of each DNA sample was measured on a NanoDrop spectrophotometer and all samples were diluted to a concentration of 50 ng/μl.

| Initial strain haplotype determination
After DNA extraction, strains were assigned using two known RFLPs in the cytochrome c oxidase subunit I (COI) mitochondrial gene (Levy et al., 2002;Nagoshi, Meagher, Adamczyk, et al., 2006). Briefly, the primer pair JM-76/JM-77 was used to amplify a 568 bp fragment of COI (Levy et al., 2002). Then 4 μl of the PCR product was added to both 2.5 μl of SacI (New England BioLabs) and 2.5 μl of MspI (New England BioLabs) diluted to their optimal working concentrations.
Reactions were incubated at 37°C for 1 h, and the products were run on a 1.8% agarose gel. The amplified C-strain mtDNA is cut once by MspI, and not by SacI while the R-strain mtDNA shows the reciprocal pattern (Nagoshi, Meagher, Adamczyk, et al., 2006). Based on the cut patterns of both restriction enzymes, each individual was assigned as having either a C-strain or an R-strain mitochondrial haplotype. After haplotype determination, DNA was stored at −20°C until sequencing.

| DNA sequencing, SNP calling, and filtering
DNA samples were sent to Texas A&M AgriLife Genomics and Bioinformatics Services (TxGen) for quality control, library preparation, and double digest restriction-site-associated DNA sequencing TA B L E 1 Collection location and date for all sequenced fall armyworm samples. The most common crop type surrounding each trap is listed as host plant. The numbers of individuals from each collection is given both as the number per predetermined strain mitochondrial haplotype (R-or C-) and the number of total individuals from each collection (R-+ C-) FastQ Screen v.0.14.0 with the BWA aligner was used to align raw reads to both the C-strain and R-strain published S. frugiperda genomes (Gouin et al., 2017). Sequences that did not match uniquely to one or both genomes were removed to clear the remaining sequences of all potential contaminant DNA (e.g., bacteria, pathogens).
Forward and reverse reads were then matched together using the repair function in BBMAP v.3.8.08 (Chaisson & Tesler, 2012). After filtering out contaminant DNA and DNA that matched multiple locations in the S. frugiperda genome, an average of 42.5% of the initial raw reads were retained for SNP calling.
Genomic loci that contained SNPs were identified using the dDocent v.2.2.16 SNP-calling pipeline (Puritz et al., 2014). In brief, dDocent removed low quality bases using Trimmomatic, and then mapped reads to the Liu et al. (2019) published chromosome map for S. frugiperda using BWA. The program FreeBayes then identified genomic loci containing SNPs and indels, and these variants were concatenated into a single VCF file. Our initial VCF file contained 441,437 variants.
Variants were filtered using VCFtools v.0.1.16 (Danecek et al., 2011). Specifically, all indels were removed and the remaining SNPs were filtered for a minimum PHRED score of 30. Only SNPs that were present in all individuals at a minimum of 3x coverage were kept in the final dataset. Finally, the dDocent_filters script (https:// github.com/jpuri tz/dDoce nt/blob/maste r/scrip ts/dDoce nt_filters) was run to complete SNP filtering. After filtering, the VCF file was manually examined and 236 SNPs did not map to a specific chromosome but rather to an 'unplaced_scaffold.' These unmapped SNPs were removed, leaving 5439 mapped SNPs in the final dataset.

| Analysis of molecular variance
The VCF file was uploaded into RStudio v.3.6.1 using R/vcfR v 1.9.0 (Knaus & Grünwald, 2017). To determine which collection factors (sampling year, location, season, and strain) significantly impacted the population structure of S. frugiperda, we carried out an analysis of molecular variance (AMOVA) using R/poppr v.2.8.3 (Kamvar et al., 2014). We used a Monte Carlo test with 1000 random permutations to determine the statistical significance of each factor in the AMOVA.

| Population structure and strain admixture analyses
To examine the population structure within fall armyworm samples, the VCF file was converted to a biallelic.bed (Plink Binary Biallelic Genotype Table) file and then to Eigenstrat format using PLINK v.1.07 (Purcell et al., 2007) and EIGENSOFT v.7.2.1 , respectively. We then used the smartpca function within EIGENSOFT 7.2.1 to conduct a smart principal component analysis that identified and removed outliers in the dataset caused by population stratification Price et al., 2006). This program calculated Tracy-Widom statistics to determine the number of significant eigenvalues, or principal components, within the PCA. The PCA results were then plotted using R/ggplot2 (Wickham, 2016). Because the smartPCA revealed two distinct SNP based population clusters that corresponded to host strains, putative hybrids were identified as individuals that either did not fit neatly into the R-strain or C-strain SNP clusters, or individuals that had a mismatch between their mtDNA and SNP-based strain assignments.
To determine if the putative hybrids showed significant evidence of interstrain admixture, outgroup f3 statistics were calculated using the 3-populations test function (qp3Pop) in AdmixTools v.5.0 (Patterson et al., 2012). In this test, pure-strain individuals from the R-strain and the C-strain were defined as the ancestral populations, and each putative hybrid individual was uniquely assessed for admixture using the model f3(C-strain, R-strain; putative hybrid individual). Only individuals that had significantly negative f3 values were considered true hybrids.
The program ADMIXTURE 1.3.0 was run using default parameters to determine the probability that individual moths were assigned to one or more genetically distinct groupings (Alexander & Lange, 2011). The K-values input ranged from 1 to 15, and the optimal value of K was determined as the run that resulted in the lowest cross validation (CV) error. The browser-based program CLUMPAK was used to visualize the population assignment of all individuals (Kopelman et al., 2015).

| Assessing genomic patterns of strain divergence
Using the smartPCA and admixture results, each individual was assigned to the R-strain or the C-strain based on their SNP group-  (Tessnow et al., 2021).
After strain determination, generalized linear models with a binomial distribution and logit link were run to determine the effects of col-

| Factors that contribute to genetic variation in sequencing data
Prior to sequencing, we determined the relative proportion of each central US collection that was comprised of the R-strain and C-strain individuals using mtDNA haplotypes. Although traps were placed F I G U R E 1 Proportion of individuals with C-strain (blue) or R-strain (orange) mtDNA in each moth collection. The dark grey overlay illustrates the proposed central US flyway for fall armyworm moths (inferred from Westbrook et al. (2016)) in and around corn and sorghum fields which are typically considered C-strain host plants (Pashley, 1988a(Pashley, , 1988b, we found that most locations contained a mix of both C-and R-strain haplotypes. Unexpectedly, several collections, especially during the fall season, were solely comprised of individuals with R-strain haplotypes ( Figure 1). Because we had collected a representative sample of both host strains during most collection times, a mix of individuals comprising both host strains were sequenced (Table 1).
Although there were several potential sources of genetic structure in our dataset (year, location, season, and strain), our AMOVA revealed predetermined host strain haplotype was the only factor that significantly accounted for genetic variability in the data ( were not significant ( Table 2).
The major effect of host strain haplotype on genetic structure was further supported by a smartPCA conducted on the SNP dataset.
In this analysis, the data clustered into two distinct groupings along PC1 that roughly corresponded with the predetermined mitochondrial haplotypes (Figure 2a). Other principal components did not identify any additional population groupings. When conducting the SmartPCA, 33 samples were removed as outliers due to cryptic relationships (genetically too similar within collections). Roughly half of these outliers were from the spring collections conducted in the Lower Rio Grande Valley, in which caterpillars were collected as opposed to moths. Since, it is reasonable that some sibs or half sibs were collected when sampling caterpillars from the same fields, we continued our analysis with only the remaining 393 unrelated individuals.
In addition to conducting a PCA, Tracy-Widom statistics were calculated to evaluate the statistical significance of each eigenvalue or principal component . We found that only the Because multiple lines of evidence indicated that S. frugiperda population structure across the central US was explained entirely by strain, further analyses focused on admixture and genetic differentiation between strains.

| Admixture between host strains
Although the two clusters in the PCA could largely be explained by strain mtDNA haplotypes, there were 34 individuals that showed a mismatch between the mtDNA haplotypes and their SNP cluster assignment. One hundred percent of these mismatched individuals contained R-strain mtDNA, but clustered within the C-strain SNP cluster. Additionally, two individuals with R-strain mtDNA and three individuals with C-strain mtDNA mapped directly in between the two SNP clusters on the smartPCA (Figure 2a). This assemblage of 39 individuals comprised of mismatches and those that did not group were classified as putative hybrids.
f3 statistics were calculated to determine if these putative hybrids exhibited significant admixture between the two strain source populations. In this analysis, f3 values significantly lower than zero indicate significant admixture between two source populations. We found that all five individuals that did not neatly group with either SNP cluster exhibited significant admixture between the two host strains.
Additionally, two individuals that grouped with the C-strain SNP cluster but carried the R-strain mtDNA (mismatches), exhibited significant strain admixture (Figure 2b). These seven individuals will henceforth be referred to as the hybrids and are indicated as such on Figure 2a.
ADMIXTURE analysis was then conducted to calculate the probability of individuals assigned to one or more (K) genotypic groups.
The lowest cross-validation error (CV) occurred when K = 2 or with two genotypic groups. K = 3 or higher did not show any additional population resolution. These two genotypes largely corresponded with the predetermined strain haplotypes (Figure 2c). Slightly more admixture was detected among individuals with the R-strain mtDNA compared to those with C-strain mtDNA.

| Genomic patterns of strain divergence
To determine the level of divergence between strains at every SNP

| Temporal differences in strain activity
In our temporal collection experiments, a total of 156 moths were  When the data were combined across both locations, we found that both factors of time (χ 2 = 46.144, df = 2, p < .001) and location (χ 2 = 20.784, df = 1, p < .001) had a significant effect on the probability of collecting C-or R-strain individuals. However, the interaction of time x location was not significant (χ 2 = 0.702, df = 2, p = .704). This indicates that although the relative proportions of each strain differed between the two locations (College Station was 82% R-strain and 18% C-strain while Lubbock was 5% R-strain and 95% C-strain), the timing of activity between strains remains similar.

| DISCUSS ION
The fall armyworm is an excellent system to study how reproductive isolation may be maintained in two sympatric populations or strains.
Although these strains are capable of intermating, molecular markers have not become homogenous across the strains since their identification several decades ago (Levy et al., 2002;Nagoshi, Meagher, Adamczyk, et al., 2006), suggesting they remain genetically distinct in the field despite occasional gene flow. Here, we used ddRADSeq data and temporal collection data to address two primary questions about fall armyworm strain divergence in the field: (1) To what degree do strains exhibit reproductive isolation? and (2) Do differences in daily phenology significantly contribute to reproductive isolation between strains in the field?

| Genomic evidence of reproductive isolation
Here, we present evidence for two genetically distinct populations of fall armyworms in the central US that correlate with the previously described host-associated strains (Pashley, 1986). Although our data clearly indicates these strains are genetically distinct, there was evidence for occasional interbreeding. Using f3 tests, we identified seven samples that exhibited significant interstrain admixture. These samples included five putative F1 hybrids that mapped neatly in the middle of the two host strains on the PCA, and two individuals that had a mismatch between their mtDNA strain haplotype and their SNP strain genotype. This is strong evidence that hybridization occurs between the two strains in the field. Since both R-strain and C-strain mtDNA was recovered among the putative F1 hybrids, we conclude that in the field, fall armyworm females of both the C-strain and the R-strain occasionally mate with males of the opposite strain. Since two mismatch individuals also showed significant signs of admixture but appeared to be closer related to the C-strain than the R-strain, these individuals are predicted to be the offspring of a hybrid female of R-strain maternal origin backcrossed to a C-strain male.
Thirty-two additional individuals showed a mismatch between the maternally inherited mtDNA strain marker and their SNP genotype.
Although these individuals were not significantly admixed according to our f3 test, we suspect they are the result of past hybridization events followed by several generations of backcrossing. Overtime, the signal of admixture has been reduced and can no longer be detected using f3 statistics, however, the mismatch between the mtDNA and SNP genotype is still evident. Interestingly, 100% of moths that have a mismatch between their maternally inherited mtDNA and their SNP genotype have the R-strain mtDNA markers. Previous studies that defined hybrids as individuals with a mismatch between the mtDNA and nuclear strain markers also found the majority of putative hybrids collected in the United States have R-strain maternal origin (Nagoshi, 2010;Nagoshi et al., 2017;Nagoshi & Meagher, 2003;Prowell et al., 2004). This indicates that in the field hybrid females with R-strain maternal origin are successfully mating with C-strain males. However, hybrids with C-strain maternal origin are not backcrossing to the Rstrain. This is strong evidence that a unidirectional barrier to reproduction exists, limiting gene flow between these strains in the field.
Although this pattern of unidirectional introgression has been routinely recovered in field data, laboratory assays have been less consistent. Some studies have found a unidirectional mating bias where R-strain females are able to mate and produce offspring with C-strain males, but the reverse is not true (Pashley & Martin, 1987), while other laboratory assays successfully conducted reciprocal crosses of both strains (Quisenberry, 1991;Whitford et al., 1988) and have even found that hybrid females of C-strain maternal origin are more fertile than hybrids of R-strain maternal origin . Still another study successfully conducted reciprocal interstrain crosses, but found that F1 hybrids of C-strain maternal origin had a drastic reduction in fitness, while F1 hybrids with R-strain maternal origin had only a minor fitness cost (Dumas et al., 2015). This last study by Dumas et al. (2015) is most consistent with our field data. Although it is not clear why variable results have been found in laboratory mating assays, these studies have consistently reported behavioral and/or physiological barriers that limit hybridization between strains, serve as barriers to gene flow, and can reinforce strain identity.
Although we detected hybridization between strains, mitochondrial markers have been used as reliable strain indicators across multiple regions of the Western Hemisphere for the past several decades Levy et al., 2002;Nagoshi, Meagher, Adamczyk, et al., 2006;Nagoshi et al., 2007). If gene flow occurs between strains and causes a mismatch between the mtDNA and SNP genotype, then these markers would be expected to homogenize across strains and become less reliable overtime. Since this has not been observed, we hypothesize that the combination of selection against unfit hybrids, a one-way migration that removes admixture occurring north of the overwintering site, and genetic drift caused by large population size fluctuations at the overwintering site (Nagoshi & Meagher, 2004b), play a role in maintaining the genetic integrity of these strains despite occasional interstrain gene flow.

| Temporal divergence between strains
Despite collecting fall armyworm moths in fields dominated by corn strain host plants, our initial trap captures from across the central US generally comprised both R-strain and C-strain individuals. This coexistence at the same location and time is consistent with previously reported trap captures across the United States Nagoshi & Meagher, 2004a), and indicates that male moths of both the C-and R-strain routinely follow female pheromone trails to search for mates in the same fields. Although fall armyworm strains have been previously described as host plantassociated, both strains have exhibited equal preference and performance when fed on a variety of crop types (Groot et al., 2010;Pashley, 1988b;Pashley et al., 1995). This frequent overlap in habitat and mating locations combined with a lack of evidence to support strong host-plant associations, suggests that differences in host plant use are unlikely to be the primary factor maintaining genetic differentiation between these strains.
Another factor that has been implicated in the genetic divergence of these strains is allochronic differences in mating time Pashley et al., 1992;Schöfl et al., 2009). After observing the mating behavior of C-and R-strain fall armyworm colonies in the lab, Pashley et al. (1992) found that they exhibited strong differences in the timing of their mating activities, with the C-strain becoming active early in the scotophase (i.e., night) and the R-strain becoming active much later in the scotophase (Pashley et al., 1992).
This temporal difference has been observed in subsequent laboratory mating assays Schöfl et al., 2009), and linked to heritable polymorphisms in the circadian rhythm modulator gene, vrille . However, this difference has never been investigated in the field. Therefore, we used temporal field collections to experimentally test the hypothesis that C-and Rstrain fall armyworm moths in the field exhibit significant temporal differences in their activity periods.
We found strong evidence across both College Station, TX and Lubbock, TX that strains exhibit differences in the timing of their nightly mating activities in the field. C-strain moths were collected in pheromone traps early in the night (0-5 h after sunset), while the R-strain was only active later in the night more than 7 h after sunset Pashley et al., 1992;Schöfl et al., 2009). Although C-and R-strain moths were both collected at each of our sampling locations, the population composition differed. In College Station, TX, the majority of the population were of the R-strain while in Lubbock, TX, most of the population was made up of the C-strain. Despite these differences, there was no interaction between location and time. This is strong evidence that regardless of the number of C-or R-strain individuals present in a population, the temporal differences in activity between strains persist. This clear shift in daily phenology between strains in wild fall armyworm populations confirms previous lab observations and provides a critical pre-zygotic isolating mechanism that likely reduces interstrain hybridization while these moths are residing in the same fields.
Although all moths collected in the pheromone traps in this study were males, laboratory studies have indicated that male nightly activity periods are more labile than those of females (Pashley et al., 1992). Thus, it is likely that females are exhibiting similar if not stronger differences in activity patterns in the field.
Our data strongly support allochrony as a critical prezygotic reproductive isolating mechanism underlying fall armyworm strain differentiation in the field. In light of our data, there are now several lines of evidence to suggest that referring to fall armyworm strains as 'host-associated' may be ecologically inaccurate. First, behavioral studies have failed to consistently find differences in host plant preference and performance between the strains (Groot et al., 2010;Juárez et al., 2014;Pashley, 1988b;Pashley et al., 1995). Second, clear phenological differences have been routinely observed in the lab Pashley et al., 1992;Schöfl et al., 2009). And lastly, differences in nightly activity patters between the two strains have now been confirmed in the field. Thus, the combined evidence indicates that it may be more accurate to refer to these two genetically distinct fall armyworm populations as allochronic strains.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All raw sequence reads have been made available through the NCBI sequence read archive (SRA). These sequences can be retrieved using the bioproject accession numbers PRJNA645462 or can be accessed at the following link https://www.ncbi.nlm.nih.gov/biopr oject/ ?term=PRJNA 645462. Additionally, all information about our SNP dataset and the moths collected to assess temporal differences in strain activity has been made available in the Dryad Digital repository and can be found using the following DOI accession number: 10.5061/dryad.4qrfj6qc1.
O RCI D TA B L E A 2 Chromosome position and NCBI BLAST matches for sequences containing SNPs that appeared as outliers (high F st values) on the Manhattan plot (Figure 3a). Many of the SNPs that showed high levels of divergence were clustered into six groups.