Comparative characterization of microbiota between the sibling species of tea geometrid moth Ectropis obliqua Prout and E. grisescens Warren

Abstract For a wide range of insect species, the microbiota has potential roles in determining host developmental programme, immunity and reproductive biology. The tea geometrid moths Ectropis obliqua and E. grisescens are two closely related species that mainly feed on tea leaves. Although they can mate, infertile hybrids are produced. Therefore, these species provide a pair of model species for studying the molecular mechanisms of microbiotal involvement in host reproductive biology. In this study, we first identified and compared the compositions of microbiota between these sibling species, revealing higher microbiotal diversity for E. grisescens. The microbiota of E. obliqua mainly comprised the phyla Firmicutes, Proteobacteria and Cyanobacteria, whereas that of E. grisescens was dominated by Proteobacteria, Actinobacteria and Firmicutes. At the genus level, the dominant microbiota of E. grisescens included Wolbachia, Enterobacter and Pseudomonas and that of E. obliqua included Melissococcus, Staphylococcus and Enterobacter. Furthermore, we verified the rate of Wolbachia to infect 80 samples from eight different geographical populations, and the results supported that only E. grisescens harboured Wolbachia. Taken together, our findings indicate significantly different microbiotal compositions for E. obliqua and E. grisescens, with Wolbachia possibly being a curial factor influencing the reproductive isolation of these species. This study provides new insight into the mechanisms by which endosymbiotic bacteria, particularly Wolbachia, interact with sibling species.


