Epigenetic variation allows for rapid changes in phenotypes without alterations to nucleotide sequences. These epigenetic signatures may diverge over time among isolated populations. Epigenetic incompatibility following secondary contact between these populations could result in the evolution of reproductive isolating mechanisms. If epigenetic incompatibility drove the evolution of species isolating mechanisms, we expect to see significant epigenetic differentiation between these species. Alternatively, epigenetic variation could be the result of predominantly environmental variables and not align along species boundaries. A methylation sensitive amplified fragment length polymorphism analysis was performed on individuals of the closely related katydid species Neoconocephalus robustus and N. bivocatus. We observed significant variation in total methylation levels between species. However, genetic differentiation remained larger than epigenetic differentiation between species groups. We measured a significant correlation between the epigenetic and genetic distance between individuals. Epigenetic differentiation is therefore likely the result of an interaction between genetic and epigenetic loci and not a mechanism for species differentiation. We therefore did not find evidence to support our hypothesis of an epigenetically mediated mechanism for speciation between N. robustus and N. bivocatus.

genetic differentiation, methylation, MSAP, Neoconocephalus, population genetics

Epigenetic variation can lead to changes in phenotypic expression without any change to the nucleic acid sequence (Berger et al. 2009). Unlike genetic variation, epigenetic changes can produce reversible, heritable phenotypic changes within a lineage (Jablonka and Lamb 1998). Epigenetically controlled phenotypic plasticity may therefore play an important role in rapid, adaptive changes (Ledón-Rettig et al. 2012, Burggren and Crews 2014). DNA methylation, the most commonly studied form of epigenetic modification, involves the addition of methyl groups, usually to CpG dinucleotides that regulate gene expression (Boyes and Bird 1991). Methylation regulated phenotypic plasticity has been documented across multiple taxonomic groups (Braam and Davis 1990, Rapp and Wendel 2005, Varriale and Bernardi 2006, Elango et al. 2009, Herrera and Bazaga 2010). Until recently, epigenetically dependent variation has been largely overlooked as an evolutionary mechanism contributing to speciation (Smith et al. 2016). Facultative phenotypes can be maintained in a population as adaptive alternatives to divergent environmental conditions. These alternative phenotypes may ultimately become fixed in a population by way of genetic assimilation (Pál and Miklós 1999, West-Eberhard 2005).

Phenotypic variation can also evolve in response to epigenetic incompatibility (= hybrid incompatibility caused by epialleles; Bente and Scheid 2017, Blevins et al. 2017) between groups. Following the epigenetic diversification of groups in isolation, intermediate forms may show reduced fitness following secondary contact (Jablonka and Lamb 1999). Aphids moved to a new host plant and allowed to reproduce asexually quickly developed a preference for the new host and produced morphological characteristics similar to conspecifics utilizing this host. Backcrossing to the parental line resulted in the production of non-viable offspring (Shaposhnikov 1966). Heritable phenotypic variation in aphids has been correlated with changes in methylation (Field et al. 1989, Walsh et al. 2010). Methylation dependent incompatibility could be reinforced by reproductive isolation (Pál and Miklós 1999).

In many animals, behavioral isolation plays a significant role in maintaining species boundaries. Acoustic communication has been studied as a mechanism for reproductive isolation in both vertebrate and invertebrate groups (Gerhardt and Huber 2002). Neoconocephalus, a New World genus of katydids, possesses a diverse range of habitat preferences (Walker et al. 1973, Frederick 2013) and call phenotypes (Schul and Patterson 2003, Beckers and Schul 2008, Deily and Schul 2009). Females utilize these calls in mate recognition and in phonotaxis, the directional movement towards a calling male (Greenfield 1990, Triblehorn and Schul 2009, Beckers and Schul 2010). The acoustic communication system of Neoconocephalus allows for reproductive isolation among the multiple species that may be found in sympatry (Schul et al. 2014).

