Evolution of a key trait greatly affects underground community assembly process through habitat adaptation in earthworms

Abstract Underground community assemblies have not been studied well compared with aboveground communities, despite their importance for our understanding of whole ecosystems. To investigate underground community assembly over evolutionary timescales, we examined terrestrial earthworm communities (Oligochaeta: Haplotaxida) in conserved mountainous primary forests in Japan as a model system. We collected 553 earthworms mostly from two dominant families, the Megascolecidae and the Lumbricidae, from 12 sites. We constructed a molecular taxonomic unit tree based on the analysis of three genes to examine the effects of a biogeographic factor (dispersal ability) and an evolutionary factor (habitat adaptation) on the earthworm community assembly process. The phylogenetic distance of the earthworm communities among sites was positively correlated with geographic distance when intraspecific variation was included, indicating that the divergence within species was affected by biogeographic factors. The community assembly process in the Megascolecidae has also been affected by environmental conditions in relation to an evolutionary relationship between habitat environment and intestinal cecum type, a trait closely related to habitat depth and diet, whereas that in the Lumbricidae has not been affected as such. Intestinal cecum type showed a pattern of niche conservatism in the Megascolecidae lineage. Our results suggest that investigating the evolution of a key trait related to life history can lead to the clear description of community assembly process over a long timescale and that the community assembly process can differ greatly among related lineages even though they live sympatrically.

timescale. Two main factors explain why a species lives in a particular place (Webb, Ackerly, McPeek, & Donoghue, 2002): biogeographic factors control whether the species is able to arrive at that place via dispersal, and evolutionary factors control whether the ancestor of the species had evolved the traits necessary to live in that particular environment.
Related species tend to share similar ecological traits that are adaptive in similar environments (a phenomenon known as ecological niche conservatism), although related species may also evolve different habitat preferences as a consequence of adaptive radiation (Emerson & Gillespie, 2008). Thus, if habitat adaptation is conserved through the phylogeny, and no competitive exclusion occurs, community assemblages at sites with similar environments tend to be composed of phylogenetically closely related species (Vamosi & Vamosi, 2007;Webb et al., 2002).
Studies of the community assembly process over evolutionary timescales have been conducted on aboveground plants and insects and in aquatic ecosystems, but few such studies have investigated organisms living underground. The underground environment is very different from those aboveground, in showing less annual and seasonal variability. In addition, soil animals have a poor dispersal ability and consequently show high genetic divergence among areas (Chang, Lin, & Chen, 2008;Minamiya, Yokoyama, & Fukuda, 2009). Therefore, the process of forming animal community assemblages underground is expected to differ greatly from that aboveground. Although several community phylogenetic studies have investigated soil animals such as ants (Smith, Hallwachs, & Janzen, 2014) and earthworms , these studies did not consider the evolutionary ecology of habitat adaptation, which can greatly affect community assembly.
In this study, we address the effects of biogeographic and evolutionary factors on the community assembly process of underground animals by using terrestrial earthworms (Oligochaeta: Haplotaxida) as a model system. Earthworms are appropriate representatives of soil animals, as they make large contributions to the decomposition of litter and soil organic matter (Darwin, 1881;Edwards & Lofty, 1977) and behave as ecosystem engineers that create and destroy habitat for other organisms by modifying the habitat where they live (Lavelle et al., 2016). We examined their dispersal ability as a biogeographic factor and their habitat adaptation as an evolutionary factor.
Numerous invasive earthworms are known worldwide (Blakemore, Ito, Kaneko, et al. 2006;Hendrix et al., 2008), so we set our 12 study sites in conserved mixed coniferous-broadleaved or beech forests at relatively high elevations (650-1,920 m), which we expected to contain undisturbed soil animal communities without invasive earthworms.
The ecological niches of megascolecid earthworms are related to their intestinal cecum type; Japanese native megascolecid earthworms have a pair of intestinal ceca connected with the intestinal tract (Ishizuka, 2001;Ishizuka & Minagoshi, 2014). The shape of the intestinal cecum is related to the digestive enzyme that it produces (Nozaki, Ito, Miura, & Miura, 2013), which in turn is related to ecological factors such as habitat depth and diet (Ishizuka, 2001;Ishizuka & Minagoshi, 2014;Uchida et al., 2004). Those species with the simple (Figure 1a) and serrate types of ceca live in the soil layer and feed on decomposed organic matter; those with serrate ceca live in the deep soil layer and are seldom collected from the topsoil layer. Species with the manicate type ( Figure 1b) live in the litter layer and feed on less decomposed organic matter. Owing to the evolution of different cecum types among members of the Megascolecidae, we expect that the community assembly process of megascolecid earthworms in Asia is more strongly affected by evolutionary factors than by biogeographic factors. We hypothesize that the evolution of intestinal cecum type in the Megascolecidae has led to an evolutionary change in habitat environment and thus affects community assembly of this family.