Introduction
Important interactions of microbiota organisms with insects are very common in nature (Mao et al., 2018). Some of the microbiota harboured in host cells are considered endosymbionts, constituting a symbiotic bacteriome. Other microbiota species opportunistically colonize different tissues of insects and the gut lumen, which can be affected by many factors, such as the host's diet and living environment (Colman et al., 2012;Engel and Moran, 2013;Huang and Zhang, 2013). Overall, the microbiota plays an important role in an insect's life activity, providing the host with essential nutrients and protection from predators, parasites and pathogens (Tsuchida et al., 2010) and affecting reproduction (Zhang et al., 2015). Indeed, an increasing number of studies are focusing on the use of sequencing technology to analyse the functions and development of microbial communities associated with termites, silkworms and other insects (Su et al., 2016;Sun et al., 2016;Chen et al., 2017;Wang et al., 2017). However, the effect of the microbiota on host insects remains incompletely understood.
Ectropis obliqua and E. grisescens are primary defoliators in tea plantations due to their wide distribution and destructive nature (Jiang et al., 2014;Zhang et al., 2014). These moths infest thousands of hectares of tea per year, severely reducing the growth and impacting tea production in the following year (Jiang et al., 2014;Zhang et al., 2016a). Ectropis obliqua and E. grisescens were named in 1894 and 1930, respectively, but they have always been treated as the same species in tea garden management in China because they have similar morphology and are difficult to be distinguished (Jiang et al., 2014). After the two species were reacquainted in 2014, differences in morphology, reproductive capacity and virus susceptibility have been reported (Jiang et al., 2014;Mao et al., 2017;Bai et al., 2018a). In particular, the technique of DNA barcoding can accurately distinguish them according to the genetic distance of approximately 3.7% based on COI sequences (Jiang et al., 2014). Correspondingly, the distribution of the two species has become clearer, with E. obliqua being distributed in China, Japan and the Korean Peninsula and E. grisescens only in China (Wehrli, 1945;Sato, 1984;Kim et al., 2001). In Zhejiang Province, in eastern China, E. grisescens is more widespread than E. obliqua, though they are both also found in some areas (Bai et al., 2018b). Moreover, morphological and phylogenetic evidence supports that E. obliqua is closely related to E. grisescens. Intriguingly, these sibling species can mate but produce infertile hybrids (Xi et al., 2014). Indeed, the hybrid F1 generation showed hatching, survival to adult stage and per cent of normal adult rates that were much lower than those with intra-species mating. Furthermore, a self-cross of F1 generation adults produced either infertile eggs or no eggs (Xi et al., 2014;Zhang et al., 2014). Reproductive interference also exists between these sibling species (Zhang et al., 2016a). As these phenomena differ from those resulting from common reproductive isolation, whereby different species cannot mate and breed, we suggest that these sibling species constitute a suitable model pair for exploring reproductive isolation.
Previous studies have shown some microbiota organisms can manipulate host reproduction and even cause reproductive isolation (Philipp and Nancy, 2013;Zhang et al., 2016aZhang et al., , 2016b. Moreover, we previously found that F1 hybrids of E. obliqua and E. grisescens showed the characteristics including unbalanced sex ratio, lower hatchability of eggs, desynchronized development of larvae and infertility that resembled cytoplasmic incompatibility which was caused by some microbiota (Bourtzis et al., 1996;Zhang et al., 2014;Wang et al., 2019). Thus, we sought to ascertain the microbiota involved in the reproductive isolation of these sibling species. In this study, we analysed and compared differences in E. grisescens and E. obliqua microbiotal composition and found that an obvious difference in the presence of Wolbachia may be a factor influencing their reproductive isolation.

Collection and preparation of samples
Ectropis obliqua larvae were collected from Yuhang (Hangzhou City, Zhejiang Province) and E. grisescens from Xinchang (Shaoxin City, Zhejiang Province). At least 200 larvae were collected in each location. The collected larvae were reared in a phytotron (temperature 24-26°C, humidity 50-70%, photoperiod L14:D10). The larvae were fed fresh leaves of the tea cultivar Yingshuang for successive three generations. Male and female individuals were separated at the pupa stage. Two days after eclosion, adult moths were randomly collected for the analysis of bacterial communities. In addition, samples from different geographical populations were collected to evaluate the rate of Wolbachia infection without rearing in the phytotron. Sampling localities of tea geometrid can be seen in fig. 1.

DNA extraction
Wings were removed, and other tissues were washed with sterilized water and ground for 2 min. The tissue homogenate was used for metagenomic DNA extraction using the DNeasy Blood and Tissue kit (Qiagen Co. Inc., Germany) according to the manufacturer's instruction. The quality of the extracted DNA was assessed by electrophoresis on a 1% (w/v) agarose gel. The concentration of DNA extracted was measured using a Nanodrop 2000, and the DNA samples were stored at −20°C for experiments.

Sample identification
Identification of E. obliqua and E. grisescens was confirmed by sequence analysis of the mitochondrial cytochrome oxidase I gene (COI) (Jiang et al., 2014). A fragment of the COI gene was amplified using the forward primer LepF1 5 ′ -ATTCAACCAAT-CATAAAGATATTGG-3 ′ and the reverse primer Enh_LepR1 5 ′ -CTCCWCCAGCAGGATCAAAA-3 ′ (Jiang et al., 2014). PCR using 2x Master Mix (TSINGKE Bio Inc., Hang zhou City, Zhe jiang, China) in a total reaction volume of 50 μl was performed using a Veriti 96 Well Thermal Cycler (Applied Biosystems, Foster City, CA, USA), as follows: an initial denaturation step of 95°C for 2 min, 5 cycles of 95°C for 30 s, 46°C for 1 min, and 72°C for 30 s, 35 cycles of 95°C for 30 s, 51°C for 1 min, and 72°C for 30 s and a final extension at 72°C for 10 min. Sequencing of PCR products was performed using an ABI377 genetic analyser. Sequences were assembled and edited with SeqMan 7.1.0. and aligned with CLUSTAL 1.83. The sequence of E. obliqua COI gene (accession number KJ704358), which was obtained from GenBank, was used as a control to calculate genetic distance based on Kimura-2-parameter model (100 bootstrap replicates) and build a cluster analysis tree using a maximum likelihood approach implemented in MEGA 5.05. Compared to the control sample (KJ704358), genetic distances of ∼0-2.5 and 3.2-4.0%, respectively, were identified for E. obliqua and E. grisescens (Jiang et al., 2014).

PCR amplification of microbial 16S rDNA genes
The V3-V4 region of 16S rDNA was amplified using the forward primer 341F 5 ′ -CCTACGGGNGGCWGCAG-3 ′ and the reverse primer 805R 5 ′ -GACTACHVGGGTATCTAATCC-3 ′ (Sinclair et al., 2015). PCR reactions were carried out with Phusion® High-Fidelity PCR Master Mix (Thermo Scientific, Waltham City, MA, USA), and the conditions of PCR amplification are as follows: pre-amplification was 94°C for 3 min, 5 cycles of 94°C for 30 s, 45°C for 20 s, and 65°C for 30 s, followed by 20 cycles of 94°C for 20 s, 55°C for 20 s, and 72°C for 30 s and a final extension at 72°C for 5 min. PCR products were examined on a 2% agarose gel. Only samples with a clear band between 400 and 450 bp were chosen for further experiments.

Sequencing of 16S rDNA gene amplicons
Sequence libraries that included 20 individuals were generated using TruSeq® DNA PCR-free Sample Preparation Kit (Illumina) according to the manufacturer's instruction. Library quality was assessed using a Qubit@ 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system (Agilent). The libraries were sequenced using the Illumina Miseq platform by Zhejiang Tianke High Technology Development Co. Ltd. (Zhejiang, China), and 300 bp paired-end (PE) reads were generated.

Detection of the Wolbachia infection rate
To determine whether the samples were infected with Wolbachia, the wsp gene was amplified from genomic DNA (Gong and Shen, 2002;Laura et al., 2006). The forward primer wsp_F1 5 ′ -GTCC-AATARSTGATGARGAAAC-3 ′ and the reverse primer wsp_R1 5 ′ -CYGCACCAAYAGYRCTRTAAA-3 ′ were used (Laura et al., 2006), and amplification was as follows: denaturation at 95°C for 3 min, 35 cycles of 95°C for 1 min, 59°C for 1 min, and 72°C for 90 s and a final extension at 72°C for 10 min. PCR products were detected by 1% agarose gel electrophoresis.

Data analysis
PE reads were merged using FLASH, and quality filtering of spliced sequences (raw tags) was performed under specific conditions to obtain high-quality clean tags according to QIIME. To detect and remove chimeric sequences, the tags were compared with reference sequences obtained from the Gold database using the UCHIME algorithm. Sequence analysis was performed using Uparse software. Sequences with ≥97% similarity were assigned to the same operational taxonomic unit (OTU). Representative sequence for each OTU was screened for further annotation.
Taxonomic assignment was achieved using the SILVA reference database (http://www.arb-silva.de/) with a threshold of 90%. The α diversity was applied when evaluating the complexity of species diversity for a sample using six indices: observed species, Chao1, Shannon, Simpson, ACE, Good's coverage. All indices were calculated using QIIME (Version 1.9.1). The β diversity analysis was applied to evaluate the differences in species complexity among the different samples, and principal coordinate analysis (PCoA) was performed based on the matrices of pairwise distances among all the microbiota. Heatmap and hierarchical cluster were built based on the relative abundance of the top 15 genera identified in the bacterial communities of samples by using pheatmap package in R program. The highly relative abundances of microbiota phyla and genera were visualized using ggplot2 package (version 3.0.0) in R program. Statistical analyses were performed using SPSS 17.0 (IBM). One-way ANOVA was followed by Tukey test for means comparison of α diversity. The level of significance was set at P < 0.05.

Species verification
A total of 20 individuals were chosen for bacterial community analysis. All samples from the Yuhang population were identified as E. obliqua and numbered O1-10 (group O); all samples from the Xinchang population were identified as E. grisescens and numbered G1-10 (group G) ( fig. 2; table S1).

Sequencing data and sequence read diversity analysis
A total of 752,452 raw PE 16S rDNA reads were generated from these 20 samples. After removal of low-quality reads, 658,249 (∼75.5%) valid reads were obtained. The length of each read was between 404 and 426 bp, with an average of 417.2 bp. Q20 values for all samples were >98%. All reads clustered into 1492 OTUs (3% distance, average neighbour clustering), covering ten phyla, 40 classes, 79 orders, 118 families, 198 genera and 125 species.
Shannon and Simpson indices were used to evaluate bacterial diversity, and Chao1 and Ace indices were employed to estimate the total number of species in the samples (Sun et al., 2016); Good's coverage was applied for sequencing results.
Significantly higher values of the Simpson index were found for E. grisescens compared to E. obliqua, whereas E. obliqua displayed
significantly higher Ace and Chao1 index values (table 1). Despite a lack of a significant difference for the Shannon index between E. grisescens and E. obliqua, higher bacterial community diversity was found for E. grisescens. Overall, the high Good's coverage (both >99%) values suggest that the OTUs covered most of the bacterial communities present and that our metagenomic data were reliable.
To compare similarity and dissimilarity between all samples, PCoA and hierarchical clustering analysis were performed. The points in fig. 3 represent samples in the PCoA plot, and the distance of each point indicates the similarity of different samples. PCoA separated the samples into two clusters, with 75.19 and 17.80% of the total variation being explained by the PCo1 and PCo2 axes, respectively. With the exception of sample O2, which was closer to group G, samples of the same species clustered together. In addition, hierarchical clustering analysis separated the samples into two clades according to species (fig. 4). Overall, the results suggest that the bacterial communities in these two species differed and that the major contribution of the difference was PCo1.
The 16S rDNA V3-V4 region was amplified to compare the differences of the microbiota of these sibling species of tea geometrid moths. A total of 286,700 valid reads and 1043 OTUs were obtained from ten of E. obliqua samples, comprising ten phyla, 38 classes, 75 orders, 111 families, 179 genera and 118 species. Among these ten phyla ( fig. 5), Firmicutes (52.80%), Proteobacteria (24.52%) and Cyanobacteria (14.05%) were highly abundant. At the genus level, 14.48% reads were not identifiable, and the remaining reads belong to the bacteria of 179 genera. The predominant genera (>1%) were Melissococcus (29.14%),  In total, 371,549 valid reads and 449 OTUs were obtained from the sequencing data for E. grisescens, covering seven phyla, 19 classes, 34 orders, 58 families, 75 genera and 40 species. Compared to E. obliqua, fewer microbiota organisms were identified in E. grisescens at all classification levels, and the composition of abundant taxa varied in these two species. The bacteria found in E. grisescens belong to seven phyla (fig. 5), with Proteobacteria (79.13%), Actinobacteria (14.06%) and Firmicutes (5.13%) being the most abundant. Comparatively, the abundances of Firmicutes and Cyanobacteria were significantly greater in E. obliqua than in E. grisescens, whereas the abundances of Proteobacteria and Actinobacteria were higher in E. grisescens ( fig. 5).

Detection of Wolbachia in sibling species of tea geometrid moths
We screened for four microbiota (Wolbachia, Cardinium, Spiroplasma and Rickettsia), which have been shown to influence the reproduction (Zhang et al., 2016b). The results showed that only Wolbachia was found in both species. Thus, we evaluated the richness of Wolbachia in all samples used for bacterial community analysis and found it to be significantly different between the sibling species of tea geometrid moths. Specifically, the average richness of Wolbachia was 28.97% in E. grisescens, whereas it was dramatically lower in E. obliqua (0.02%). Wolbachia was not detected in female E. obliqua samples (O1-O5), but were detected in 60% of male samples (O6-O10). The richness of Wolbachia ranged from 0.03 to 0.15% in three male E. obliqua samples. In E. grisescens, infection rates varied among samples, and no

688
Zhibo Wang et al. significant difference between female and male moths was found (table 2). Furthermore, we detected the infection rate of Wolbachia for samples from eight different geographical populations by amplifying wsp gene. A total of 86 samples were randomly selected for the identification of species using the CO1 gene, which were deposited in GenBank (Supplementary material table S2). The samples from Yuhang (Zhejiang) and Liyang (Jiangsu) were identified as E. obliqua and those from Xinchang (Zhejiang), Guiyang (Guizhou), Nanchang (Jiangxi), Yingde (Guangdong) and Enshi (Hubei) as E. grisescens. Additionally, samples from Langxi (Anhui) were identified as the two species (ten samples of E. obliqua and six samples of E. grisescens). The intra-specific genetic distances between those individuals were 0-0.2%, while the interspecific genetic distances were 3.4-3.8% (Supplementary material  table S3). The results of Wolbachia detection showed a Wolbachia infection rate of 0 and 100% for E. obliqua and E. grisescens, respectively (figs S1-S5). In addition, we also detected the infection rate of Wolbachia for 20 samples used for bacterial community analysis by wsp gene marker. The result showed Wolbachia were not detected in the samples with low Wolbachia infection rates (such as O6, O9 and O10) ( fig. S6).

Discussion
The sibling species of tea geometrid moths E. obliqua and E. grisescens both feed on tea leaves. In our study, we controlled the rearing conditions to eliminate the influence of food and environment to explore microbiotal differences in these species. The results showed higher microbiotal diversity for E. grisescens than E. obliqua, which may offer clues for understanding why E. grisescens has a wider distribution and greater adaptability than does E. obliqua. In general, predominant microbiota differed at the genus level between these sibling species. The predominant microbiota were Melissococcus, Staphylococcus and Enterobacter in E. obliqua but Wolbachia, Enterobacter and Pseudomonas in E. grisescens. Regarding microbiota species related to reproduction, Wolbachia abundance was significantly different between these sibling species.
Wolbachia is Gram-negative bacterium first found in the oophoron of Culex pipiens (Hedges et al., 2008). Wolbachia is mainly harboured in the cytoplasm of a host germ cell, with two transmission routes: common maternal vertical transfer to progeny (Hoffmann et al., 1990) and less common horizontal transfer among different hosts, widely broadening the host range (Huigens et al., 2000;Ahmed et al., 2016). Previous studies have reported the presence of Wolbachia in nematode species and many arthropods, with widespread infection in insecta (Hilgenboecker et al., 2010). Indeed, Wolbachia has garnered intense interest because of its ability to alter the biology of its host, especially with regard to reproduction. This bacterium can manipulate insect-host reproduction in various ways including parthenogenesis (PI), feminization, male killing (MK) and CI (Lus, 1947;Terry et al., 1997;Zhang et al., 2015;Lindsey et al., 2016). To date, PI has been found in mites, hymenopterans and thrips (Arakaki et al., 2001;Werren et al., 2008). Generally, they have a specific sex-determination system where unfertilized eggs develop into haploid males while fertilized eggs develop into diploid females (Weeks and Breeuwer, 2001). Wolbachia can induce PI, whereby females develop without fertilization (Lindsey et al., 2016). In MK, male individuals are killed during embryonic development, which is induced by multiple factors, and this process has been reported in more than 20 insects (Hurst et al., 1997), though with Wolbachia as a known cause in only Adalia bipunctata and Acraea encedon (Jiggins et al., 2010). Among all reproductive impacts, CI is the most typical and conspicuous (Zhang et al., 2015). In the simplest form, CI can be described as embryonic mortality that occurs when uninfected females mate with Wolbachia-infected males (Bourtzis et al., 1996). Most CI embryos exhibit defects in paternal chromosome condensation, resulting in paternal 'diffused chromatin' that cannot be normally distributed to the zygote during metaphase I (Uyen et al., 2006). Previous studies have reported CI in Acari, Coleoptera, Diptera, Hemiptera, Hymenoptera, Isopoda, Lepidoptera and Orthoptera (Bourtzis et al., 1996;Werren et al., 2008;Chevalier, 2012;Pinto et al., 2013;Zhang et al., 2015).
In our study, E. grisescens was found to be infected with Wolbachia, but E. obliqua showed little infection by this genus, which is in line with the CI requirement of mating between uninfected and infected individuals. Recent research has indicated that the hatching rate of the filial generation was notably decreased when female E. grisescens mated with male E. obliqua and that it was even lower when female E. obliqua mated with male E. grisescens (Xi et al., 2014;Zhang et al., 2014). Overall, our results are in accord with CI and suggest that Wolbachia might be an important factor causing reproductive isolation between these species.
Generally, the fundamental rule of distinguishing species is reproductive isolation which may result from prezygotic or postzygotic barriers (Dobzhansky, 1970). In our study, the kind of phenomenon, can mate but produce sterile offspring between the two tea loopers, belongs to postzygotic isolation (Haldane, 1922). The postzygotic isolation is a rare phenomenon relative to prezygotic in nature and exists in those species which have a close genetic relationship such as Spodoptera frugiperda (fall armyworm). Spodoptera frugiperda is polyphagous and a major agricultural pest in the North and South American continent and Caribbean (Gouin et al., 2017). It consists of two sympatric host-plant strains, C strain feeding mostly on maize cotton and sorghum and R strain mostly associated with rice and various pasture grasses (Gouin et al., 2017). These two strains are morphologically indistinguishable but estimated to be 2.09% on average in the COI gene (Kergoat et al., 2012). Compared to genetic distances of other species, Dumas et al. proposed that C and R strains were pairs of differentiated species (sister-species) in the Figure 5. Relative abundances of microbiota phyla in sibling species of tea geometrid moths. Columns of different colour represent the abundance of the top ten microbiota. Others represent the sum abundance of those microbiota not among the top ten.
Spodoptera genus (Dumas et al., 2015a(Dumas et al., , 2015b. More important, C and R strains also showed the phenomenon of postzygotic isolation (Groot et al., 2016). Though bacterial endosymbionts can cause genetic incompatibilities in hybrids, Dumas et al. investigated the presence of Wolbachia and several other bacteria in both C and R strains of S. frugiperda, but did not detect any of them, accordingly, eliminated the factor that bacteria manipulate their host reproduction (Dumas et al., 2015a). However, in our study, E. obliqua and E. grisescens have identical dietary habits and differential bacterial endosymbiont in Wolbachia. Therefore, this pair of insects constitute suitable material for research that Wolbachia mediating host reproduction isolation. The results of this study reveal the diversity of microbiota taxa between sibling species of tea geometrid moths. In particular, the notable difference in Wolbachia may be the major factor influencing the reproductive isolation of these sibling species. Overall, Wolbachia is an important microbiota genus manipulating insect-host reproduction in various ways. In addition, more functions of Wolbachia can be explored in these model sibling species of tea geometrid moths, such as Wolbachia interaction with pheromones and EoNPV. Nonetheless, further research is required to explore the mechanism by which Wolbachia is involved in reproductive isolation between sibling species of tea geometrid moths.