The sibling species (Snyder et al. 2009) N. robustus (Scudder, 1862) and N. bivocatus Walker, Whitesell, & Alexander, 1973 differ both in their habitat preferences and call types. Neoconocephalus robustus is a grassland generalist and utilizes grasslands with a wide range of flora and environmental conditions. Neoconocephalus bivocatus, on the other hand, is a prairie obligate preferring drier habitats composed primarily of tall prairie grasses (Walker et al. 1973). Both species co-occur in some areas, e.g. at the edges of prairie remnants or along creeks within prairies. This species pair is largely morphologically cryptic but maintains divergent call characters and call preferences (Deily and Schul 2004, 2006). While the two species can be distinguished using genetic markers (e.g. mitochondrial haplotypes or AFLP markers (Snyder et al. 2009, Ney and Schul 2017)), genetic differentiation is weak and evidence for hybridization (e.g. intermediate calls) occur frequently. Population-level genetic variation was shown to be very low, with almost no geographically associated genetic variation between sites as distant as 448 km apart (Ney and Schul 2017). This suggests that these two species diverged very recently, most likely during the current glacial cycle (Frederick 2013).

A likely scenario of this divergence event includes an ancestral population living in both habitat ranges. In each habitat, different epigenetic patterns would be expressed, which then became fixed within each subpopulation, resulting in decreased hybrid fitness. Genetic differences leading to the differences in male calls and female preferences evolved later, in part driven by reproductive reinforcement (Deily and Schul 2004). Thus we hypothesize here that epigenetic differentiation, driven by the different environmental conditions experienced in their core habitats, drove divergence of the species. The lack of substantial genetic differentiation between N. robustus and N. bivocatus could point to an underlying epigenetic incompatibility driving species differentiation. This hypothesis predicts that the epigenetic differentiation between the species would be larger and more stereotyped than the genetic differences.

DNA methylation has been described in many insect groups including multiple orthopteran species (Sarkar et al. 1992, Robinson et al. 2011). Methylation sensitive amplified fragment length polymorphism (MS-AFLP) analysis is one of the few techniques that allows for the quantification of genome-wide patterns of cytosine methylation without any previous knowledge of genome sequences (Xiong et al. 1999). The MS-AFLP technique uses two isoschizomer restriction enzymes (MspI and HpaII) that recognize the same restriction site, cutting in the same location (5’-C^CGG-3’), but have different sensitivities to the presence of cytosine methylation. By comparing the presence/absence of fragments produced by both enzymes, the methylation status at each restriction site can be evaluated. In addition, genetic polymorphisms can be evaluated between individuals that lack the restriction site in both digestions as an indication of a change in nucleotide sequence.

To distinguish whether divergence of these two species was initiated by epigenetic, rather than genetic variation, we addressed two questions using a MS-AFLP technique. First, we asked to what extent epigenetic and/or genetic differentiation predicted species assignment. Second, we analyzed whether the patterns of genome-wide DNA methylation were correlated with genetic variation.

Specimen collection.—We utilized a total of 94 males collected in the summers of 2006, 2013, and 2014 from eight grassland sites around the state of Missouri (Suppl. material 1: Table S1). We used males’ calls to localize individuals in the field and collected them by hand after sunset. We identified males in the field as belonging to members of the target group through their call and morphological features, including cone pigmentation and body size, prior to collection (Walker et al. 1973). As the two species are morphologically cryptic, we identified individuals as either N. bivocatus or N. robustus based on their divergent call characteristics as described below. Individuals were collected from among eight grassland sites within Missouri and from within three different years (Suppl. material 1: Table S1). Thus, a range of epigenetic patterns, expressed within variable environmental conditions, was included in each species sample. Fixed differences observed between the two species are therefore likely to be species specific rather than the result of differences in environmentally mediated methylation patterns.

Call recordings.—We recorded male calls in 2006 and 2013 within three days of collection using an Audiotechnica ATR 55 microphone and a Marantz PMD-671 solid-state recorder (16 bit, 48 kHz sampling rate). Recordings were made outdoors with males placed in individual mesh cages (approximately 10×20×10 cm) spaced at least 3 m apart. In 2014, we recorded male calls in the field immediately preceding collection using a Tascam DR-40 linear PCM recorder (16 bit, 48 kHz sampling rate). Ambient temperatures ranged from 22–28°C.