| Study sites and methods
We conducted our study at 12 sites (20 m × 20 m) in primary mixed coniferous-broadleaved forests or beech forests in mountainous areas in the central part of mainland Japan (Figure 2, Appendix S1).
We collected earthworms in early summer and in early autumn (sites

| Morphology of earthworms
We identified earthworm specimens to the family level from their morphology. Because most earthworm species in Japanese primary forests have not been described, no morphological identification keys of earthworms at the species level have been established (Ikeda et al., 2012;Ishizuka, 2001;Minamiya et al., 2009). In addition, most collected samples were juveniles, which can be identified only to the family level from morphology. Thus, we used molecular operational taxonomic units (MOTUs) in this study (see "Phylogenetic analysis" section for details). We also examined the intestinal cecum of all intact megascolecid earthworms.
We measured the wet weight of collected litter samples and the dry weight after oven-drying them at 70°C for more than 3 days. Large branches, which are not consumed by earthworms, were removed.
Dried litter samples were ground, and carbon and nitrogen contents were measured with a Sumigraph NC-22F CN analyzer (Sumika Chemical Analysis Service, Ltd., Tokyo, Japan). The measured data were calibrated against an acetanilide standard.
Soil samples were oven-dried at 70°C for more than 3 days and then passed through a sieve with 2-mm-diameter openings to remove gravel, stones, and roots. The gravel and stones were oven-dried at 105°C overnight, and roots were oven-dried at 105°C for a day to measure their dry weight. To calculate water content, we oven-dried 3 g of each soil sample at 105°C overnight. We calculated the soil bulk density per 400 ml from these data. We put 10 g of oven-dried soil into 25 g of ion-exchanged water to measure soil pH on a portable pH meter (Horiba, Kyoto, Japan). Dried soil was ground, and carbon and nitrogen contents were measured with a Sumigraph NC-22F CN analyzer.
We conducted principal component analysis for the average value of each environmental parameter at each site in R v. 3.2.4 software with standardization by the "prcomp" command (R Core Team, 2016).
We considered the difference in PC scores between sites to represent the difference in habitat environment for earthworms between sites.
3.1 Cycle Sequencing Kit (Applied Biosystems). The products were electrophoresed in an ABI 3130XL sequencer (Applied Biosystems).
We obtained at least one gene region from all specimens. Sequence data are deposited in GenBank (GenBank accession numbers: LC199973-LC200394). We used MAFFT v. 7.273 software (qinsi method) for the 16S alignment (Katoh & Standley, 2013). The alignments were inspected by eye for obvious misalignments. Haplotypes were detected using the "pgelimdupseq" command in Phylogears2 v. 2.0 software (Tanabe, 2008), and individuals with the same haplotypes were excluded from phylogenetic analyses.

| Phylogenetic analysis
We LC228479-LC228482). First, we constructed a haplotype tree with all haplotypes to include the divergence within species. We used the log-normal relaxed clock model for the analysis. The estimated clock rate for earthworms (Chang et al., 2008) was used as the prior for the clock rate of COI (initial value = 0.024). We performed three Markov chain Monte Carlo runs for 300 million generations, with trees sampled every 10,000 generations; the first 25,000 trees were discarded as burn-in. We combined the remaining trees from three runs and used them for constructing a tree.
The node ages at the divergence points of the Megascolecidae + Lumbricidae and the Moniligastridae were used for setting the node We split potential species (MOTUs) from the phylogenetic trees that include the divergence within species (haplotype tree), because most earthworm species have not been described in Japanese primary forests, and thus morphological identification keys have not been established (Ikeda et al., 2012;Ishizuka, 2001;Minamiya et al., 2009). We used a multilocus species delimitation method implemented in the program "tr2" (Fujisawa, Aswad, & Barraclough, 2016). The haplotype tree was used as a guide tree for analysis.
We also constructed each gene tree for analysis in BEAST v. 1.8.4 software (Drummond et al., 2012). We used the same substitution models and settings for calibration that were used for the haplotype tree. We performed a Markov chain Monte Carlo run for 200 million generations, with sampling every 10,000 generations for each gene.
The first 10,000 trees were discarded as burn-in. We used the remaining trees for constructing each gene tree. The separated potential species were considered as species in the analysis to construct the MOTU tree.
We constructed the MOTU tree, in which the divergence among haplotypes within species was considered as the variation within species, in the program *BEAST (Heled & Drummond, 2010) in BEAST v. 1.8.4 software (Drummond et al., 2012). We used the same substitution models and settings for calibration that were used for the haplotype tree except for the random local clock model. We performed three Markov chain Monte Carlo runs for 400 million generations, with trees sampled every 10,000 generations, and the first 35,000 trees were discarded as burn-in. We combined the remaining trees from the other three runs and used them for constructing trees. This MOTU tree was used for calculating the weighted Unifrac distance described below. We calculated rarefaction curves of MOTUs in EstimateS v. 9 software (Colwell, 2013) to compare species richness among families.

