While our recent study indicated that echinoderms have a comparable evolutionary rate in their genomes compared with that of chordates (Li et al., 2020), it still remains to be tested if echinoderms have experienced larger modifications to their development than their sister groups. Meanwhile, it is tempting to know if developmental stages that establish pentaradial symmetry and its later stages are more highly derived than the embryos of their sister groups. To answer this question, we developed a method for quantifying the derivedness of embryos using gene expression profiles (Figure 1A). Specifically, we used whole embryonic developmental transcriptomes as a phenotype of embryos, and compared them between echinoderms (Anneissia japonica, Apostichopus japonicus, Lytechinus variegatus, Strongylocentrotus purpuratus) and chordates (Mus musculus, Gallus gallus, Pelodiscus sinensis, Xenopus laevis, Danio rerio, Oryzias latipes, Ciona intestinalis, and Branchiostoma floridae).

In contrast to approaches that compare the expression of only conserved 1:1 orthologs, our approach encompasses expression of paralogs and potentially lost genes to cover as many evolutionary changes as possible. As indicated in a previous study, 1:1 ortholog-based comparisons encounter a major barrier when multiple species are compared (Hu et al., 2017). In brief, the number of 1:1 orthologs that can be identified comes up to a very small ratio of genes in the entire genome (e.g., only 1,704 1:1 orthologs could be identified in their analysis with eight chordate species (Hu et al., 2017)). This situation is intensified in our analysis when both chordates and echinoderms are involved (13 species in total), where only 271 1:1 orthologs could be identified. These not only account for ∼1% of all genes in a typical vertebrate genome (∼20,000 genes) but also overlook the changes (such as those by gene duplication and/or gene loss) that occur during evolution, which leads to possible underestimation of how derived each embryo is from their common ancestor.

To identify genes in orthologous groups of distantly related species, we first compared protein-coding genes and identified 22,699 ortholog-groups using the PorthoMCL software for eight chordates, four echinoderms, and one outgroup species (Li et al., 2003; Tabari and Su, 2017) (Supplementary Figure S1). We then calculated the normalized expression (see also Table 1 in Methods for the normalization methods tested) for each ortholog-group by 1) taking the mean expression of paralogs and 2) giving “zero” expression value to potentially lost genes (taking sum-expression of paralogs also provided similar results, as in a previous study (Hu et al., 2017); see also Figure 1B and Supplementary Figure S3A). However, species-specific genes were not included, as the analyses including these genes did not meet the criteria we utilized (described in detail below).

Finally, transcriptomic distance was calculated with various combinations of normalization, distance, and tree inference methods and scanned for a suitable combination (Figure 2A) based on the following three criteria (Figure 2B): 1) developmental stages cluster by species (Kalinka et al., 2010; Hu et al., 2017); 2) tree topology is consistent with known phylogeny inferred from genomic sequences (suggested by previous studies (Kalinka et al., 2010; Hu et al., 2017)); and 3) within-species transcriptomic distances show gradual changes along development (high smoothness in the distance image).

Criterion 1, a topology of developmental stages clustered by species, can be expected for the species covered in this study. For example, even with six closely related Drosophila species that split less than 40 million years ago and share very similar morphological features, their whole embryonic transcriptomes still cluster by species rather than by stages (Kalinka et al., 2010). This is presumably because despite sharing many conserved developmental features, larger differences including those with respect to species-specific characteristics exist, such as ovariole number and genome size (Markow and O’Grady, 2007), and this may have clustered samples by species. Development-related genes account for only 10–25% of the whole genome of chordate species (Hu et al., 2017). Therefore, in this study, we selected representative species that should be distantly related enough to avoid the tree topology of clustering by stage. Criterion 2 was based on the assumption that the derivedness of the tree based on overall transcriptomic profiles does not differ from the known phylogenetic topology deduced from the genome, and this is supported by a study that analyzed more closely related Drosophila species (Kalinka et al., 2010). For criterion 3, we deployed image texture analysis to score the degree of smoothness of distance matrix images generated by different distance methods, as transcriptomic changes are expected to be continuous along developmental stages.