Temporal call analysis.—Call temporal patterns were analyzed and species assignments made as described in Ney and Schul (2017). In short, we marked each sound pulse produced during a single closing movement of the wings (Walker 1975) using custom software. The recordings were rectified and the envelope extracted with a temporal resolution of 0.125 ms. Pulses were detected automatically using a threshold-based algorithm and manually checked before saving the data to text-files. We analyzed about 2 s of each male’s call containing 150–250 pulses. Further analyses were conducted on the text files in MS Excel.

Male calls of N. robustus have a single pulse rate of about 200/s, equivalent to a pulse period of about 5 ms (Walker 1975). Females recognize this pattern by the absence of silent gaps longer than about 2 ms (Deily and Schul 2004). In contrast, male calls of N. bivocatus have two alternating pulse periods of about 4.5 ms and 7 ms (Deily and Schul 2004), resulting in a ‘galloping rhythm’ or ‘pulse pairs’. Female N. bivocatus recognize this pattern by the rate of the pulse pairs: pulse pair rates around 87/s are attractive and attractiveness decreases towards faster and slower rates. For the pulse pairs to be detectable, the two alternating pulse periods must differ sufficiently (J. Schul, unpublished work).

We quantified the ratio of the means of the alternating pulse periods (longer pp / shorter pp). For the single pulse calls of N. robustus, this should result in values close to one, while in N. bivocatus significantly larger values result (Bush and Schul 2010). Among the calls of 94 males analyzed, there were two distinct groups of values, one between 1 and 1.18 and a second group with ratios >1.38. A natural break occurred in the data between 1.18 and 1.38. We classified ratios <1.18 as ‘N. robustus’, and >1.38 as ‘N. bivocatus’ in call type.

Spectral call analysis.—Call spectra were analyzed and species assignments made as described in Ney and Schul (2017). In short, we analyzed the amplitude spectrum of male calls using a fast Fourier transformation (FFT, Hamming window, frame length 256) in Audacity v.1.3 (Audacity Team 2008). Spectra were averaged over 2 s of the call. The main energy in the spectrum of Neoconocephalus calls is concentrated in a narrow band between 7 and 15 kHz with the peak frequency differing among species (Schul and Patterson 2003). We measured the center frequency of this low frequency band as the geometric mean of the upper and lower cut-off frequencies at -3 dB from the peak amplitude.

In N. bivocatus, the center frequencies of this band have a mean of about 10 kHz among individuals, ranging from 7 to 15 kHz. Females have little spectral selectivity in this frequency range (Deily and Schul 2006). The center frequencies of calls with a pulse period ratio of >1.38 that also fell into this range of center frequencies were classified as N. bivocatus. One individual had a pulse period ratio greater than 1.38 and a center frequency lower than 7 kHz; this male was classified as an intermediate caller, as its call characteristics differed from that of both species.

In N. robustus, the low frequency band is narrower than in N. bivocatus and is typically limited to 10 kHz and below (Schul and Patterson 2003). Indeed, frequencies well above 10 kHz have an inhibitory effect on female phonotaxis (Deily and Schul 2006). Of the individuals with pulse period ratios <1.18, center frequencies of all but five individuals clustered between 6 and 9 kHz. The five remaining individuals had center frequencies ranging from 9.7 to 11.1 kHz, which would be less attractive for N. robustus females. These individuals were classified as possessing an intermediate call type. All intermediate callers were removed from the analysis of epigenetic and genetic differentiation.

Molecular analysis.—We aim here to analyze the differences of trans-generational methylation patterns between two species. Therefore, we selected tissue that is likely to have little tissue-specific methylation differences between the two species. We collected tissue from hind femurs (mostly muscle and cuticle) that should be under similar selection in both species and thus little differential tissue-specific methylation. Using leg tissue also avoids the risk of contamination from GI tract microbiota.

We removed the hind femurs of collected males and placed them in 95% EtOH for DNA preservation. We later extracted DNA from the hind femurs using the DNeasy Blood & Tissue Kit (Qiagen Inc., Valencia, CA, USA). DNA quantification was performed on each sample by spectrophotometry (NanoDrop 1000, Thermo Scientific, Wilmington, DE). Genomic DNA was stored at -80°C prior to molecular analysis.

