The RTBV sequence (NCBI accession: NC_001914.1) was used as an outgroup for several analyses. For eRTBVL families, the consensus sequences of eRTBVL-A, -B, and -C were obtained from NCBI (NCBI accessions BR000029.1, BR000030.1, and BR000031.1, respectively) and the consensus sequence of eRTBVL-X used in this study was as described by Chen et al. (2018). According to Chen et al. (2018), collections of eRTBVL-D were employed with a BLASTn search of the O. sativa japonica genome using the consensus sequences of six other eRTBVL groups as queries. The segments with the highly reliable hits (e-values < 1 × 10⁻³ and lengths > 100 bp) and showing <85% sequence identities to non-eRTBVL-D consensus sequences were collected as eRTBVL-D”.

The eRTBVL-D loci mapped to the Nipponbare genome (O. sativa ssp. japonica) were used as the references. The genome sequences are available for the seven AA-genome species O. sativa ssp. japonica and indica, O. rufipogon, O. glaberrima, O. barthii, O. glumaepatula, O. meridionalis, and O. longistaminata ( Supplementary Table 1 ). These genome sequences were mainly obtained through Gramene (http://gramene.org/) and NCBI (https://www.ncbi.nlm.nih.gov/Traces/wgs/?view=wgs) except for indica (V2plus: https://ars.els-cdn.com/content/image/1-s2.0-S1674205217300424-mmc9.zip) (Zhang et al., 2016). The genome sequence of O. punctata, which is an Oryza BB-genome species, was also used to examine whether eRTBVL-D sequences are present in the databases of Gramene and NCBI ( Supplementary Table 1 ). The BB genome species are considered the closest group genetically to the AA genomic species, with O. punctata being the BB genome species, and O. punctata is the plant species whose genome has been analyzed the most among the BB genome species (Kim et al., 2007). We selected the candidates for eRTBVL-D orthologous sequences in each genome of the AA-genome species when the corresponding eRTBVL-D in the Nipponbare genome possessed identical contiguous flanking sequences. The selected candidate sequences were confirmed with a highly reliable hit (e-values < 1 × 10^-3 and lengths > 100 bp) with ones in the Nipponbare genome and were judged as orthologous (Chen et al., 2018). To retrieve the sequences of flanking genes for each locus, the genes closest to an eRTBVL-D locus in the Nipponbare genome and the orthologous genes in the AA-genome species were collected from Gramene ( Supplementary Table 2 ). The second flanking gene was used for analysis when the closest ortholog was found in less than five AA-genome species.

The phylogenetic relationships of five eRTBVL families, A, B, C, D, and X, were analyzed by constructing phylogenetic trees for three predictable regions in the eRTBVL sequence, a protease (PR) region, a reverse transcriptase/RNase H (RT/RH) region, and an ORFz region ( Figure 1A ). A set of the orthologous sequences within the AA-genome species, such as an eRTBVL-D locus, its flanking sequences, and the closest genes, were used for phylogenetic analysis. The five RTBVL-D sequences from d2, d8, d9, d12, and d13 were omitted from the phylogenetic analysis data sets due to insufficient length or a lack of species carrying the segment ( Supplementary Table 1 ). Each of three eRTBVL-D-linked sequences, d4 to d5 (d4-5), d6 to d7 (d6-7), and d10 to d11 (d10-11), which were contained on intervals shorter than 1000 bp, was analyzed as a single locus. In these three loci, the adjacent eRTBVL-D sequences resided on different DNA strands ( Supplementary Table 1 ). The sequence alignments of seven data sets (d1, d3, d4-5, d6-7, d10-11, d14, and d15) were generated in MUSCLE (Edgar, 2004), and all gaps were removed. The alignments of the flanking gene sequences were generated with E-INS-i algorithm of MAFFT v.7.490 (Katoh et al., 2019), which is suitable for alignments containing large gaps. Maximum likelihood (ML) phylogenetic analysis was performed, and the best-fitting substitution models for each locus were determined by model selection analysis. ML trees were established based on the support of 1000 bootstrap replicates. All analyses related to phylogenetic analysis were completed using MEGA version 7.0 (Kumar et al., 2016). Evolutionary analyses were conducted in MEGA X (Kumar et al., 2018; Stecher et al., 2020). The numbers of base substitutions per site and variances (analytical method option) were calculated using the Jukes-Cantor model (Jukes and Cantor, 1969; Kimura and Ohta, 1972). All positions containing gaps and missing data were eliminated (complete deletion option).

Fifteen eRTBVL segments found in the Oryza sativa ssp. japonica (cv. Nipponbare) genome were identified as members of the eRTBVL-D family. These sequences were distinguished from those of previously identified eRTBVL-A, -B, -C, and -X families (Chen et al., 2018). The 15 segments did not contain the entire sequence of the eRTBVL genomes, but corresponded partly or fully to the functional regions in eRTBVL as follows: d2 to ORFx; d13 to coat protein (CP); d3, d14, and d15 to PR; d3, d10, d11, d14, and d15 to reverse transcriptase, RNase H (RT/RH), and ORF4/ORFz; and d1, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d14, and d15 to the intergenic region (IGR) ( Figure 1A ). To examine genetic distances among the eRTBVL families, we constructed phylogenetic trees for the PR, RT/RH, and ORF4/ORFz regions using the consensus sequences from the five families and the RTBV representative sequence ( Figure 1B ). Each phylogenetic tree in these three functional regions showed that eRTBVL-D was present in the Oryza genomes prior to the integration of the other four eRTBVL families. As Chen et al. (2014; Chen et al., 2018) reported, the phylogenetic tree structures of the eRTBVL-A, -B, -C, and -X families were topologically closely related. The similar structures of the three phylogenetic trees in the PR, RT/RH, and ORFz regions supported the idea that eRTBVL-D sequences existed in the ancestral genome prior to the insertions of the other eRTBVL family members. At four of the loci containing the 15 eRTBVL-D segments (on chromosomes 1, 2, 4, and 7), two or three segments were located within about 1000 bp of the japonica genome (cv. Nipponbare) ( Supplementary Figure 1 , Supplementary Table 1 ). We regarded each of these four clusters as an eRTBVL-D locus. Ten positions in the Nipponbare genome were designated as eRTBVL-D loci ( Supplementary Figure 1 ).

Figure 2 ; Supplementary Table 1 summarize characteristics of the orthologous segments of the 15 eRTBVL-D sequences in the nine AA- and BB-genome species and subspecies studied (O. sativa ssp. japonica and indica, O. rufipogon, O. glaberrima, O. barthii, O. glumaepatula, O. meridionalis, O. longistaminata, and O. punctata). The d2 sequence was missing in O. glaberrima, O. barthii, O. glumaepatula, and O. longistaminata. In O. longistaminata, d4 and d5 were lost, and it was unclear whether d10 and d11 were present ( Figure 2 ). O. meridionalis also contained four ambiguous sequences possibly corresponding to d4, d5, d9, and d15. O. punctata, which is the Oryza BB-genome species most closely related to the AA-genome species, contained 3 of the 15 eRTBVL-D segments ( Figure 2 ; Supplementary Table 1 ). No corresponding sequences were present in a more distant species, O. brachyantha (Oryza FF-genome) (Chen et al., 2018). Taken together, the patterns of eRTBVL-D distribution in the Oryza species strongly support the idea that eRTBVL-D integrations occurred immediately before the AA-genome speciation.

The putative eRTBVL-D sequence was considered to comprise genes or open reading frames (ORFs) similar to those of the other eRTBVL families, which were found in the PR, RT/RH, ORFz, and IGR regions ( Figure 1A ). To examine the phylogenetic relationships of the orthologous eRTBVL-D segments in the AA-genome species, we constructed phylogenetic trees of the eRTBVL-D segments in each of the four representative regions containing such segments. The phylogenetic trees from the PR (d14 and d15), RT/RH (d3, d11, d14, and d15), ORFz (d3, d10, d11, d14, and d15), and IGR (d1, d3, d5, d6, d7, d8, d9, d10, d11) regions displayed mostly clusters consisting of a single orthologous segment ( Supplementary Figure 2 , regions a-f), whereas in the phylogenetic trees of the IGR region, d1/d3 ( Supplementary Figure 2D ) and d10/d11 ( Supplementary Figures 2D–F ) clusters contained different segments. In the d10/11 cluster, the d11 segment from O. meridionalis was located in the d10 cluster in the phylogenetic trees for regions d to f in the IGR. The d1/d3 cluster was a mixture of both segments, suggesting that the structural alterations during the AA-genome speciation occurred in the vicinity of the d1-3 locus ( Supplementary Figure 2D ). The five eRTBVL-D segments were excluded from the phylogenetic analyses because they were short and missing the orthologous segments from four or more species.

Zhu et al. (2014) estimated the genetic distances among the AA-genome species based on the sequences of 53 nuclear genes and 16 intergenic regions. Here, we used these genetic distances and phylogenetic relationships in the eight AA-genome species as the standard speciation pattern (SSP). We compared phylogenetic profiles of eRTBVL-D segments with the SSP of the AA-genome species (Zhu et al., 2014). Two eRTBVL-D loci, d4-5 and d6-7, contained less than 1000 bp of sequence between two eRTBVL-D segments, and each locus contained eRTBVL-D segments on different strands ( Supplementary Table 1 ). Each phylogenetic analyses of these two eRTBVL-D loci was performed with the two eRTBVL-D segments containing these interspersed sequences between the segments. We constructed phylogenetic trees of eight eRTBVL-D segments, d1, d3, d4-5, d6-7, d10, d11, d14, and d15, and determined that these trees were classified into two types. Trees constructed from Type A segments (d4-5, d10, d11, and d15) were topologically consistent with the SSP for the eight AA-genome species ( Figures 3A, B , left). In contrast, analysis of Type B segments (d1, d3, d6-7, and d14) gave rise to a variety of phylogenetic trees for AA-genome species that differed from the SSP ( Figures 3A, B , right) and that also varied among the four eRTBVL-D segments. The phylogenetic trees derived from d1 and d3 showed that Asian species were in the same clade as the Oceania species O. meridionalis ( Figures 3B, d1, d3 ), whereas the tree from d6-7 identified a species of African origin, O. longistaminata, as a sister branch of O. sativa ssp. japonica with a high bootstrap value ( Figures 3B, d6-7 ).

The d14 tree consisted of two clades, each of which contained distantly related species: one for O. barthii, O. glumaepatula, O. glaberrima, O. sativa ssp. japonica, and O. meridionalis, and the other for O. sativa ssp. indica, O. longistaminata, and O. rufipogon ( Figures 3B, d14 ). Unlike the d6-7 trees, those for d14 indicated low bootstrap values for the relationships among the orthologous sequences due to small numbers of nucleotide polymorphisms ( Figures 3B, d14 ). The length of the sequences that could be analyzed was approximately 350 bp for d6-7, which has a sufficient bootstrap value, whereas the d14 sequence was 106 bp but had a lower bootstrap value due to the presence of fewer nucleotide substitutions ( Supplementary Table 1 ). The phylogenetic trees from the four Type B eRTBVL-D segments formed specific diverging branches that were apparently inconsistent with the SSP ( Figure 3 ). The structure of these four trees caused us to question whether these varied forms were dependent on each specific eRTBVL-D sequence or on the individual chromosomal regions surrounding those sequences.

Because it is likely that all the eRTBVL-D segments were inserted into the ancestral genome immediately before AA-genome speciation, we observed variations in the genetic distances of different eRTBVL-D segments located on diverse chromosomal regions. We used the phylogenetic tree for the SSP for relative comparison of genetic distances of eRTBVL-D segments between AA-genome species. Type A eRTBVL-D segments (d4-5, d10, d11, and d15) were considered appropriate for evaluating genetic distances ( Figure 3B , left) because, unlike Type B segments, they did not have aberrant phylogenetic relationships. In each of the four Type A segments, we determined the genetic distances between O. sativa ssp. japonica and each of the six other species (O. sativa ssp. indica, O. rufipogon, O. glaberrima, O. barthii, O. glumaepatula, and O. meridionalis) ( Figure 4 ).

Figure 4 shows the genetic distances of the segments d4, d5, d10, d11, and d15 of Type A eRTBVL-D for each AA-genome species relative to japonica. The genetic distances in the AA-genome species analyzed by Zhu et al. (2014) were considered as the standard genetic distances of the SSP. The genetic distances among the five segments detected in japonica vs. glaberrima and japonica vs. barthii were different relative to Zhu’s distances. All the genetic distances of the five segments in japonica vs. rufipogon were smaller than Zhu’s distances. In the remaining three species comparisons vs. japonica, the genetic distances were more or less similar to Zhu’s distances. When comparing the genetic distances between each segment and Zhu’s distance, those for d10 and d11 differed from Zhu’s distance, with the exception of japonica-glumaepatula. In contrast, the genetic distance of d15 was not detectably different from Zhu’s distance. These results indicate that the genetic distances vs. japonica tend to be characteristic of each species to some extent. We found segments such as d15 in eRTBVL-D that evolved with a similar genetic rate as SSP, whereas others, such as d10 and d11, had genetic distances that differed from Zhu’s distances. This demonstrated the existence of chromosomal sites with specific, distinct evolutionary rates of change.

To test whether the adjacent genes exhibited the same evolutionary patterns as the eRTBVL-D sequences, we constructed phylogenetic trees using the genes flanking each eRTBVL-D sequence. Because the eRTBVL sequences of d1, d2, and d3 were clustered on chromosome 1 ( Supplementary Figure 1 ), we analyzed all three as a single locus. The flanking gene sequences were obtained from the regions upstream of d1 and downstream of d3 in japonica. This was also the case for d10 and d11 on chromosome 7 ( Supplementary Figure 1 ), which we considered as a single region for the phylogenetic analysis. Detailed information for the flanking genes is shown in Supplementary Table 2 .

The phylogenetic trees for the genes flanking d4-5, d10-11, and d15, classified in Type A, showed a similar phylogenetic relationship with the SSP for the AA-genome species ( Supplementary Figure 3 ). The flanking genes surrounding the d1-3 locus in Type B exhibited phylogenies different from those of the corresponding eRTBVL-D segments. For example, the d1-3 segment may have exclusively differentiated during AA-genome speciation relative to the neighboring regions ( Figures 5, d1-3 ). The eRTBVL-D segments of d1-3 contained sequences homologous to each other, and these might represent paralogous relationships. The flanking genes of the Type B loci d6-7 and d14 exhibited similar phylogenetic trees with respect to their eRTBVL-D segments, but these also were not consistent with that of the SSP ( Figures 5, d6-7, d14 ). This suggests that each of the chromosomal regions containing the d6-7 and d14 segments is likely to have evolved independently from the SSP. As observed for the eRTBVL-D segments themselves ( Figure 3 ), the bootstrap values were relatively high for the flanking genes of d6-7 but were relatively low for the genes near d14 ( Figure 5 ).

Chromosomal rearrangements such as insertions/deletions (indels) and translocations have been found to be associated with some eRTBVL-D segments, and these are presumed to have occurred during speciation of the AA-genome species. We compared the chromosomal structures around eRTBVL-D loci of the japonica (Nipponbare) genome with those of the other AA-genome species. Chromosome 1 in japonica retained three eRTBVL sequences, d1, d2, and d3, and these have three different alignments among the AA-genome species. Although the d2 sequence between d1 and d3 was present in O. sativa, O. rufipogon, and O. meridionalis, it was absent from the genomes of O. glaberrima, O. glumaepatula, and O. longistaminata ( Figures 2 , 6 ). In O. barthii, only d3 is present, and the d1 and d2 segments have been lost ( Figures 2 , 6 ). Notably, O. glaberrima shared phylogenetic structure not with its direct ancestor, O. barthii, but instead with O. longistaminata. Both O. barthii and O. longistaminata originated from Africa, even though they did not originate from the same lineage according to the SSP ( Figure 3A ). These rearrangements of the d1-3 segments could have caused an alteration of the phylogenetic lineage from the SSP of the AA-genome species.

Most of the AA-genome species possessed a similar structure in the region surrounding the d10-11 segment present in chromosome 7. An interchromosomal rearrangement, i.e., a translocation of the region containing d10-11 and a gene adjacent to chromosome 1, was observed for O. glumaepatula ( Figure 6 ). This translocation contained the 350- to 400-kbp region including the two eRTBVL-D segments ( Figure 6 ). The presence of the eRTBVL-D segments substantially facilitated our detection of the genomic rearrangements. These local chromosomal rearrangements occurred in the eRTBVL-D sequences and might have occurred in other chromosomal regions in the AA-genome species. They can be discovered unless the sequences containing them share the same origin. In regard to phylogenetic relationship, the d10-11 segments were grouped into Type A along with the SSP, suggesting that the translocation did not affect the phylogenetic relationships of d10-11 segments among the AA-genome species ( Figures 3B, d10, d11 , 6 ).

eRTBVL-D is a group of viruses that were active in the early stages of differentiation of the AA-genome species of the genus Oryza, as evidenced by the common distribution of EVEs in the genomes of all AA-genome species (Chen et al., 2014; Chen et al., 2018). The O. punctata genome, which is the Oryza BB-genome species most closely related to the AA-genome group, contains three possible eRTBVL-D segments (d4, d5, and d12) of the 15, but we found no clear homologous sequences for the other 12 segments ( Figure 2 ). As it is likely that the integration of eRTBVL-D occurred in an active period of virus infection (Chen et al., 2018), each eRTBVL-D segment may have been inserted into the ancestral AA-genome species just prior to its differentiation into distinct species and thus was inherited by all current AA-genome species.

Because the other eRTBVL families are present only in the genomes of O. sativa and O. rufipogon, eRTBVL-D is probably the ancestral virus of the entire eRTBVL family (Chen et al., 2018). eRTBVL-A, -B, -C, -X, and -D are more closely related to each other than to the current active virus, RTBV ( Figure 1B ) (Chen et al., 2014). RTBV infects rice via the green leafhoppers. If eRTBVL-D was a similar virus, it was possible that it was transmitted via vector-borne insects. Previous analyses of the eRTBVL families have shown that insertions of each family into the genome occurred almost simultaneously during their active periods (Chen et al., 2014; Chen and Kishima, 2016). Based on their structural features, Chen et al. (2014) identified the evolutionary timing and order of insertion of each family of eRTBVL into O. sativa ssp. japonica, indica, and O. rufipogon. Their results indicated that each family was inserted during a limited period of the speciation process. Geering et al. (2014) newly classified Florendovirus from the Caulimoviridae, to which the eRTBVL family belongs, and associated the occurrence of Florendoviruses with the lineages of the host plant genus. Relationships between speciation and virus endogenization have been reported in Musa (Ndowora et al., 1999), tomato (Strasburg et al., 2012), yam (Seal et al., 2014), Dahlia (Almeyda et al., 2014), and Citrus (Yu et al., 2019) species, with virus generation and endogenization being correlated with speciation in each case.

Multiple genome databases are available for japonica and indica, while for the other AA-genome species, few sequences other than those used here are publicly available. The present study aims to investigate whether the local chromosomal regions had specific evolutionary movements. Therefore, it is conceivable that using multiple genomic data from japonica and indica would provide similar answers for the respective eRTBVL-D copies in speciation.

As shown in Figure 4 , the genetic distances between the Type A segments of japonica and other species varied among AA-genome species. These results also implied that the evolutionary rates varied depending on the different positions of the segments in the respective chromosome. As reported by Chen et al. (2018), eRTBVL segments are assumed to have undergone neutral evolution after insertion; therefore, the genetic distances observed can be considered to directly reflect the mutation rates. Figure 4 and Supplementary Table 3 show the genetic distances between japonica and the other AA-genome species for the five Type A segments, d4, d5, d10, d11, and d15. For comparison with neutrally evolving sequences of corresponding species, standard genetic distances for synonymous, noncoding, and intergenic regions analyzed by Zhu et al. (2014) are also shown in Figure 4 . The different genetic distances among Type A loci in the same combination with japonica suggested the existence of variation in mutation rates for the different genomic regions. Relative to Zhu’s distances (2014), the genetic distances for each segment between japonica and other species also showed inconsistency, implying a variation in mutation rates for different lineages. In this way, eRTBVL segments can be utilized as evolutionary indicators of variations in local chromosomal regions.

Phylogenetic relationships of the Type B segments have different topologies from the SSP ( Figure 3B , right). Phylogenetic incongruence can be caused not only by biological factors such as violations of orthology due to lineage sorting, hidden paralogy, and horizontal gene transfer but also by analytical factors such as violation of assumed evolutionary models (Som, 2015). Genomic factors in the violations were often related to major biological changes associated with evolutionary movements, such as introgression (Cronn et al., 2003; Eaton and Ree, 2013), genetic drift (Tomiuk et al., 1998), meiotic drive (Presgraves, 2010), adaptation (Strasburg et al., 2012; Bamba et al., 2019), and horizontal transmission (MacLeod et al., 2005; Bailly et al., 2007). The phylogenetic data collected by Zhu et al. (2014) has revealed more than 20 topologies in cases where they used 53 nuclear genes and 16 intergenic regions among the AA-genome species. Thus, in the process of analyzing phylogenetic trees, stochastic incongruence of a topology might be caused by rapid radiation of genes or intergenic regions. These factors could explain what we observed in the phylogenetic relationships of Type B segments (d1, d3, d6-7, and d14) ( Figure 3B , right), as each of the Type B segments exhibited some of the characteristics associated with phylogenetic incongruence.