| Testing phylogenetic distance of community assemblages among sites
To examine the phylogenetic distance of community assemblages among sites, we calculated weighted Unifrac distance using the MOTU tree by using the "unifrac" command in the phyloseq package (2) node M2 (2) node M1 (1) Megascolecidae node (1)  distance is expected. We tested the positive correlation between weighted Unifrac distance and geographic distance or environmental differences (PC1 and PC2 scores). To remove the correlation between geographic distance and environmental differences (geographic distance and PC1 score: p = .042; geographic distance and PC2 score: p = .152 by Spearman's rank correlation test), we conducted a onetailed partial Mantel test with Spearman's rank correlation coefficient by using the "mantel" command with 1 million permutations in the "ecodist" package in R. Because the strength of the effect of each factor can vary over evolutionary timescales, we investigated the correlation at the following time points: (1) Megascolecidae node and Lumbricidae node; and (2) nodes within the Megascolecidae (see Figure 3). We also conducted the same analysis using the haplotype tree including the nodes within the Lumbricidae (see Figure 4).

| Ancestral state reconstruction
To examine the evolutionary relationship between habitat environment and intestinal cecum type and the niche conservatism in these traits in the Megascolecidae, we reconstructed the ancestral states for habitat environments of megascolecid earthworms (PC1 and PC2 scores) by using parsimony ancestral state reconstruction with the squared change assumption in Mesquite v. 3.11 (Maddison & Maddison, 2016). PC1 and PC2 scores at the site where each earthworm sample was collected were set as the habitat environments for that sample. PC1 and PC2 scores of the species collected at more than one site were weighted by the number of samples collected at each site. Ancestral state reconstruction was conducted for the megascolecid MOTU tree.
The ancestral states for intestinal cecum types of megascolecid earthworms were also reconstructed by parsimony ancestral state reconstruction in Mesquite v. 3.11 (Maddison & Maddison, 2016).
To examine the difference in habitat environments between the simple cecum type lineage and the manicate cecum type lineage, we used the phylogenetic generalized least squares method to control for phylogenetic signal by using the "pgls" command in the caper package (Orme et al., 2013) in R. We tested the correlation between the evolutionary change of intestinal cecum type and that of habitat environment using the MOTU tree. We compared AIC scores between the model with the evolutionary change in cecum type and the null model.

| RESULTS
We collected 553 earthworms (230 Megascolecidae, 282 Lumbricidae, 41 Moniligastridae) in total from the litter and topsoil layers (Appendix S2) at the 12 sites (Figure 2, Appendix S1). Significantly more megascolecid earthworms with manicate ceca than with simple ceca were collected from the litter layer (Fisher's exact test, p < .001; Figure 3).  (2) Relationships of weighted Unifrac distance with (d) geographic distance and (e, f) environmental differences (PC1 and PC2, respectively) at two nodes within the Megascolecidae (M1, n = 15; M2, n = 21) and two nodes within the Lumbricidae (L1, n = 45; L2, n = 45). Statistical results of Spearman's rank correlation tests are also shown (see Table 3 for details). Significant P values are in bold italic type. Solid lines show the regression lines for significant relationships   was highly positively correlated with soil water content, proportion of carbon in soil, and depth of the A 1 layer, and it was highly negatively correlated with soil bulk density (Table 1, Appendix S1). The contribution of PC2 was 25.0%. PC2 was highly positively correlated with C:N ratio in litter and soil, and it was highly negatively correlated with the proportion of carbon in litter and litter depth (  Tables 2 and 3).
Our phylogenetic tree indicated that the manicate cecum type evolved from the simple cecum type once, and each cecum type was clustered in the tree (Figure 9a). The PC1 score was higher and the PC2 score was lower in the lineage after the evolutionary change from simple to manicate cecum type than in the lineage with simple cecum type, according to the phylogenetic generalized least squares method should be suitable habitat for underground decomposers owing to sufficient habitat space with moist, soft, carbon-rich soil, which would provide abundant food resources for earthworms (see Appendix S1).
The sites with a lower PC2 score should be suitable habitat for litter decomposers owing to a deep and nutritionally rich litter layer, as a litter layer with a low C:N ratio has a high nitrogen concentration.

| DISCUSSION
The Unifrac distance was positively correlated with geographic distance in the haplotype tree, and none of the haplotypes were shared among sites. These results indicate a very low rate of migration among primary forests due to the earthworms' poor dispersal ability. However, this pattern was not detected in the MOTU tree or at the Megascolecidae node in the haplotype tree, suggesting that the low rate of migration in earthworms causes geographic isolation at a relatively short timescale. Such a trend is likely common in soil-inhabiting animals, because their dispersal ability is considered to be poorer than that of aboveground animals (Chang et al., 2008;Minamiya et al., 2009 Simple evolutionary change in a key trait has greatly affected the process of community assembly over an evolutionary timescale. In contrast, our findings indicate that evolutionary factors have not affected the process of community assembly in the lumbricid lineage over an evolutionary timescale, whereas biogeographic factors have affected it (Figures 7 and 8; Tables 2 and 3). Although some lumbricid species live in the soil layer and others in the litter layer (Bouché, 1977;Domínguez et al., 2015;Perel, 1977), we detected no evolutionary pattern of habitat preference in our study. This pattern would be partly caused by low species diversity in the Lumbricidae in Japanese primary forests. In addition, the fundamental niche width in lumbricid species may be wider than that in megascolecid species.
Using earthworms as a model system, we demonstrated that the underground animal community of each area was unique and had a long divergence history affected by biogeographic factors owing to the animals' poor dispersal ability. Although previous studies for soil animals have not investigated the evolutionary ecology of habitat adaptation, our findings also show that the process of community assembly in megascolecid earthworms has been driven by the evolution of intestinal cecum morphologies through the change of habitat preference. We found a clear difference in the community assembly process between the Megascolecidae and the Lumbricidae, even though species from both lineages live sympatrically. Our results suggest that investigating the evolution of a key trait, such as the intestinal cecum, that is closely related to life history can lead to the clear description of the community assembly process over an evolutionary timescale.

ACKNOWLEDGMENTS
We thank Y. Igarashi and K. Ohmura