The results of our scanning show that both rank and logarithmic transformations of the expression data met criterion 1, namely, clustering samples by species (shown in Figure 2C and in the visualized expression images in Supplementary Figure S2). In addition, we found that calculating pairwise distances of gene expression profiles using “1—Spearman’s correlation coefficient ρ” (abbreviated as Spearman distance below) meets all three criteria and is the most suitable method among the tested distance methods (Figures 2C, 2D, 3A). Finally, we tested several commonly used tree inference methods, where BIONJ outperformed the other tested methods with its speed and algorithmic design to cluster samples with the lowest distances together (Saitou and Nei, 1987; Gascuel, 1997), although most methods generated similar results except for the Fitch-Margoliash method (Supplementary Figure S3B). A simplified scheme of the most suitable combination based on our criteria is summarized in Figure 2E. As mentioned above, our analyses did not cover the expression of species-specific genes. Intuitively, inclusion of species-specific genes may allow us to more comprehensively cover and evaluate the changes made during evolution; however, we excluded these genes as the analyses including species-specific genes did not meet criterion 2 (consistent with known phylogeny). X. laevis became an outgroup of other vertebrates in this tree (Supplementary Figure S4). This may be due to the low accuracy of gene prediction, which led to an excessive number of species-specific genes which biases the distance calculation (Supplementary Figure S1D). Notably, 7,879 (40%) of these non-paralogous, specific-specific genes in X. laevis are lowly expressed (max TPM <1) throughout the developmental stages examined, suggesting that many of these genes might be inaccurately annotated.

We defined the derivedness index of each embryo as the total branch length from the putative common ancestral node on the inferred tree (Figure 3A), which is the common ancestor of chordates and echinoderms. As has been assumed, some of the echinoderm species, especially embryos of two sea urchins, showed much higher derivedness indices than those of invertebrate chordates. Notably, the tendency was more obvious for the developmental stages with pentaradial symmetry (shaded stages in Figure 3A). Meanwhile, the other echinoderms, such as the feather star (an early-diverged echinoderm species (Cannon et al., 2014; Li et al., 2020)) and the sea cucumber showed rather less-derived indices when compared to those of vertebrates (Figure 3A). These tendencies were corroborated by summarizing the range of the derivedness indices of all embryos for each species [Figure 3B; changes in the indices among different species were statistically significant (p < 2e-16 by Kruskal-Wallis test)], and that for all chordates and echinoderm species. (Figure 3C). Nevertheless, when we compared the derivedness index of the penta-radial phase of echinoderm development (including metamorphosis, juvenile, and adult stages) against the bilateral phase (early stages to larval stages, before penta-radial structures start to appear) in each species, the penta-radial phase was indeed more derived in the feather star and the two sea urchin species (Figure 3D, statistically significant, tested with the Wilcoxon signed-rank test). Considering that the two sea urchin species, which split around 180 million years ago (Li et al., 2020), share similar developmental characteristics, the reason why L. variegatus showed a relatively higher derivedness index than S. purpuratus remains to be clarified (L. varietagus became even more highly derived when species-specific genes were considered; Supplementary Figure S4). However, this could partially be due to L. variegatus having a faster genomic evolutionary rate than S. purpuratus (Li et al., 2020; see also Supplementary Figure S19). These results suggest that although echinoderms may not have experienced larger modifications to their molecular developmental systems than chordates since their split from the deuterostome common ancestor, sea urchins appear to have accumulated more changes in their developmental systems, especially in the later embryonic stages when pentaradial structures become evident. However, this may not be solely due to gene expression from pentaradial structures or their source structures, as the sample “Lv_8wpfRudiment” (the forming rudiment in 8-weeks-post-fertilization embryo of L. variegatus) is almost as highly derived as “Lv_8wpfLarva” (the remaining larval body with the rudiment removed), which does not contain pentaradial structures. The differences in transcriptomic derivedness indices of the pre-metamorphosis and the pentaradial phases could at least partly be attributed to how genes are expressed at different levels in the two phases (rather than the two developmental phases deploying different sets of genes) (Supplementary Figure S17). Further studies delineating the molecular mechanisms of the metamorphosis process may help explain why stages with pentaradial structures tend to show highly derived indices.