We used a MS-AFLP assay modified from Xu et al. (2000). DNA (55 ng) from each sample was digested in two separate double digest reactions (EcoRI/HpaII and EcoRI/MspI). EcoRI selectively cuts the sequence 5’-GAATTC-3’. HpaII and MspI are isoschizomers, meaning they selectively cut at the sequence (5’-CCGG-3’), but differ in their sensitivity to cytosine methylation at those sites. HpaII will only cut if the external cytosine is hemimethylated. MspI cuts when the internal cytosine is either hemi or fully methylated. Both enzymes will cut if the target site is completely unmethylated. Using this method, we evaluated the CpG methylation of restriction sites by comparing the fragments produced by both digests.

Digestion and ligation were carried out together to prevent regeneration of restriction sites. Synthetic double stranded DNA adaptors (Xu et al. 2000) were ligated to the cleaved ends of restriction sites. The EcoRI/HpaII and EcoRI/MspI digestion/ligation reaction (11 ul final volume) is comprised of 1.1 μL 10X CutSmart^TM buffer (New England Biolabs, USA), 0.55 μg/μL of Bovine Serum Albumin (BSA) (New England Biolabs), 0.3 μL water, 5 U of EcoRI HF, 1 U of either HpaII or MspI, 1 U T4 ligase (New England Biolabs), 1 μL ATP (10 μM), EcoRI adaptors (5 μM), either HpaII or MspI adaptors (50 μM), and 5.5 μL genomic DNA (10 ng/ μL). The reaction was incubated at 37°C for 2 hours. Preselective PCR was conducted with 1:10 dilute digestion/ligation products, the Eco+A (5’-GACTGCGTACCAATTCA-3’), and the HpaII/MspI+A (5’-GATGAGTCTAGAACGGA-3’) primers using thermocycler settings as described in Snyder et al. (2009). We performed selective PCR independently with three primer pairs (Table 1) and 1:100 dilute preselective PCR products from both the HpaII and MspI reactions. Fluorescently labeled Eco primers (Eco+AAC (6FAM), Eco+AGC (PET)) were used in selective PCR (as described in Snyder et al. 2009) and the products multiplexed and diluted to produce a 1:10 final dilution of each product. Fragments were separated in an ABI 3730 genetic analyzer at the DNA Core Facility, University of Missouri. We called MS-AFLP bands using GeneMarker v.1.6 (Hulce et al. 2011) using an automated peak-calling scheme (as described in Holland et al. 2008) and called alleles between 75–500 bp with a minimum peak intensity of 50.

Data analysis.—We obtained presence/absence fragment data for both EcoRI/HpaII and EcoRI/MspI datasets from GeneMarker (Hulce et al. 2011). Presence of a fragment in both the EcoRI/HpaII and EcoRI/MspI digestions (1/1) was described as unmethylated. The presence of a fragment in one digestion but not the other (0/1 or 1/0) was defined as methylated (either hemimethylated or internal cytosine methylated, respectively). If fragments were absent in both digestions (0/0) the loci were considered uninformative as this state could be due to either full methylation at the target site (methylation of both the inner and outer cytosine) or the absence of the fragment due to variation in the nucleotide sequence between individuals (Schulz et al. 2013, Fulneček and Kovařík 2014).

The methylation sensitivity of each locus was identified using the MSAP Package (Pérez-Figueroa 2013) in R v.3.1.2 (R Core Team 2015). Loci were classified as methylation sensitive loci (MSL) if there was evidence of methylation in at least 5% of the sampled individuals at that locus. Genetic markers were extracted from the MS-AFLP loci. Fragments that were present in one or both MspI and HpaII analyses were scored as present. Fragments that were absent in both reactions were scored as absent. This method of using MS-AFLP loci to estimate genetic parameters has been found to produce similar results to standard AFLP markers (Smith et al. 2016). Only polymorphic MSL and genetic loci that differed among sampled individuals were used in further analyses. We evaluated variation for each epigenetic (MSL) and genetic locus individually then calculated the mean diversity of MSL and genetic loci using Shannon’s diversity index (I). We compared the relative frequency of total methylation (internal cytosine methylated and hemimethylated fragments) and unmethylated MSL between species. Significant variation between species was tested using a Mann-Whitney U test.