Consistent with previous studies that supported the developmental hourglass model with 1:1 orthologs (Irie and Kuratani, 2011; Levin et al., 2012; Wang et al., 2013; Uesaka et al., 2019), our results show a similar trend in vertebrates, with conserved mid-embryogenesis (especially organogenesis stages) being the least derived within each species from the common ancestor of chordates (Figures 4A–F). Meanwhile, in some species, the identified least derived stages showed a slight shift from the previously reported conserved stages; for instance, while Prim-5-6 of the zebrafish (D. rerio) was previously suggested to be the most conserved stage (or vertebrate-phylotypic period), but the least derived stage was at a slightly earlier, 14-somites stage, and the Prim-5-6 stage was the second least derived embryo (see also Supplementary Figure S7). Larger differences were observed in Ciona and amphioxus. While the most conserved stages, identified by 1:1 orthologs, in Ciona were around st.24–29, the least derived stages were st.1–10. In amphioxus, while the mid-neurula stage (Hu et al., 2017; Marlétaz et al., 2018) was the most conserved, stage L2 (open-mouth larva) was found to be the least derived. These results imply that “conserved” phylotypic period are in general less derived than the earlier/later stages, but these may show larger differences especially when more number and/or distantly related species are being compared.

Nevertheless, among the species analyzed, for the range of derivedness indices and position of the common ancestral node of all the embryos of each species, the Chinese soft-shell turtle tends to show more derived developmental systems than the other vertebrates. This is consistent with the view that turtles possess many highly derived morphological diapsid features (Nagashima et al., 2009; Gilbert and Corfe, 2013). However, genomic analyses indicated that turtles have a slow evolutionary rate (Bradley Shaffer et al., 2013; Green et al., 2014), and further studies are required to explain this discrepancy. Meanwhile, given that the fossil records of turtles remained morphologically conserved since their appearance (Li et al., 2008; Benton, 2014), it is possible that the turtle embryos are transcriptomically more derived (expressing orthologous genes in a different way) than the other vertebrate species while maintaining a slow genomic evolutionary rate. In contrast, our results indicate that the mouse appears to possess the least derived molecular developmental system among the vertebrates compared (Figures 3B, 6B-iii). Another unexpected result was that the overall stages of Ciona showed the least derived indices when compared with those of the embryos of other chordates, including the early diverged amphioxus. This is in contrast to both morphological and genomic studies, which suggested that tunicates could be the most derived species in chordates (Holland, 2015, 2016). Moreover, a similar tendency was also corroborated by the tree that covered species-specific genes (Supplementary Figure S4) and the tree drawn by 1:1 orthologs (Supplementary Figure S5). The exact reason for this is unclear; however, a potential reason would be due to the vast number of potentially lost genes in Ciona, as it may bias the distances against different species by giving higher transcriptomic similarity to the rest of the species. Given that our method correctly captured transcriptomic derivedness, the results imply that their orthologous genes used during embryonic development remain rather ancestral, despite their highly divergent genomic sequences (Berná and Alvarez-Valin, 2014) after the split from the Olfactores common ancestor.

While there were a few reports that hint how the body plan-establishing phase became conserved (Bogdanović et al., 2016; Hu et al., 2017; Zalts and Yanai, 2017; Uchida et al., 2018), we further sought for potential hints toward the mid-embryonic conservation. Specifically, we asked which ortholog-groups could potentially characterize this conserved phase, given that the least derived stage only differed slightly from the reported conserved phase in vertebrates. For this purpose, we sought ortholog-groups that were highly/lowly expressed during the stages with less derivedness, by calculating the correlation (Spearman’s correlation coefficient) between the expression of each ortholog-group and derivedness indices during development. For example, ortholog-groups with strong negative correlations across species (dark blue in Figure 5A, abbreviated as “negative DCOs,” derivedness-correlative ortholog-groups) tended to show higher expression around the less derived developmental phase and lower expressions in more derived stages. We especially looked for negative DCOs across the six vertebrate species, as they may represent transcriptomic features of the phylotypic period and could potentially provide a hint about the mechanism of its conservation.