We tested for differences in MSL and genetic differentiation using two-way analyses of molecular variance (AMOVA; Excoffier et al. 1992) that grouped individuals again by species. Significance of the test statistic’s (Φ_ST) deviation from zero was estimated based on 1000 permutations of individuals among groups. Principal coordinate analyses (PCoA) of both epigenetic and genetic loci were performed using the R stats package v.3.1.2 (cmdscale; R Core Team 2015) to visualize the Euclidean distance between species. A population structure analysis was implemented for all individuals using the program STRUCTURE v.2.3.3 (Pritchard et al. 2010). The admixture model was used, allele frequencies correlated, with a run length of 100,000 (Burnin = 10,000) for 10 replicates each of K = 1–10 (genetic clusters). The estimate of the most well supported K was calculated as described in Evanno et al. (2005) and implemented in Structure Harvester v.0.6.94 (Earl 2012). The program Clumpp v.1.1.2 (Jakobsson and Rosenberg 2007) was used to align the 10 repetitions of the most well supported number of clusters.

If similar signals exist for the MSL and genetic structure then epigenetic and genetic distance may be significantly correlated. We estimated the Euclidean distance between individuals for both epigenetic and genetic datasets using the R stats package v.3.1.2 (R Core Team 2015) and compared the distance matrices using a Mantel test (Mantel 1967; ade4 package v.1.7.2; Dray et al. 2007). We utilized 10,000 permutations of the design matrix to determine the significance of the correlation coefficient.

Genome-wide variation in methylation.—We included 94 males in our analysis. We determined the males’ species assignments based on the temporal and spectral frequency call preferences of N. robustus and N. bivocatus. The call analysis classified the 94 males as 31 N. bivocatus, 57 N. robustus, and 6 having an intermediate call type; five of the intermediates had the N. robustus temporal pattern and one the N. bivocatus pattern, with frequency spectra that fell outside of these respective species’ specific pattern (Fig. 1). Individuals possessing intermediate call phenotypes were removed from further analyses. The MS-AFLP analysis yielded 227, 277, and 318 fragments in each of three selective PCR reactions (Table 1). Of these fragments, a total of 364 were polymorphic for their methylation status among sampled individuals and were classified as MSL. A loci-calling scheme, utilizing the MS-AFLP dataset, allowed for the inference of the genetic state of fragments among individuals, similar to a traditional AFLP analysis. In total, 668 polymorphic genetic loci were produced from 822 MS-AFLP fragments.

Figure 1. 

Species assignment based on call pulse period ratio and center frequency. Labeled boxes indicate the calls classified as N. robustus and N. bivocatus. Individuals that fall outside of species classifications were removed from further epigenetic and genetic analyses (as described in Ney and Schul 2017).

We compared whether species differed significantly in their proportion of genome-wide methylated sites. The relative frequency of genome-wide methylation (combined hemimethylation and internal cytosine methylation) showed low (Fig. 2; N. robustus = 0.272; N. bivocatus = 0.240) but significant variation between species (Mann-Whitney U, W = 587, p = 0.001243).

Figure 2. 

Comparison of genome-wide methylation levels between species. Mann-Whitney U test; p<0.005 (**). Between-species significant variation is in total methylation (internal cytosine methylated and hemimethylated fragments).

Epigenetic and genetic structure.—We examined epigenetic and genetic diversity as they relate to species assignment. The within-species epigenetic Shannon diversity index was 5.0616 ± 0.2054 and 5.0094 ± 0.2087, within N. robustus and N. bivocatus, respectively. N. robustus and N. bivocatus’ genetic Shannon diversity indexes were 5.2802 ± 0.2242 and 5.2179 ± 0.1943, respectively. Diversity in the genetic loci, as measured in this study, was significantly greater than that of the MSL (epigenetic) diversity measured in both species (Wilcoxon rank sum test; N. robustus, W = 2760, p < 0.0001; N. bivocatus, W = 779, p = 0.0001).

Epigenetic and genetic structure was evaluated using a two-level AMOVA. The between-species epigenetic divergence was Φ_ST = 0.0504 (P < 0.0001). This is slightly less than one-third of the genetic divergence observed between species, Φ_ST = 0.1591 (P < 0.0001). Greater genetic variation between species would suggest that genetic mechanisms are underlying species differentiation (Table 2).