To avoid unexpected bias, we excluded samples obtained by Quartz-Seq (2-cells to blastocyst in mouse and all stages in sea cucumber) for this analysis, as these were obtained by different RNA-seq protocols than the others (discussed in detail in the next section). Among the 22,699 ortholog-groups analyzed, we focused on 7,775 ortholog-groups that had at least one gene counterpart in each vertebrate species (Figure 5B). Among these, we found 695 ortholog-groups showing negative correlation in all six vertebrate species concomitantly (category 6 in Figure 5B; Supplementary Table S19), including 18 HOX ortholog-groups, 201 development-related ortholog-groups, and 161 ortholog-groups involved in signaling transduction (Figure 5C, Supplementary Figures S8–10; the ratio of development-related ortholog-groups is significantly higher in category 6 than in any other category; Fisher’s exact test), consistent with the actively proceeding organogenesis in these stages. Along with whole embryonic expression, some of aforementioned genes are also known to show conserved spatial expression patterns (Thisse et al., 2004; Darnell et al., 2007; Yokoyama et al., 2009; Bowes et al., 2010; Richardson et al., 2014; Smith et al., 2018) (Figure 5C). Of note, these negative DCO genes that show spatial conservation contain many genes that are reported to be involved in neural patterning. For example, Nkx6.1, Pax2, Lmx1b, Wt1, Evx1/2, En1, Mnx1, Sox1, and Olig2/3 are involved in the dorsal-ventral patterning of neurons in the spinal cord (Catela et al., 2014). Additional examples are the two ortholog-groups with the strongest mean negative correlation. Patched1/2 and Rfx4, which are receptors of sonic hedgehog (Shh) and a downstream target of Shh signaling, respectively, whose mutants show severe defects in the neural tube (Goodrich et al., 1999; Ashique et al., 2009; Murdoch and Copp, 2010; Sedykh et al., 2018) (see also Supplementary Figure S11; Shh also showed negative correlation across species). Genes involved in other conserved structures of the phylotypic embryo were also identified as negative DCOs (Figures 5C,D), such as those involved in somitogenesis (Patched1/2, Nkd2), heart (Nkx2.5, Tbx20, Gata4), and mesonephros (Pou3f3). Besides, this list of negative DCOs also included genes that were not previously considered to be directly related to embryonic development, such as Fidgetin (the ortholog-group showing the strongest mean negative correlation of −0.798 among all ortholog-groups analyzed). Fidgetin has been recently suggested to regulate the cytoskeleton in the spinal cord and somites in the mouse embryo (Leo et al., 2015). More genes that are involved in organogenesis could potentially be identified from these DCOs.

In addition, consistent with a recent report suggesting that developmental stages with more pleiotropic (repetitively used) genes tend to be evolutionarily conserved, ortholog-groups showing negative correlation in more species tended to be expressed in more developmental stages (temporal pleiotropy, Figure 5E), suggesting that negative DCOs are potentially useful for characterizing the least derived stages. However, we note that this method may overlook ortholog-groups that have significant contribution to the phenotype of the least derived stages but show consistently high expression throughout development because they may not show a strong negative correlation with the derivedness index.

Finally, we extended the analysis to the entire chordate and echinoderm clades, which have been difficult to perform with 1:1 ortholog-based or conservation-oriented methods (as only ∼1% protein-coding genes of the entire genome could be compared). Among the 4,026 chordate ortholog-groups, 230 were identified as negative DCOs in all eight chordate species analyzed (Figure 5F, Supplementary Table S20). As expected, these contained a variety of genes that are known to be important for chordate embryogenesis, such as Tbx20 (Belgacem et al., 2011), which showed a strong average negative correlation. This implies that negative DCOs might provide hints for the identification of genes retaining conserved functions shared in chordates in future studies, however, it also has to be noted that sequence-based similarity does not always indicate their similar functions, and roles of negative DCOs have to be analyzed carefully in the future studies. Similarly, in echinoderms, we found 2,414 DCOs showing negative correlation across all three echinoderm species (left-most in Figure 5G, Supplementary Table S21), which contained genes mainly expressed in gastrula in sea urchins and 2-cell to 32-cell stage in feature star (the least derived stages shown in Figure 4I–L). Since the functions of genes in echinoderm species remain largely unexplored, we could not deduce much information from the negative DCOs in echinoderms; however, some implications were obtained for smaller groups in echinoderms, namely, sea urchins (Figure 5G). In contrast, genes involved in the establishment and growth of the pentaradial phase were included in the cluster showing positive correlation in all three echinoderm species since this developmental phase was shown to be more derived.

While many of our results correspond with previous EvoDevo studies, it has to be noted that there are several technical issues to overcome. First, the choice of outgroup species may affect the topology of the tree, as is the case in most evolutionary studies. In contrast to Figure 3A, the tree with Drosophila melanogaster as the outgroup (Supplementary Figure S13) showed a tree topology that deviated more from the known phylogeny. This could potentially be because fewer genes in D. melanogaster, (which is often considered a highly derived species in arthropods (Andrioli, 2012), could be identified as orthologous genes of deuterostomes. Second, differences in read depth may cause biases in expression quantification, as samples with deeper sequencing are expected to detect more lowly expressed genes than shallower samples. To test how these may affect the derivedness index, additional trees were plotted from expression data with depth adjusted uniformly to 10 million (randomly picked up 10 million mapped reads, 10 M; which corresponds to the depth of the shallowest sample), gradual down-sampling of read depths, or adjusted proportionally to the exome size of each species. In the tree with the 10 M expression data, the range of derivedness indices for six species, notably Ciona, increased significantly (Figures 6A–i, see also Supplementary Figures S14, S16), while the tree with exome sizes-adjusted depths showed significant changes in fewer species (Supplementary Figure S15). Notably, the topologies of both trees were inconsistent with known phylogeny. This suggests that although the derivedness index tends to be influenced by read depth, the influence may be dependent on the species or the samples; thus, more comprehensive studies are needed to determine which read depth normalization method would be more suitable for measuring transcriptomic derivedness. Finally, as we included samples acquired by different library preparation protocols where the Quartz-Seq protocol involved an additional whole-transcript amplification step (Sasagawa et al., 2013), we asked whether these samples were comparable to each other. To answer this, we collected new RNA-seq data from E9.0 and E15.5 mouse embryos by performing both TruSeq and Quartz-Seq from the same starting total RNA sample. A tree with all samples, including biological replicates, was plotted. Unexpectedly, samples tended to cluster by protocols rather than by stages or biological replicates (Figures 6B-iv,v), suggesting that samples obtained by these two protocols may not be completely comparable. To avoid this bias, we excluded samples acquired by Quartz-Seq for the analysis of derivedness-correlative ortholog-groups.

Experimental procedures and animal care were conducted in strict accordance with the guidelines approved by the University of Tokyo (approval ID: 14-03, 16-2). The animals were sacrificed with minimal suffering. Individual embryos were blindly selected from the wild-type population.

The RNA-seq data utilized were from published datasets (Wang et al., 2013; Hu et al., 2017; Ichikawa et al., 2017; Li et al., 2018, 2020) and three other studies [purple sea urchin (Tu et al., 2012); oyster: (Zhang et al., 2012); Drosophila: (Nègre et al., 2011):]. These include major early to-late developmental stages of eight representative chordate species (Mus musculus, Gallus gallus, Pelodiscus sinensis, Xenopus laevis, Danio rerio, Oryzias latipes, Ciona intestinalis, Branchiostoma floridae), four echinoderms (Anneissia japonica, Apostichopus japonicus, Lytechinus variegatus, and Strongylocentrotus purpuratus), and two outgroup species (Crassostrea gigas, Drosophila melanogaster), with 2-3 biologically independent replicates of each stage to represent the general population (except S. purpuratus, C. gigas, and D. melanogaster). All samples were sequenced using Illumina platforms. Details of the included datasets are summarized in Supplementary Tables S1–13 (developmental stages covered) and Supplementary Table S14 (e.g., sample accession numbers, sequencing platforms).

To compare datasets collected by TruSeq and Quartz-Seq (Sasagawa et al., 2013) protocols (the latter involves an additional whole-transcript amplification step), additional samples of mouse E9.0 and E15.5 (C57BL/6J strain, CLEA Japan) were prepared following the same procedures reported in the mouse and sea cucumber datasets (Hu et al., 2017; Li et al., 2018). Both TruSeq and Quartz-seq were performed on the same total RNA extracted from pooled embryos. An average of nine embryos from three independent parents were pooled for each E9.0 sample, and an average of three embryos, E15.0. Library preparation for sequencing (TruSeq: single-end, 100 bp, non-strand-specific; Quartz-Seq: paired-end, 100 bp, non-strand-specific) was performed according to the manufacturer’s instructions using the Illumina TruSeq RNA Library Prep Kit v2 or the standard protocol of Quartz-Seq (Sasagawa et al., 2013) with Illumina Nextera XT Library Prep Kit. The normalized libraries were sequenced using an Illumina HiSeq 1,500. The raw reads data are available in SRA through the accession number PRJNA749373.

Adapter sequence trimming and low-quality read filtering in RNA-seq raw reads were performed with trimmomatic (version 0.36) with default parameters (ILLUMINACLIP:2:30:10, LEADING:3, TRAILING:3, SLIDINGWINDOW:4:15, and MINLEN:36) (Bolger et al., 2014). Quality checks were performed using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were aligned to the respective genome using HISAT2 (version 2.1.0) (Kim et al., 2019), with mitochondrial scaffolds and NUMT sequences manually removed from the genomic sequence files to avoid multi-copy quantification biases. The genome versions are summarized in Supplementary Table S15. The expression levels of coding genes were quantified in TPM (transcripts per million) (Wagner et al., 2012) using StringTie (version 1.3.4d) (Pertea et al., 2015).

Gene-based expression profiles of embryos of all species were summarized into an ortholog-group-based expression table by the method explained in the main text and Figure 1B. Our method was modified from the approach previously reported by the EXPANDE Project Consortium (Hu et al., 2017). Genes were first grouped into ortholog-groups using PORTHOMCL (Li et al., 2003; Tabari and Su, 2017) with a BLASTP (version 2.7.1) e-value threshold of 1e-5, as an approach considering paralogs and potentially lost genes. The expression level per ortholog-group was then estimated by the mean (or sum) expression level of all genes in the ortholog-group in each species. The expression level of “0” was assigned to an ortholog-group with no predicted gene in that species.

The branch length of each embryo from the common ancestral node on the tree (putative common ancestor of chordates and echinoderms) was extracted and defined as the derivedness index. In calculating derivedness indices for species as a whole (for whole life cycles), the range of derivedness indices of all the developmental stages in the inferred tree (Figure 3A) was represented as box plots (whether the differences among species were statistically significant were examined by Kruskal-Wallis test; Figure 3B). To compare the derivedness indices of penta-radial phase against pre-metamorphosis stages in each echinoderm species (Figure 3D), the range of the derivedness indices of the corresponding developmental phase was represented as a box plot (Mann-Whitney U test). Penta-radial phase for each species was defined as: (feather star) doliolaria onwards; (sea cucumber) metamorph-1 onwards; (green sea urchin) 7wpf onwards; (purple sea urchin) vestibular-invagi onwards. To examine the derivedness index of each developmental stage (Figure 4), the results from the 100 random BRI trees were utilized, and for each developmental stage, the range of derivedness indices in the 100 trees was summarized for each embryo (whether the differences among developmental stages were statistically significant were examined by Friedman test).

Samples collected by Quartz-Seq were omitted from the analysis to avoid any unexpected bias. For each species and ortholog-group, Spearman’s correlation coefficient between its expression levels and derivedness indices along early to late embryos was calculated. Ortholog-groups showing negative correlations across certain groups of animals (vertebrates, chordates, and echinoderms) were further analyzed. A total of 7,775 DCOs, with at least one gene counterpart in each vertebrate species, were classified into seven categories based on the number of vertebrate species that showed a negative correlation.

To normalize read depth between samples of different species, the maximum possible number of reads was selected while keeping the same reads-to-exome-size ratio in all species. Best-hits were first selected using samtools (version 1.12; single-end samples: samtools view -F 4 | samtools view -F 256; paired-end samples: samtools view -f 2 | samtools view -F 256) (Li et al., 2009). Random pick-up of reads was performed from the BED file containing all mapped best-hit reads. For paired-end samples, only one read was reported for the BED file using the bamtobed function with the “-bedpe” option from bedtools (version 2.30.0) (Quinlan and Hall, 2010). The ratio of all best-hit reads to exome size in each sample was calculated, and the smallest ratio was selected as the baseline (Lv_60cell sample) to calculate the number of reads needed for each species (see Supplementary Table S18 for the calculated results). We tried to calculate exome size using the original genome annotation file as well as a version with all the untranslated regions (UTRs) removed to minimize unintentional bias, since only some of the species were annotated with UTRs (Supplementary Figure S15). To remove UTRs from annotation files, the script “gff3_file_UTR_trimmer.pl” from PASApipeline (Haas et al., 2003) was used. To calculate exome size, the command line “awk '{if ($3 = = "exon"){print $0}}' $gffname | gff2bed - | bedops -m - | awk ‘BEGIN{FS = "\t"; count = 0}{count = count + ($3-$2)}END{print count}'” with the BEDOPS tool (version 2.4.39) (Neph et al., 2012) was used.

Bioinformatics analyses were performed using in-house R (4.0.3) (R-Core-Team, 2020), Python (3.7), and shell scripts, together with the software and packages summarized in Supplementary Table S17.