The principal coordinate analyses (PCoA), calculated using MSL and genetic profiles, showed that genetic variance was smaller than epigenetic variance within species (Fig. 3). The epigenetic PCoA between species elucidated little meaningful differentiation among groups, with ellipses (95% confidence intervals) overlapping almost entirely. The PCoA of genetic data between species explained 17.7% of the genetic variance in the first two coordinates (Fig. 3A) and showed moderate variation between individuals based on call phenotype, with species showing divergence in Euclidean space. Three N. robustus individuals clearly fell into the N. bivocatus cluster within the genetic PCoA (Fig. 3A), indicating a possible mismatch between phenotypic and genotypic assignments.

The Bayesian analysis of genetic structure revealed the best-supported number of genetic clusters to be K = 2, based on ΔK values (Fig. 4B; Evanno et al. 2005). These two genetic clusters aligned closely to species boundaries (Fig. 4A). As in the PCoA, the same three individuals identified as N. robustus were assigned primarily to the N. bivocatus genetic cluster, indicating that they possess the N. robustus call type but a N. bivocatus genotype. As in the PCoA, the epigenetic structure analysis did not show support for significant epigenetic (MSL) structure between species groups. ΔK values did not rise above 30 for any permutation in the number of clusters (Fig. 4D), an indication of low support for epigenetic structure. In addition, the best supported number of clusters, K = 4, showed no grouping of epigenetic clusters by species assignment (Fig. 4C).

While significant epigenetic variation was observed between species, significant genetic variation was also detected. We investigated the correlation between inter-individual genetic and epigenetic Euclidean distance. Between-individual epigenetic and genetic distance showed a strong positive correlation (Fig. 5, Mantel R = 0.8448, p = 0.0001).

Figure 3. 

PCoA of N. robustus and N. bivocatus utilizing genetic (A) and epigenetic (B) data. Plotted are the two most informative principal components calculated for the genetic and epigenetic loci datasets, as derived from the MS-AFLP fragment analysis. A. Genetic Euclidean distance with individuals grouped by species assignment. B. Epigenetic Euclidean distance with individuals grouped by species assignment. Group labels show the centroid of the points for each group. The long axis of the ellipse represents the direction of maximum dispersion and the short axis the direction of minimum dispersion.

Figure 4. 

Consensus shared ancestry population structure for epigenetic and genetic loci. A. and C. Bar plots using MS-AFLP loci to estimate genetic (A) and epigenetic (C) structure among N. robustus and N. bivocatus using the software package STRUCTURE. B. and D. Delta K graphs for K = 1–10 genetic clusters showing moderate support for K = 2 genetic clusters (B) and low support for K = 4 epigenetic clusters (D).

We found significant differentiation among species in both epigenetic and genetic markers. Genetic differentiation, however, was larger than epigenetic differentiation between species. The AMOVA and PCoA analyses indicated that methylation patterns varied among individuals but showed little differentiation between species. Epigenetic distance among individuals correlated with inter-individual genetic differentiation, suggesting that the low levels of epigenetic differentiation between species have been pulled along by genetic differentiation.

Genetic differentiation.—Genetic differentiation between species was low but significant, as would be expected between two closely related taxa (Snyder et al. 2009, Ney and Schul 2017). Genetic structure aligned with species boundaries, with the exception of the three phenotypically N. robustus individuals that showed a majority assignment to the primarily N. bivocatus genotypic cluster. These mismatch genotype/phenotype assignments were confirmed for these three individuals in both the genetic structure analysis (Fig. 4A) as well as the genetic PCoA (Fig. 3A). Similar occurrences of mismatched genotype/phenotype individuals have been found between N. robustus and N. bivocatus previously (Frederick 2013, Ney and Schul 2017). These mismatched genotype/phenotype individuals are likely the result of recent hybridization events. A severe drought during the sampling period (Ney and Schul 2017) led locally to sharp population declines of one or both species. Both reductions in population sizes and environmental disturbances can increase rates of hybridization between closely related species (Lamont et al. 2003, Seehausen 2006).

Epigenetic variation.—N. robustus showed a significantly higher level of genome wide methylation than N. bivocatus (Fig. 2). This may be due to variation in the species’ habitat preference. Variable environmental conditions can cause shifts in epigenetics, observed between diverging natural habitats as well as between natural and altered habitat types (Gao et al. 2010). N. robustus, unlike N. bivocatus, is a grassland generalist and may show epigenetic patterns resulting from exposure to a more variable set of environmental conditions during a lifetime and across generations. Because of this, N. robustus may possess a larger repertoire of adaptive genes for the varied environmental conditions it may utilize. The higher level of methylation found in N. robustus may therefore be the result of higher levels of methylation required to silence the large repertoire of adaptive genes when not in use.

While epigenetic differentiation between species was significant, it remained lower than genetic differentiation. This study did not show support for the epigenetic regulation of species-specific phenotypes; however, this does not eliminate a possible epigenetic mechanism underlying phenotypic differentiation. While the MS-AFLP technique has many benefits, there are also some inherent limitations to its application. For example, MS-AFLPs can underestimate genome-wide levels of methylation (Fulneček and Kovařík 2014). As a genome-wide scanning method MS-AFLP only detects methylation at 5′-CCGG-3′ sites. In addition, it is unable to discriminate between full methylation states (hypermethylation of both cytosines) and sequence variation at or near the restriction site (Xiong et al. 1999). Phenotypic variation controlled by a smaller number of methylated sites is unlikely to be detected. A similar investigation into genome-wide methylation found no significant differentiation in methylation between the solitarious and gregarious phases of Locusta migratoria (Robinson et al. 2011), despite divergent phenotypes and the identification of differentially expressed genes between morphs (Kang et al. 2004). Thus, this technique is most valuable as a first step in investigating genome-wide methylation patterns.

Correlation between epigenetic and genetic diversity.—MS-AFLP variation often shows correlations with genetic variation (Liu et al. 2012). In our study, individuals showing greater genetic variation on average showed greater epigenetic variation as well (Fig. 5). Changes in DNA methylation over time are correlated with genetic relatedness and suggest that DNA methylation maintenance may be under genetic control (Bjornsson et al. 2008). In addition, genetic variation in retrotransposons (i.e. their presence/absence) could affect the methylation state of retroelements, requiring more methylation if more repetitive elements are present (Michaud et al. 1994). Repetitive elements have been estimated to make up thirty percent of the genome of the orthopteran L. migratoria (Wilmore and Brown 1975). In addition, the retrotransposon SINE is differentially methylated between the solitarious and gregarious phases of the species (Guo et al. 2010).

Mechanisms of neutral evolution could also account for the correlation between epigenetic and genetic variation. In the event of substantial gene flow between species, strong divergent selection would be needed to maintain divergent epigenetic variation between species. Gene flow between groups would reduce differentiation accumulated via drift. Evidence from this study, however, suggests that genetic differentiation between species is significant and gene flow therefore relatively low (Table 2). Stochastic processes of drift then could allow epigenetic patterns to diverge between species in parallel, resulting in correlations between epigenetic and genetic variation without a causal link (Richards et al. 2010).

Figure 5. 

Scatterplot of between-individual Euclidean genetic and epigenetic distance showing significant positive correlation between genetic and epigenetic differentiation. The correlation was tested using a Mantel test and 10,000 permutations of the design matrix to determine significance.

We hypothesized that the phenotypic variation observed between N. robustus and N. bivocatus was the result of an epigenetic-mediated mechanism of species differentiation that led to genetic isolation. While we found clear evidence of genomic methylation in Neoconocephalus, our findings did not support a role for methylation in species isolation. Both the lower level of epigenetic differentiations and the correlation between inter-individual epigenetic and genetic diversity support the alternative hypothesis, i.e. that differences in methylation patterns between species evolved in response to genetic variation. Epigenetics may still play a key role in phenotypic differentiation within Neoconocephalus katydids through the differential regulation of a relatively small number of genes of large effect; a mechanism not detectable with the methods used here. Further work identifying differentially expressed genes between call types could allow for the targeted analysis of methylation patterns at these sites.

We thank Katy Frederick and Nathan Harness for assistance with specimen collection and recording males. We thank Kim Hunter, Katy Frederick, and Megan Murphy for valuable feedback on the manuscript.