The impact of local genomic properties on the evolutionary fate of genes

  1. Yuichiro Hara  Is a corresponding author
  2. Shigehiro Kuraku
  1. Research Center for Genome & Medical Sciences, Tokyo Metropolitan Institute of Medical Science, Japan
  2. Molecular Life History Laboratory, Department of Genomics and Evolutionary Biology, National Institute of Genetics, Japan
  3. Department of Genetics, Sokendai (Graduate University for Advanced Studies), Japan
  4. RIKEN Center for Biosystems Dynamics Research, Japan

Abstract

Functionally indispensable genes are likely to be retained and otherwise to be lost during evolution. This evolutionary fate of a gene can also be affected by factors independent of gene dispensability, including the mutability of genomic positions, but such features have not been examined well. To uncover the genomic features associated with gene loss, we investigated the characteristics of genomic regions where genes have been independently lost in multiple lineages. With a comprehensive scan of gene phylogenies of vertebrates with a careful inspection of evolutionary gene losses, we identified 813 human genes whose orthologs were lost in multiple mammalian lineages: designated ‘elusive genes.’ These elusive genes were located in genomic regions with rapid nucleotide substitution, high GC content, and high gene density. A comparison of the orthologous regions of such elusive genes across vertebrates revealed that these features had been established before the radiation of the extant vertebrates approximately 500 million years ago. The association of human elusive genes with transcriptomic and epigenomic characteristics illuminated that the genomic regions containing such genes were subject to repressive transcriptional regulation. Thus, the heterogeneous genomic features driving gene fates toward loss have been in place and may sometimes have relaxed the functional indispensability of such genes. This study sheds light on the complex interplay between gene function and local genomic properties in shaping gene evolution that has persisted since the vertebrate ancestor.

Editor's evaluation

The study provides a fundamental understanding of the driving forces behind gene losses in genome evolution and connects the propensity for gene losses to local genomic features like mutation rate and expression pattern. The methodology is compelling, as it identifies "elusive human genes" through independent gene losses in at least two mammalian lineages. The comparative genomics and statistical analyses are thorough and rigorous, making this study appealing to readers interested in exploring the global patterns and underlying mechanisms of gene fate evolution across the phylogenetic tree.

https://doi.org/10.7554/eLife.82290.sa0

Introduction

In the course of evolution, genomes continue to retain most genes with occasional duplications, while losing some genes (Blomme et al., 2006; Fernández and Gabaldón, 2020; Shen et al., 2018). This retention and loss can be interpreted as gene fate; genes are stably retained in the genome, but some factors may cause them to transition to a state where deletion occurs. Accordingly, identification of the factors allowing gene loss may facilitate our understanding of gene fate. Gene retention or loss has generally been considered to depend largely on the functional importance of the particular gene from the perspective of molecular evolutionary biology (Albalat and Cañestro, 2016; Bartha et al., 2018; Blanc et al., 2012; Liu et al., 2015; Olson, 1999; Sharma et al., 2018; Shen et al., 2018). Genes with indispensable functions have usually been retained with highly conserved sequences in genomes through rapid elimination of alleles that impair gene functions (Hirsh and Fraser, 2001; Krylov et al., 2003; Miyata et al., 1980; Pál et al., 2006). On the contrary, genes with less important functions are likely to accept more mutations and structural variations, which can degrade the original functions, leading to gene loss through pseudogenization or genomic deletion (Jordan et al., 2002; Yang et al., 2003). To date, gene loss has been imputed to the relaxation of functional constraints of individual genes. Gene loss has further been revealed to drive phenotypic adaptation in various organisms (Albalat and Cañestro, 2016; Olson, 1999), as well as in a gene knockout collection of yeasts in culture (Giaever and Nislow, 2014; Maclean et al., 2017).

To uncover the association between fates and functional importance of the genes, molecular evolutionary analyses have been conducted at various scales, from gene-by-gene to genome-wide. A number of studies have revealed that the genes with reduced non-synonymous substitution rates (or KA values) and ratios of non-synonymous to synonymous substitution rates (KA/KS ratios) are less likely to be lost (Jordan et al., 2002; Yang et al., 2003). A genome-wide comparison of duplicated genes in yeast revealed larger KA values for those lost in multiple lineages than those retained by all the species investigated (Byrne and Wolfe, 2007). Other comprehensive studies of gene loss across metazoans and teleosts revealed that the genes expressed in the central nervous system are less prone to loss (Fernández and Gabaldón, 2020; Roux et al., 2017). These observations again suggest that gene fate depends on the functional constraints of a particular gene.

Besides functional constraints, several studies have identified the genes lost independently in multiple lineages, revealing that the genomic regions containing these genes ‘prefer’ particular characteristics associated with structural instability (Cortez et al., 2014; Hughes et al., 2012; Lewin et al., 2021; Maeso et al., 2016). In mammals, tandemly arrayed homeobox genes derived from the Crx gene family were lost in multiple species (Lewin et al., 2021; Maeso et al., 2016). The findings suggest that genomic features containing tandem duplications facilitate unequal crossing over, leading to frequent gene loss. Mammalian chromosome Y, which contains abundant repetitive elements and continues to reduce in size, has lost a considerable number of genes (Cortez et al., 2014; Hughes et al., 2012). In the stickleback genome, a Pitx1 enhancer was independently lost in multiple lineages inhabiting freshwater due to its genomic location in a structurally fragile site, leading to recurrent loss of pelvic fins (Xie et al., 2019). Genes and genomic elements in such particular regions may be prone to loss in a more neutral manner than the relaxation of functional importance or via functional adaptations. Accordingly, these studies focusing on the particular genomic regions led us to search for the common features in genomes that potentially facilitate gene loss. Genome-wide scans have revealed heterogeneous distributions of a variety of sequence and structural features so far, for example, base composition (Bernardi and Bernardi, 1986; Cohen et al., 2005; Katzman et al., 2011), the frequency of repetitive elements (Korenberg and Rykowski, 1988; Medstrand et al., 2002), and DNA-damage sensitivity induced by replication inhibitors (Debatisse et al., 2012; Helmrich et al., 2006). However, the extent to which these characteristics are associated with gene fates has not been understood well at a genome-wide level.

The accumulation of near-complete genome assemblies for various organisms facilitates comprehensive taxon-wide analysis of gene loss (Fernández and Gabaldón, 2020; Guijarro-Clarke et al., 2020; Rice and McLysaght, 2017). Along with this motivation, we recently performed a comprehensive analysis on the fate of paralogs generated via the two-round whole-genome duplications in early vertebrates (Hara et al., 2018a). The results revealed that the genes retained by reptiles but lost in mammals and Aves rapidly accumulated not only non-synonymous but also synonymous substitutions in comparison with the counterparts retained by almost all the vertebrates examined, indicating that those genes prone to loss show increasing mutation rates. Furthermore, these loss-prone genes were located in genomic regions with high GC contents, high gene densities, and high repetitive element frequencies. These findings suggest that the fates of those genes are influenced not only by functional constraints but also by intrinsic genomic characteristics. Because the findings were restricted to a set of particular genes, they prompted us to examine whether this trend is associated with gene fates on a genome-wide scale.

In this study, we inferred molecular phylogenies of vertebrate orthologs to systematically search for the genes harboring different fates in the human genome. We previously referred to the nature of genes prone to loss as ‘elusive’ (Hara et al., 2018a; Hara et al., 2018b). In this study, we define the elusive genes as those that are retained by modern humans but have been lost independently in multiple mammalian lineages. As a comparison of the elusive genes, we retrieved the genes that were retained by almost all of the mammalian species examined and defined them as ‘non-elusive,’ representing those persistent in the genomes. We conducted a careful search for gene loss to reduce the false discovery rate (FDR), which is usually caused by incomplete sequence information (Botero-Castro et al., 2017; Deutekom et al., 2019). By comparing the genomic regions containing these genes, we uncovered genomic characteristics relevant to gene loss. We associated the elusive genes with a variety of findings from deep sequencing analyses of the human genome, including transcriptomics, epigenomics, and genetic variations. These data assisted us to understand how intrinsic genomic features may affect gene fate, leading to gene loss by decreasing the expression level and eventually relaxing the functional importance of ‘elusive’ genes.

Results

Identification of human ‘elusive’ genes

We defined an ‘elusive’ gene as a human protein-coding gene that existed in the common mammalian ancestors but was lost independently in multiple mammalian lineages (Figure 1; see ‘Materials and methods’ for details). We searched for such genes by reconstructing phylogenetic trees of vertebrate orthologs and detecting gene loss events within the individual trees. To search for elusive genes, we paid close attention to distinguishing true evolutionary gene loss from falsely inferred gene loss caused by insufficient genome assembly, gene prediction, and orthologous clustering (Botero-Castro et al., 2017; Deutekom et al., 2019), as described below.

Detection of ‘elusive’ genes.

(a) Pipeline of ortholog group clustering and gene loss detection. (b) Definition of an elusive gene schematized with ortholog presence/absence pattern referring to a taxonomic hierarchy. Red and orange crosses denote the gene loss in the common ancestor of a taxon and the loss specific to a single species, respectively. (c) A representative phylogeny of the elusive gene encoding Chitinase 3-like 2 (CHI3L2). Taxa shown in the tree were used to investigate the presence or absence of orthologs. The Sciuromorpha, Hystricognathi, Eulipotyphla, Carnivora, and Chiroptera are absent from the tree, indicating that the CHI3L2 orthologs were lost somewhere along the branches framed in gray in the tree. In addition, the orthologs of many members of the Myomorpha were not found, suggesting that gene loss occurred in this lineage.

We first produced highly complete orthologous groups comprised of nearly complete gene sets. We merged multiple gene annotations of a single species followed by assessments of the completeness of the gene sets (Figure 1a). Using these gene sets, we then created two sets of ortholog groups with different methods and merged them into a single set (Figure 1a). In searching for gene loss events, we restricted our study to those that occurred in the common ancestors of particular taxonomic groups. This procedure relieved false identifications of gene loss in a species or an ancestor of a lower taxonomic hierarchy caused by incomplete genomic information (Figure 1b).

We integrated gene annotations from Ensembl, RefSeq, and the sequence repositories of individual genome sequencing projects to produce gene annotations for 114 mammalian and 132 non-mammalian vertebrates. From these, we selected the annotations of 101 and 90 species, respectively, that exhibited high completeness in the BUSCO assessment (Simão et al., 2015; Supplementary Table S1 in Supplementary file 1a). Using these gene sets, clustering of ortholog groups was conducted by OrthoFinder, and these groups were integrated into the ortholog groups provided by the Ensembl Gene Tree. This integration resulted in 50,768 vertebrate ortholog groups. Phylogenetic tree inference of the integrated ortholog groups and pruning of the individual trees based on gene duplications resulted in 17,495 mammalian ortholog groups that contained human genes. We classified the mammalian species into 15 taxonomic groups ranging from order to family (listed in Table S1; Supplementary file 1a). For the individual mammalian orthologs, we searched for the taxa in which the gene was absent in all the species examined (Figure 1b). We interpreted this gene absence as an evolutionary loss that occurred in the common ancestor of the taxon. Validating the gene loss through an ortholog search in genome assemblies and synteny-based ortholog annotations, we extracted the ortholog groups that were retained by humans but were lost independently in the common ancestors of at least two taxa (Figure 1c). Hereafter we call the human genes belonging to these ortholog groups ‘elusive genes.’ To compare these, we also selected the ortholog groups that contained all of the mammals examined including single-copy human genes. We called these ‘non-elusive genes.’ This comprehensive scan of gene phylogenies resulted in 813 elusive and 8050 non-elusive genes (Supplementary Table S2; Supplementary file 2).

Genomic signatures of the human elusive genes

The loss-prone nature of the elusive genes suggests a relaxation of their functional constraints. To uncover the molecular evolutionary characteristics associated with each elusive gene, we computed synonymous and non-synonymous substitution rates in coding regions, namely KS and KA, respectively, between human and chimpanzee and mouse orthologs for the elusive and non-elusive genes. In addition, we computed nucleotide substitution rates for introns (KI) between human and chimpanzee (Pan troglodytes) orthologs and compared them between the elusive and non-elusive genes. The results showed larger KA values in the ortholog pairs of the elusive genes than in those of the non-elusive genes (Figure 2a, Figure 2—figure supplement 1). This indicates a rapid accumulation of amino acid substitutions in the elusive genes, potentially accompanied by the relaxation of functional constraints. Our analysis further illuminated larger KS and KI values for the elusive genes than in the non-elusive genes (Figure 2b and c, Figure 2—figure supplement 1). Importantly, the higher rate of synonymous and intronic nucleotide substitutions, which may not affect changes in amino acid residues, indicates that the elusive genes are also susceptible to genomic characteristics independent of selective constraints on gene functions.

Figure 2 with 1 supplement see all
Genomic and evolutionary characteristics of elusive genes.

Distributions of non-synonymous, synonymous, and intronic nucleotide substitution rates, namely KA (a), KS (b), and KI (c) values, respectively, between the human–chimpanzee orthologs of the elusive and non-elusive genes. Distribution of gene length (d) and GC content (e) of the human elusive and non-elusive genes. (f) Distribution of gene density in the genomic regions where the human elusive and non-elusive genes are located. The plots consist of 249 elusive and 5145 non-elusive genes that retained chimpanzee orthologs (a, b), 473 and 4626 of those which harbored introns aligned with the chimpanzee genome (c; see ‘Materials and methods’), and all of the 813 elusive and 8050 non-elusive genes (d–f). Diamonds and bars within violin plots indicate the median and range from the 25th to 75th percentile, respectively.

To further scrutinize the characteristics reflecting the genomic environment rather than gene function, we analyzed genomic characteristics that may distinguish the elusive from non-elusive genes. A comparison between these two categories revealed shorter gene-body lengths and higher GC contents of elusive rather than non-elusive genes (Figure 2d and e). Furthermore, a scan of intragenomic gene distribution revealed that the elusive genes were located in the genomic regions with high gene density compared with the non-elusive genes (Figure 2f). Our findings indicate that such elusive genes have distinct characteristics in the human genome. These genomic characteristics, as well as high nucleotide substitution rates, were consistent with the findings in our genome analyses using the amniote and elasmobranch genomes (Hara et al., 2018a; Hara et al., 2018b).

Tracing elusiveness back along the vertebrate evolutionary tree

The origins of the human elusive genes can be traced back along the evolutionary tree, at least to the mammalian common ancestor. To investigate possible antiquities of the genomic properties associated with elusive genes, we investigated their orthologs in non-mammalian vertebrates by scrutinizing the ortholog groups used for elusive gene identification. We found that 152 out of 813 elusive genes originated in mammalian lineages, and this proportion was larger than those of the elusive genes (65 out of 8050, p=2.50 × 10-110), indicating that the elusive genes are more abundant in recently born genes than non-elusive genes. We then selected 517 elusive and 7900 non-elusive genes that originated in the common ancestors of jawed vertebrates or earlier. These subsets allowed us to examine the degree of retention of non-mammalian vertebrate orthologs in the elusive and non-elusive genes. On average, approximately 40% of these elusive genes were found to be retained by non-mammalian vertebrates, while this proportion increased up to 90% for the non-elusive genes. (Figure 3—figure supplement 1a). In the coelacanth, gar, and shark, the orthologs of the elusive genes were less frequently retained by all the species than those of the non-elusive ones (Figure 3—figure supplement 1b). The results suggest that the origins of the loss-prone propensity of the elusive genes potentially date back to the period long before the emergence of the Mammalia.

We further examined the genomic characteristics associated with the human elusive genes in the vertebrate orthologs. In all the species examined, orthologs of the elusive genes exhibited high GC content and compact gene bodies. Additionally, in most of these species, the orthologs of elusive genes were located in genomic regions with high gene density compared with orthologs of the non-elusive genes (Figure 3, Figure 3—figure supplement 2). In addition, we computed KS and KA values between the orthologs of the vertebrate species and their close relatives for elusive and non-elusive genes. In any of the species pairs except for avians, the orthologs of the elusive genes were found to harbor higher KA and KS values than those of the non-elusive gene orthologs (Figure 3, Figure 2—figure supplement 1). These observations indicate that these genomic characteristics probably originated before the emergence of gnathostomes, a monophyletic group of chondrichthyan and bony vertebrates, and have been retained for approximately 500 million years.

Figure 3 with 2 supplements see all
Long-standing characteristics of elusive genes.

Retention of the genomic and evolutionary characteristics of the human elusive genes across vertebrates. The individual round squares with arrows indicate significant increases or decreases of the distribution of particular characteristics in the orthologs of the human elusive genes and their flanking regions compared with those of the non-elusive genes in these selected vertebrate genomes. For the chimpanzee and mouse genomes, KA and KS values were computed between the human elusive genes and the orthologs of these mammals. For non-mammalian species, these values were computed with ortholog pairs for the elusive/non-elusive genes between the corresponding species and their closely related species: turkey for chicken, green anole for central bearded dragon, and whale shark for bamboo shark. Distributions of these metrics for non-human species are shown in Figure 2—figure supplement 1 and Figure 3—figure supplement 2. Species name: mouse, Mus musculus; chicken, Gallus gallus; central bearded dragon, Pogona vitticeps; Western clawed frog, Xenopus tropicalis; coelacanth, Latimeria chalumnae; spotted gar, Lepisosteus oculatus; bamboo shark, Chiloscyllium plagiosum.

Abundant polymorphism in elusive genes

The observation of large KS and KA values in the elusive genes prompted us to examine the extent to which these genes have accommodated genetic variations in modern humans. Large-scale human genome resequencing projects have identified a huge number of genetic variations, from rare to common, and from single-nucleotide variants (SNVs) to chromosome-scale structural variants, facilitating tackling this issue. We retrieved copy number variants (CNVs) and rare SNVs in the human genome from the Database of Genomic Variants, release 2016-08-31 (MacDonald et al., 2014) and dbSNP release 147 (Sherry et al., 2001), respectively, and computed their densities in the individual genic regions. We found that the genic regions of the human elusive genes contained abundant rare SNVs, as well as deletion and duplication CNVs, compared with those of the non-elusive genes (Figure 4a–c). This result suggests that genomic regions containing the elusive genes are not only prone to loss but also to duplication.

Genetic variations of the elusive and non-elusive genes within human populations.

Comparison of the density of rare single-nucleotide variants (SNVs) (a), deletion copy number variants (CNVs) (b), duplication CNVs (c), and Z-scores of synonymous (d), missense (e), and loss-of-function variants (f). We used opposite numbers of the Z-scores in d–f so that the elusive genes have higher values than non-elusive genes as in Figure 2a, b, c, e, f and Figure 3a–c. (a–c) 813 elusive genes and 8050 non-elusive genes were used. (d–f) 544 elusive genes and 7303 non-elusive genes for which genetic variants were available in GnomAD were used. Diamonds and bars within violin plots indicate the median and range from 25th to 75th percentile, respectively.

To evaluate the functional consequences of abundant genetic variants in the elusive genes, we investigated genetic variations stored in the gnomAD v. 2.1 database, a repository containing >120,000 exome and >15,000 whole-genome sequences of human individuals (Karczewski et al., 2021). This database classifies SNVs in coding regions into three categories—synonymous, missense, and loss-of-function—and the loss-of-function category contains nonsense mutations, frameshift mutations, and mutations in splicing junctions. The gnomAD site computes a Z-score, an index representing the abundance of SNVs for individual genes; positive and negative values denote fewer or more mutations in a coding region than expected, respectively (Figure 4d–f). Accordingly, the Z-score for nonsense mutations and loss-of-function mutations of the individual genes indicates the degree of natural selection: larger values demonstrate genes subjected to purifying selection, while smaller ones suggest functional relaxation. We found lower Z-scores of missense and loss-of-function mutations (higher opposite numbers of Z-scores in Figure 4e and f) in the human elusive genes than in the non-elusive genes, suggesting that the elusive genes are more functionally dispensable and potentially tolerable to harmful mutations. Additionally, opposite numbers of Z-scores of synonymous mutations of the human elusive genes were higher than those of the non-elusive genes (Figure 4d). This confirms the high mutability of genomic regions containing elusive genes, as observed in the KS values.

Transcriptomic natures of elusive genes

To further investigate how the human elusive genes have decreased functional essentiality, we examined their expression profiles. For this purpose, we compared gene expression profiles of the 54 adult tissues from the GTEx database v. 8 (The GTEx Consortium et al., 2020) between the elusive and non-elusive genes. For individual genes, we computed the maximum transcription per million (TPM) values among these tissues as the expression quantity level. For expression diversities, we employed Shannon’s diversity index H′, which is often utilized as an index of species diversity in the ecological literature, based on the proportion of TPM values across the 54 tissues.

As shown in the density scatter plots of the individual genes displaying these two indicators in Figure 5, most of the non-elusive genes possessed large maximum TPM and H′ values. Thus, most non-elusive genes are ubiquitously expressed at certain levels. By contrast, the density plot of the elusive genes displayed an additional high-density spot with small TPM and H′ values, indicating that the genes in this spot were not expressed, at least in adult tissues. The plot also showed another broad dense area of small H′ values, which contained the genes expressed in a single or a few tissues. A similar analysis was performed with the fetal single cell RNA-seq data (Cao et al., 2020), revealing that the averaged expression profiles of the elusive and non-elusive genes for the 172 cell types were concordant with those of the adult tissues (Figure 5). Our findings demonstrate that some elusive genes harbor low-level and spatially restricted expression profiles, that is, less pleiotropic states, which are rarely observed in the non-elusive genes.

Figure 5 with 1 supplement see all
Expression profiles of elusive and non-elusive genes.

The figure shows density scatter plots of the expression quantity and divergence of elusive and non-elusive genes. The numbers of the elusive/non-elusive genes and those for which the expression quantities were available are indicated in each panel. p-values were computed via 2 × 2 contingency tables presenting numbers of elusive and non-elusive genes with H′ < 1 and H′ ≥ 1. The median transcription per million (TPM) value of each of the adult tissue across individuals was retrieved from the GTEx database (The GTEx Consortium et al., 2020), and normalized TPM values of the fetal cell types were retrieved from the Descartes database (Cao et al., 2020). For the individual genes, maximum TPM and Shannon’s H′ values were computed using these processed TPM values.

Epigenetic nature of elusive genes

Our finding of the low-level and spatially restricted expression patterns of elusive genes prompted us to explore epigenetic properties involved in this transcriptional regulation. Therefore, we retrieved epigenetic data on a variety of human cell lines from a few regulatory genome databases including ENCODE, a repository that stores the comprehensive annotations of functional elements in the human genome (The ENCODE Project Consortium, 2012). Using this information, we characterized the epigenetic features of the genomic regions containing elusive genes (Figure 6).

Figure 6 with 5 supplements see all
Epigenetic features of the elusive genes.

Comparison of the distribution of ATAC-seq peak density (a), length of the topologically associating domains (TADs) including the elusive or non-elusive genes (b), the replication timing indicator based on Repli-seq (c), and overlap with the lamina-associated domains (LADs) computed from Lamin B1 ChIP-seq data. All of the analyses were performed by using the processed sequencing data publicly available (Table S3; Supplementary file 1b). (d) ATAC-seq and Hi-C were performed with A549 cells, Repli-seq was performed with HepG2 cells, and Lamin B1 ChIP-seq was performed with HAP-1 cells. In the elusive gene panels (orange bar), purple bar indicates the elusive genes with restricted expressions (H′ < 1; Figure 5). p-values for individual panels indicate the comparison between the elusive (813) and non-elusive (8050) genes and the one between the elusive genes with H′ < 1 (150) and those with H′ ≥ 1 (589). The results for other cells are shown in Figure 6—figure supplements 14 For the individual epigenetic characteristics, correction for multiple testing was performed for comparison in each cell cultures.

We compared peak densities based on the Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq), an indicator of accessible chromatin regions in the genome, in gene bodies and flanking regions between the elusive and non-elusive genes. In all of the eight cell lines examined (11 samples in total), the results showed fewer ATAC-seq peaks in the genomic regions including the elusive genes than in those including non-elusive genes, indicating that the elusive genes are likely to reside in inaccessible genomic regions (Figure 6a, Figure 6—figure supplement 1). We also searched for topologically associating domains (TADs), genomic elements with frequent physical self-interaction potentially acting as promoter-enhancer contacts (Rao et al., 2014) that included either the elusive or non-elusive genes. The result showed that a higher fraction of the elusive genes resided outside of the TADs than the non-elusive genes for all the eleven cell lines investigated (Figure 6b, Figure 6—figure supplement 2). Furthermore, the elusive genes were located in shorter TADs. These observations suggest that the elusive genes are unlikely to be regulated by distant regulatory elements compared with the non-elusive genes (Figure 6b, Figure 6—figure supplement 2).

Our investigations extended to the association of the elusive genes with further global regulation of genomic structures. We compared the percentage normalized signal of Repli-seq (Hansen et al., 2010), a high-throughput sequencing for quantifying DNA replication time as a function of genomic position, between the elusive and non-elusive genes. The results showed that elusive genes were prone to late replication in all of the 15 cell lines examined (Figure 6c, Figure 6—figure supplement 3). Late-replicating regions are frequently located at the nuclear periphery and often interact with the nuclear lamina. Therefore, we examined the nuclear position of the genomic regions including the elusive genes by referring to the lamina associating domains (LADs) that were identified by the ChIP-seq reads for Lamin B1 (van Schaik et al., 2020; Zheng et al., 2018). Compared with the non-elusive genes, the elusive genes were found to be enriched in LADs for all of the four cell lines examined (Figure 6d, Figure 6—figure supplement 4), consistent with their late replication timings (van Steensel and Belmont, 2017).

We further investigated the association of the restricted expressions of the elusive genes with epigenetic features. From 739 elusive genes whose expressions were quantified in the GTEx database, we classified the elusive genes into two groups based on the pleiotropy in terms of gene expressions: that is, 589 elusive genes with Shannon’s diversity index H′ ≥ 1 were ubiquitously expressed, that is, more pleiotropic, and 150 of those with H′ < 1 were expressed in only a few or none of the tissues examined, that is, less pleiotropic (Figure 5). Importantly, all of the four epigenetic features of the elusive genes with H′ < 1 were more pronounced than those with H′ ≥ 1: sparse ATAC-seq peaks, short TADs, late replication timings, and significant overlaps with LADs (Figure 6, Figure 6—figure supplements 14). This observation suggests that low-level and spatially restricted expressions of the elusive genes are associated with epigenetic features of these genomic regions.

High GC contents in genomic regions potentially hinder identifying an epigenetic feature by short-read sequencing because of the underrepresentation of sequence reads by amplification-based sequencing libraries. This bias might lead to sparse distributions of the ATAC-seq peaks and Hi-C contacts in the genomic regions that contain the elusive genes. However, only 3.00 and 9.00% of the elusive genes with H′ < 1 and H′ ≥ 1 were located in regions of extremely high GC content (>60%), respectively, showing that the elusive genes H′ ≥ 1 rather tend to contain more genes with high GC content (p=0.0176). Thus, the depleted epigenomic features in the genomic regions containing elusive genes are unlikely to be false discoveries caused by a technical issue, namely the underrepresentation of the sequencing reads.

Elusive gene orthologs in the chicken microchromosomes

The heterogeneous locations of the elusive genes can also be examined from a chromosome-scale viewpoint (Figure 7, Figure 7—figure supplement 1). The visualization via chromosome ideograms indicated an overlap of the elusive genes with the genomic regions enriched for the genes whose chicken orthologs are on the microchromosomes (chromosomes 11–38 and W), providing a statistical support for this trend (p=0.0175; Figure 7a). Indeed, microchromosomes of the chicken and other vertebrate exhibit genomic features including high GC content, high gene density, and rapid nucleotide substitutions in comparison with their macrochromosomes (Groenen et al., 2009; International Chicken Genome Sequencing Consortium, 2004; Schield et al., 2019; Waters et al., 2021), which also characterize genomic regions containing elusive genes. On the contrary, previous studies revealed that the chicken microchromosomes are preferentially located in the A compartments of the nucleus (Perry et al., 2020) and are early replicating (McQueen et al., 1998). These characteristics associated with the microchromosomes were opposite characteristics to the human genomic regions preferentially containing the elusive genes.

Figure 7 with 1 supplement see all
Chromosomal distribution of human elusive genes.

Red and dark blue horizontal bars beside the chromosome ideogram represent the location of elusive genes with restricted expression (Shannon’s H′ < 1) and the other elusive genes, respectively. (a) The chromosome diagrams are colored according to the density of the genes that harbor chicken orthologs in microchromosomes (number of genes/Mb). 93 and 68 elusive genes were orthologous to the chicken genes in macro- and microchromosomes, respectively, and 4211 and 2078 non-elusive genes were orthologous to the chicken genes in macro- and microchromosomes, respectively. This indicates that the chicken orthologs of the elusive genes are abundant in microchromosomes compared with those of the non-elusive genes (p=0.0175). (b) Gray regions in the diagram indicate orthologous regions of microchromosomes in the ancestors of gnathostomes (Nakatani et al., 2021). 395 and 296 elusive genes were located in the genomic regions corresponding to ancient macro- and microchromosomes, respectively, and 5950 and 1929 non-elusive genes were located in the genomic regions corresponding to these ancient chromosomes. The result recapitulates the biased localization of the elusive genes on microchromosomes (p=9.50 × 10-24). The chromosome diagrams were drawn using RIdeogram (Hao et al., 2020).

We further analyzed the ATAC-seq peaks in the chicken genome and found more peaks in the genomic regions including the elusive gene orthologs than in those containing non-elusive gene orthologs in four samples out of eight and no significant differences in the peak density in the four remaining samples (Figure 6—figure supplement 5). These observations indicate that, in an epigenetic manner, the chicken orthologs of the elusive genes are not regulated to reduce their expression level. This idea was further supported by a comparison of the expression profiles between the chicken orthologs of the elusive and non-elusive genes, showing no significant differences between them (Figure 5—figure supplement 1). Our analyses indicate that the genomic features of the elusive genes such as high GC and high nucleotide substitutions do not always correlate with a reduction in pleiotropy of gene expression that potentially leads to an increase in functional dispensability in the course of vertebrate evolution. In addition, avian orthologs of the elusive genes did not show higher KA and KS values than those of the non-elusive genes (Figure 3, Figure 2—figure supplement 1), likely consistent with not significant difference in gene expression levels between them in the species (Figure 5—figure supplement 1; Cherry, 2010; Zhang and Yang, 2015). We further compared expression profiles between the orthologs of the human elusive and non-elusive genes in several non-mammalian vertebrates and found that the orthologs of the elusive genes tend to exhibit low pleiotropy in green anole, coelacanth, and gar but not in Western clawed frog. The result suggests that the low pleiotropy of the elusive genes has persisted at least since the bony vertebrate ancestors. With respect to the chicken genome, the ‘elusive’ features for the genes orthologous to human elusive genes might have been relaxed—functional importance of the orthologs has increased—during evolution leading to chicken.

Discussion

Here we identified elusive genes that were lost in multiple lineages during mammalian evolution using a comprehensive scan of gene phylogenies. To identify gene loss events, absence of evidence (i.e. missing genes caused by incomplete genome assemblies and gene annotations) should be reviewed meticulously (Deutekom et al., 2019). Additionally, gene loss might be detected erroneously because of failure in similarity searches for orthologs of rapidly evolving genes (Moyers and Zhang, 2015). In this study, we aimed to reduce these false discoveries through our multifaceted approaches (Figure 1). We selected those species with highly complete gene annotations through integration of multiple gene annotations. Using these improved gene annotations, we created orthologous groups by employing a highly sensitive homology search with MMSeqs2 (Steinegger and Söding, 2017) and merged them into those identified in the Ensembl database. Furthermore, we restricted the loss events that were observed as gene absence in all species examined within all hierarchical levels of the selected taxonomic groups (Figure 1b). This absence is likely to have occurred as a gene loss in the common ancestor of the particular taxon rather than as a false discovery of gene loss in the individual species independently. Genuine continuous (e.g. telomere-to-telomere) genome assemblies are now available using modern sequencing technologies (Nurk et al., 2022). These genomic assemblies may help relieve the labor of examining for information losses, thereby facilitating the identification of genuine gene loss in any given species.

In the human genome, the elusive genes and their flanking regions harbor particular characteristics, including high GC content and high gene density, that may have originated long before the emergence of mammals (Figure 3). Frequent synonymous variations across modern humans in the elusive genes, consistent with higher synonymous substitution rates between the vertebrate orthologs, suggest that the genomic regions including elusive genes have been subject to rapid evolution for approximately 500 million years (Figures 2 and 4). Our findings indicate that heterogeneous genomic characteristics potentially affect the fate of genes at the latest period of vertebrate evolution. Analyses with large numbers of germline mutations in the human genome have illustrated the heterogeneity of mutation rates (Campbell and Eichler, 2013; Seplyarskiy and Sunyaev, 2021; Terekhanova et al., 2017). High GC content in the elusive genes may have facilitated an elevation of the mutation rate, as observed in the enrichment of rare variants in high-GC regions in the human genome (Schaibley et al., 2013). In addition, some of the elusive genes appear to have retained particular epigenetic marks including sparse ATAC-seq peaks, late replication timings, and location within LADs (Figure 6—figure supplements 14); these epigenetic marks are relevant to an increase in the mutation rate. Genomic regions with late replication timing exhibit increased mutation rates because of their unstable structure during the S-phase of the cell cycle (Koren et al., 2012; Stamatoyannopoulos et al., 2009). LADs retain more G-to-T mutations because of their susceptibility to oxidative damage in the nuclear periphery resulting in high levels of 8-oxoguanine (Yoshihara et al., 2014). Close coordination of the studies on gene evolution with germline mutation repertoires and spectra, which can be approximated from the collection of de novo mutations obtained by trio sequencing, may further facilitate our understanding of gene fates driven by heterogeneous genomic features—this would be viewed as ‘mutation-driven’ evolution (Nei, 2013).

The epigenetic marks of elusive genes are relevant to the suppression of gene expression (van Steensel and Belmont, 2017), and indeed, these genes harbor weakened and spatially restricted expression profiles (Figures 5 and 6 and Figure 6—figure supplements 14). However, the genomic features associated with these epigenetic marks usually exhibit lower GC contents and reduced gene density (Gilbert et al., 2004; Rao et al., 2014; van Steensel and Belmont, 2017). This discrepancy may be caused in part by a gain of local heterochromatin accompanied by suppression of the expression of transposable elements, as observed in various eukaryotic genomes (Choi and Lee, 2020; Fiston-Lavier et al., 2007; Grewal and Jia, 2007; Rangasamy, 2013; Slotkin and Martienssen, 2007; Underwood et al., 2017). Previous analyses showed frequent heterochromatinization of the human genomic regions where KRAB zinc finger genes colocalize with L1 retrotransposons (Imbeault et al., 2017; O’Geen et al., 2007; Vogel et al., 2006). One of the genomic regions found in human chromosome region 19p12 also contains many elusive genes (Vogel et al., 2006; Figure 7). Closer attention to the local gene and repeat contents including repetitive elements and tandem gene clusters might facilitate our understanding of heterochromatinization in restricted genomic regions, although we excluded such gene clusters in our search for elusive genes (Figure 1a).

A chromosomal-scale view of the distribution of elusive genes illuminated their significant correlation with the genes whose chicken orthologs are located on microchromosomes (Figure 7a). More importantly, genomic regions rich in elusive genes were traced back to the microchromosomes of the ancestral gnathostomes by reconstructing chromosomes of the ancestral genomes (Figure 7b). This inference of ancestral karyotypes augments our observations that some elusive natures of genomic sequences have been retained for hundreds of millions of years (Figure 3). In other words, the result suggests that the disparity of genomic regions that allows the ‘elusiveness’ for the genes has been retained during vertebrate evolution. On the other hand, comparisons of the expression profiles between the orthologs of the elusive and non-elusive genes for non-mammalian vertebrates suggest that the orthologs of the elusive genes have been associated with a reduction in pleiotropy of gene expression since vertebrate ancestors but acquired the diverse expressions in chicken and frog (Figure 5—figure supplement 1). Additionally, in the chicken genome, the diverse expressions of the chicken orthologs of the human elusive genes may be correlated with the abundance of ATAC-seq peaks (Figure 6—figure supplement 5). These findings again suggest that the chicken orthologs of the human elusive genes have increased pleiotropy of gene expression, which may lead to a lineage-specific acquisition of functional indispensability. It should be noted that the choices of tissues used in these analyses were largely different between the human and non-mammalian vertebrates (Tables S3 and S4; Supplementary file 1b and c). The chicken ATAC-seq data could be obtained only from developing embryos, while the human ATAC-seq in ENCODE were performed with cell lines. Therefore, the aforementioned interpretation should be treated carefully.

Finally, we note the potential evolutionary courses that facilitate the transition of gene fate from retention to loss. One possible course is a decrease in essential functions because of rapid sequence evolution in local genomic regions. The elusive genes located in those genomic regions with rapidly evolving characteristics are likely to accumulate neutral or even moderately harmful mutations in coding regions frequently, resulting in impaired essential functions. Another factor is the spatiotemporal suppression of gene expression via epigenetic constraints. Previous studies showed that lowly expressed genes are associated with low functional essentiality (Cherry, 2010; Gout et al., 2010), as shown for elusive genes in our study. Elusive genes with reduced pleiotropy may have limited opportunities to function, potentially leading to loss of their important roles. The extent of these evolutionary forces may have varied with time and lineages, resulting in a patchy loss of elusive genes phylogenetically. Interestingly, a recent large-scale scan of de novo mutations in Arabidopsis indicates the association of mutation rates with epigenetic features and functional essentiality of genes (Monroe et al., 2022). Further investigation of the association of genes with the surrounding genomic regions in various taxa may provide a common understanding of genomic and epigenomic features that potentially alter the fate of genes. Although epigenetic features are plastic, our findings indicate that the disparities of genomic regions are reflected in the heterogeneity of evolutionary forces and have been retained for hundreds of millions of years. This idea prompts us to explore evolutionary constraints on more global genomic regions that are potentially associated with structural characteristics including chromosomal composition and locations within the nucleus.

Materials and methods

Sequence retrieval

Request a detailed protocol

We retrieved genome assemblies and gene annotations of 114 mammals and 132 non-mammal vertebrates from RefSeq (accessed on April 9, 2018), Ensembl release 92, and the repositories of the individual genome projects (Supplementary Table S1 in Supplementary file 1a). Gene annotations for a single species from multiple repositories were integrated into one as follows. When gene annotations of multiple repositories were referring to the same version of the genome assembly, the annotation GTF files were merged with the ‘cuffcompare’ tool (Trapnell et al., 2012). Otherwise, translated amino acid sequences were clustered by CD-HIT v. 4.6.4 (Fu et al., 2012) with 100% sequence similarity, and the representative sequence for each cluster was retrieved by assuming that each cluster represented a single locus. Subsequently, we selected the canonical amino acid sequence for each locus: canonical peptides of the Ensembl genes were retrieved from the Ensembl database; for other resources, the longest amino acid sequence from the isoforms of a locus was chosen. The completeness of the gene annotations was performed on the gVolante web server with assessments by BUSCO v.2 (Simão et al., 2015) by referring to the vertebrate ortholog sets provided by BUSCO and CVG (Hara et al., 2015). The gene annotations of mammals, birds, and ray-finned fishes that had fewer than 1% missing genes, as well as those of the other vertebrates with fewer than 3% missing genes, were selected. Exceptionally, the gene annotations of Gavialis gangeticus (Reptilia; CVG missing ratio 3.86%), Paroedura picta (Reptilia; BUSCO vertebrate ortholog missing rate 3.25%), and Scyliorhinus torazame (Chondrichthyes; BUSCO vertebrate ortholog missing rate 4.45%) were added. Finally, the amino acid sequence set of 90 mammals and 101 non-mammalian vertebrates was subjected to t ortholog clustering. We also retrieved coding nucleotide sequences of the canonical amino acid sequences.

Ortholog clustering and tree inference

Request a detailed protocol

We retrieved gene trees of human protein-coding genes and their homologs from Ensembl Gene Tree release 92. From these gene trees, we constructed an amino acid sequence set of the homologs consisting of the species selected in the above section. This sequence set, restricted to Ensembl sequences only, was used as the ‘backbone’ of the ortholog set of all the selected species. In addition, we generated ortholog groups for all the species used by employing OrthoFinder2 v. 2.3.3 (Emms and Kelly, 2019) based on the similarity of amino acid sequences: a sequence similarity search was performed using MMSeqs2 v. 2339462c06eab0bee64e4fc0ebebf7707f6e53fd (Steinegger and Söding, 2017). The Ensembl and OrthoFinder ortholog sets were then merged to create the united set of ortholog groups, yielding 50,768 vertebrate ortholog groups.

The integrated ortholog groups were then subjected to molecular phylogenetic analysis. Amino acid sequences of the individual groups were aligned with MAFFT v. 7.402 (Katoh and Standley, 2013), and ambiguous alignment sites were removed with trimAl v1.4 (Capella-Gutiérrez et al., 2009). Phylogenetic trees were inferred with IQ-Tree v. 1.6.6 (Nguyen et al., 2015) by selecting the optimal amino acid substitution model with ModelFinder (Kalyaanamoorthy et al., 2017) implemented in the IQ-Tree tool for each sequence alignment. In the inferred phylogenetic trees, ambiguously bifurcated nodes—those with branch lengths less than 0.0025—were collapsed into a multifurcational node by the ‘di2multi’ function implemented in ape v. 5.5 (Paradis and Schliep, 2019). The trees were then rooted with the automatic rooting function ‘get_age_balanced_outgroup’ implemented in ete3 v. 3.1.1 (Huerta-Cepas et al., 2016) to minimize any discrepancy of tree topologies with the taxonomic hierarchy of the species included. Using the ortholog groups, the age of individual genes was estimated by inferring the oldest evolutionary lineage in the gene trees. We also adopted the gene age instructed by the Ensembl Gene Tree, wherever it shows an older age.

Identification of elusive genes in the human genome

Request a detailed protocol

For the individual trees, orthologs of the human genes were detected by the ‘get_my_evol_events’ function in ete3 (Huerta-Cepas et al., 2007). This function inferred gene duplication nodes in the rooted trees, resulting in separation of the trees into 17,495 subtrees of mammalian ortholog groups containing human genes. The ortholog information was referenced to extract the species with no orthologs to humans. This absence was further assessed by the ortholog annotation of human genes in the Ensembl Gene Tree database.

We selected taxonomic groups for the individual mammalian ortholog groups in which the orthologs were missing in all the species examined (Table S1; Supplementary file 1a). We restricted our study to gene losses that were likely to have occurred in the common ancestor of particular taxonomic groups, rather than those arising from the incompleteness of gene annotations. When a gene was missing in all the taxonomic groups in the same hierarchy, we considered that the gene was lost in the common ancestor of these groups. Finally, we found 1233 human genes belonging to the ortholog groups that were absent in two or more taxonomic groups and defined them as elusive genes. The gene loss events inferred by molecular phylogeny were further assessed by synteny-based ortholog annotations implemented in RefSeq, as well as a homolog search in the genome assemblies (Table S1; Supplementary file 1a) with TBLASTN v2.11.0+ (Altschul et al., 1997) and MMSeqs2 (Steinegger and Söding, 2017) referring to the latest RefSeq gene annotations (last accessed on December 2, 2022). This procedure resulted in the identification of 813 elusive genes that harbored three or fewer duplicates. Similarly, we extracted 8050 human genes whose orthologs were found in all the mammalian species examined and defined them as non-elusive genes. Because these elusive and non-elusive genes were identified in the GRCh38 human genome assembly, we performed the following analyses using this assembly.

Extraction of genomic and molecular evolutionary characteristics

Request a detailed protocol

We calculated the GC content of a gene by using its genomic region including introns and untranslated regions (UTRs). To calculate individual gene densities, we extracted genomic regions containing the genes and their flanking three genes at both ends and divided them by seven. The orthologs of the elusive and non-elusive genes were retrieved from the aforementioned gene trees. We computed KA and KS values of the ortholog pairs of human–chimpanzee, human–mouse, chicken–turkey (Meleagris gallopavo), central bearded dragon-green anole (Anolis carolinensis), and bamboo shark-whale shark (Rhincodon typus). To achieve this, we extracted ortholog groups that contained at least three of these ortholog pairs. Amino acid sequences of the human and the orthologs were aligned using MAFFT. Nucleotide sequence alignments of the coding regions were generated by ‘back-translation’ of the amino acid sequence alignments by trimAl, simultaneously removing ambiguous alignment sites. By employing coding nucleotide sequence alignments, numbers of synonymous and non-synonymous substitutions per site were computed using PAML v. 4.9a (Yang, 2007). To compute nucleotide sequence differences of the individual introns, we extracted 473 elusive and 4626 non-elusive genes that harbored introns aligned with the chimpanzee genome assembly. The nucleotide differences were calculated via the whole genome alignments of hg38 and panTro6 retrieved from the UCSC genome browser.

Multiomics analysis

Request a detailed protocol

Common and rare SNVs of the human populations were retrieved from dbSNP release 147 (Sherry et al., 2001), and human CNVs were obtained from the Database of Genomic Variants (DGV) release 2016-08-31 (MacDonald et al., 2014). The CNVs were classified into duplication and deletion variants, according to the annotation in DGV. The density of these variants in a gene was computed by dividing the number of variants identified in a gene region by its sequence length. Z-scores, indices of the tolerance against mutations, of synonymous, missense, and loss-of-function mutations of the individual genes were retrieved from gnomAD v. 2.1.1 (Karczewski et al., 2021).

Gene expression quantifications of adult and fetal tissues were retrieved from public databases. Expression profiles of adult tissues were obtained from the GTEx database v. 8 (The GTEx Consortium et al., 2020), computed by averaging TPM values across individuals. Expression profiles of fetal tissues were obtained from the Developmental Single Cell Atlas of gene Regulation and Expression (Descartes) portal (Cao et al., 2020) by calculating averaged TPM values of single cells. The maximum TPM values of the individual genes among the tissues were taken as the representative gene expression levels. As a proxy of the spatial diversity of gene expression, Shannon’s species diversity index (H′ values) was computed for each gene using the following equation:

Hi`=-k=1Rpi,klnpi,k

where Hi′ represents the Shannon’s index of ith gene in the list of the human genes, pi,k represents the proportion of the TPM values of the ith gene in the kth tissues/cell types, and R denotes the total number of tissues/cell types examined.

The ATAC-seq peaks and TAD boundaries of the human primary cells and culture strains were retrieved from the ENCODE 3 repository (Accession ID listed in Table S3; Supplementary file 1b; The ENCODE Project Consortium, 2012). Wavelet-smoothed signals of the ENCODE Repli-seq data were obtained from the UCSC genome browser (Hansen et al., 2010). The 20 kb bin-associated domains of LAD-seq that employed Lamin B1 antibodies (van Schaik et al., 2020) were retrieved from the 4D Nucleome Data Portal.

We also compared expression profiles and ATAC-seq peak densities between the orthologs of the elusive and non-elusive genes in non-mammalian vertebrates in a similar way as we did with the human datasets. Normalized gene expression profiles from RNA-seq data of normal adult tissues and early embryos for chicken, green anole, Western clawed frog, coelacanth, and spotted gar were obtained from the Bgee version 15 database (Bastian et al., 2021; Table S4; Supplementary file 1c). ATAC-seq narrow peak signals of chicken tissues were retrieved from NCBI GEO (Table S4; Supplementary file 1c) followed by coordination of the genome assembly with galGal5 with the UCSC liftOver tool (Hinrichs et al., 2006) as needed.

Code availability

Request a detailed protocol

The scripts for inferring gene presence and absence from gene trees were deposited in GitHub (https://github.com/yuichiroharajpn/ElusiveGenes, copy archived at Hara, 2022).

Statistical tests

Request a detailed protocol

Comparisons of the genomic characteristics between the elusive and non-elusive genes were tested statistically with the nonparametric Mann–Whitney U test and Fisher’s exact test implemented in R. Correction of multiple testing was performed using the Benjamini–Hochberg FDR approach. We considered p < 0.05 to be statistically significant.

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. Data from ENCODE Project was used and is available at https://www.encodeproject.org, IDs are shown in Table S3 in Supplementary File 1.

The following previously published data sets were used
    1. van Schaik T
    2. Vos M
    3. Peric-Hupkes D
    4. Hn Celie P
    5. van Steensel B
    (2020) 4D Nucleome Data Portal
    ID f1218a92-1f37-4519-85d6-ccedd5f7ad39. Dara from: Cell cycle dynamics of lamina-associated DNA.

References

  1. Book
    1. Nei M.
    (2013)
    Mutation-Driven Evolution
    Oxford University Press.
    1. The GTEx Consortium
    2. Aguet F
    3. Anand S
    4. Ardlie KG
    5. Gabriel S
    6. Getz GA
    7. Graubert A
    8. Hadley K
    9. Handsaker RE
    10. Huang KH
    11. Kashin S
    12. Li X
    13. MacArthur DG
    14. Meier SR
    15. Nedzel JL
    16. Nguyen DT
    17. Segrè AV
    18. Todres E
    19. Balliu B
    20. Barbeira AN
    21. Battle A
    22. Bonazzola R
    23. Brown A
    24. Brown CD
    25. Castel SE
    26. Conrad DF
    27. Cotter DJ
    28. Cox N
    29. Das S
    30. de Goede OM
    31. Dermitzakis ET
    32. Einson J
    33. Engelhardt BE
    34. Eskin E
    35. Eulalio TY
    36. Ferraro NM
    37. Flynn ED
    38. Fresard L
    39. Gamazon ER
    40. Garrido-Martín D
    41. Gay NR
    42. Gloudemans MJ
    43. Guigó R
    44. Hame AR
    45. He Y
    46. Hoffman PJ
    47. Hormozdiari F
    48. Hou L
    49. Im HK
    50. Jo B
    51. Kasela S
    52. Kellis M
    53. Kim-Hellmuth S
    54. Kwong A
    55. Lappalainen T
    56. Li X
    57. Liang Y
    58. Mangul S
    59. Mohammadi P
    60. Montgomery SB
    61. Muñoz-Aguirre M
    62. Nachun DC
    63. Nobel AB
    64. Oliva M
    65. Park Y
    66. Park Y
    67. Parsana P
    68. Rao AS
    69. Reverter F
    70. Rouhana JM
    71. Sabatti C
    72. Saha A
    73. Stephens M
    74. Stranger BE
    75. Strober BJ
    76. Teran NA
    77. Viñuela A
    78. Wang G
    79. Wen X
    80. Wright F
    81. Wucher V
    82. Zou Y
    83. Ferreira PG
    84. Li G
    85. Melé M
    86. Yeger-Lotem E
    87. Barcus ME
    88. Bradbury D
    89. Krubit T
    90. McLean JA
    91. Qi L
    92. Robinson K
    93. Roche NV
    94. Smith AM
    95. Sobin L
    96. Tabor DE
    97. Undale A
    98. Bridge J
    99. Brigham LE
    100. Foster BA
    101. Gillard BM
    102. Hasz R
    103. Hunter M
    104. Johns C
    105. Johnson M
    106. Karasik E
    107. Kopen G
    108. Leinweber WF
    109. McDonald A
    110. Moser MT
    111. Myer K
    112. Ramsey KD
    113. Roe B
    114. Shad S
    115. Thomas JA
    116. Walters G
    117. Washington M
    118. Wheeler J
    119. Jewell SD
    120. Rohrer DC
    121. Valley DR
    122. Davis DA
    123. Mash DC
    124. Branton PA
    125. Barker LK
    126. Gardiner HM
    127. Mosavel M
    128. Siminoff LA
    129. Flicek P
    130. Haeussler M
    131. Juettemann T
    132. Kent WJ
    133. Lee CM
    134. Powell CC
    135. Rosenbloom KR
    136. Ruffier M
    137. Sheppard D
    138. Taylor K
    139. Trevanion SJ
    140. Zerbino DR
    141. Abell NS
    142. Akey J
    143. Chen L
    144. Demanelis K
    145. Doherty JA
    146. Feinberg AP
    147. Hansen KD
    148. Hickey PF
    149. Jasmine F
    150. Jiang L
    151. Kaul R
    152. Kibriya MG
    153. Li JB
    154. Li Q
    155. Lin S
    156. Linder SE
    157. Pierce BL
    158. Rizzardi LF
    159. Skol AD
    160. Smith KS
    161. Snyder M
    162. Stamatoyannopoulos J
    163. Tang H
    164. Wang M
    165. Carithers LJ
    166. Guan P
    167. Koester SE
    168. Little AR
    169. Moore HM
    170. Nierras CR
    171. Rao AK
    172. Vaught JB
    173. Volpi S
    (2020) The gtex Consortium atlas of genetic regulatory effects across human tissues
    Science 369:1318–1330.
    https://doi.org/10.1126/science.aaz1776

Decision letter

  1. Wenfeng Qian
    Reviewing Editor; Chinese Academy of Sciences, China
  2. George H Perry
    Senior Editor; Pennsylvania State University, United States
  3. Wenfeng Qian
    Reviewer; Chinese Academy of Sciences, China

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Gene fate spectrum as a reflection of local genomic properties" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Wenfeng Qian as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by George Perry as the Senior Editor.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) The identification of "elusive genes" in the current manuscript requires additional scrutinization, given it is the foundation of the whole study. Please check published pipelines in the identification of gene losses (e.g., TOGA – https://github.com/hillerlab/TOGA) and use additional tools such as BLASTX search to test for known technical artifacts when calling genes (homology detection failure, refer to https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3000862). Please give a reason for the parameters used in the analysis (e.g., CD-Hit clustering) or examine if the conclusions remain supported by various parameters in the computational pipeline. Also, please take a look at the enrichment of elusive genes in human chr19, and use the synteny-based age estimation of the elusive genes (Shao et al. 2019).

2) Please also present the features of other genes (the genomic background other than elusive and non-elusive genes). Are these genes show intermediate patterns between those of elusive and non-elusive genes? Please edit the manuscript accordingly if the definition of non-elusiveness is actually equivalent to the genomic background.

3) It would be informative to test the links between recombination rate / LD and the genomic locations of elusive genes (compared against randomly sampled genes).

4) Please control confounding factors such as gene expression level and confirm whether the proxy of mutation rate (i.e., Ks) is actually confounded by gene importance.

5) Please consider extending the analyses on fish and birds to other genomic features.

6) The authors should reconcile the findings in this study with previous reports about microchromosomes.

7) Please consider improving the clarity/presentation of Figures 5 and 7 and examine whether the pattern remains robust using various parameter sets.

8) Please think of a better term than "elusive gene" to describe the genes that were lost independently in different lineages. Please also clearly define other terms in the manuscript (e.g., functionally indispensable vs. importance, are they the same concept?)

9) Please consider presenting the current study in the framework of mutation-selection balance, and better explain the novelty of the study over previous tremendous studies about gene losses.

Reviewer #1 (Recommendations for the authors):

Line 18. "neutral factor" is better replaced by "factors independent of gene dispensability".

Line 47. "However" should be "on the contrary"?

Lines 113-114. As indicated in the weaknesses part, I am not fully convinced these genomic, epigenomic, and transcriptomic features are completely independent of gene function.

Figure 1b. Define the red and orange crosses in the legend.

Figure 7 appeared first in the Discussion section. Can this part be moved to the Results section?

Reviewer #2 (Recommendations for the authors):

Overall, I believe that this is an interesting study. However, this version of the manuscript could be significantly improved in terms of logical depth and methodological stringency.

1. Authors actually support the concept of mutation-driven evolution, i.e., the high mutation rate in genomic regions harboring the elusive genes would predispose their fate toward death. To increase the significance of their work, I suggest authors cite (Nei 2013; Xie et al. 2019) and put their work in a bigger context.

2. Authors mentioned that elusive genes are less important and thus more prone to loss. In my view, pleiotropy is a better term compared to importance. That is, elusive genes are less pleiotropic [e.g. narrowly expressed, Figure 5] and thus their loss are more tolerable or easily compensated by other genes. Actually, narrow expression breadth has been observed to be correlated with gene loss in both humans and flies (MacArthur et al. 2012; Yang et al. 2015).

3. I am generally convinced that authors reliably identified elusive genes by identifying gene loss events in the common ancestor of multiple descendant species (to control for errors induced by assembly or annotation, Figure 1). However, Figure 7 shows the enrichment of elusive genes in human chr19. This chromosome is well known to be enriched with tandemly duplicated Krueppel-associated box C2H2 zinc-finger protein family (KZNF), many of which are primate-specific (Shao et al. 2019). I suspect that tree-based strategy implemented in Figure 1 could not be able to dissect the evolution of this super complex gene family. I am proposing two specific analyses: how many elusive genes encoded by chr19 are KZNFs? how many of them have Ensembl one-to-one orthologs across mammals?

4. With the patterns in Figure 3 and 7, authors argued that features of elusive genes are deeply ancient and could be inherited from the microchromosomes of early vertebrates. This statement has multiple problems.

a) Figure 3 only show genomic level features (e.g., high GC content) conserved in multiple vertebrates including shark and chicken. Epigenetic features analyzed in Figure 5 to 6 were only based on human data. I suggest authors to extend these analyses to shark or chicken. Although some epigenetic data could not be available for these species, transcriptome data analyzed in Figure 5 should be available for at least some species.

b) In Figure 7, authors propose the concept related with microchromosomes. These chromosomes have been extensively studied, especially in birds. Some features of microchromosomes are consistent with that of elusive genes [e.g., high GC, (Bravo et al. 2021)]. However, microchromosomes are conserved in terms of gene order and their genes generally show high protein-level constraints as shown by low Ka/Ks (Waters et al. 2021; Li et al. 2022). Authors need to reconcile their discovery with the previous rich literature.

c) Line (L) 204, among 982 human elusive genes, only 540~390 are shared by other species (e.g., shark). I suggest taking advantage of genome-level synteny based age data generated in Shao et al. 2019 to examine the age distribution of human elusive genes. If a high proportion of them are dated as being old (e.g., shared by jawed vertebrates), the statement that these genes have an ancient origin could be better supported.

References

Bravo GA, Schmitt CJ, Edwards SV. 2021. What have we learned from the first 500 avian genomes. Annu Rev Ecol Evol Syst 52: 611-639.

Li M, Sun C, Xu N, Bian P, Tian X, Wang X, Wang Y, Jia X, Heller R, Wang M et al. 2022. de novo Assembly of 20 Chicken Genomes Reveals the Undetectable Phenomenon for Thousands of Core Genes on Microchromosomes and Subtelomeric Regions. Mol Biol Evol 39.

MacArthur DG, Balasubramanian S, Frankish A, Huang N, Morris J, Walter K, Jostins L, Habegger L, Pickrell JK, Montgomery SB et al. 2012. A systematic survey of loss-of-function variants in human protein-coding genes. Science 335: 823-828.

Nei M. 2013. Mutation-driven evolution. OUP Oxford.

Shao Y, Chen C, Shen H, He BZ, Yu D, Jiang S, Zhao S, Gao Z, Zhu Z, Chen X et al. 2019. GenTree, an integrated resource for analyzing the evolution and function of primate-specific coding genes. Genome Res 29: 682-696.

Waters PD, Patel HR, Ruiz-Herrera A, Alvarez-Gonzalez L, Lister NC, Simakov O, Ezaz T, Kaur P, Frere C, Grutzner F et al. 2021. Microchromosomes are building blocks of bird, reptile, and mammal chromosomes. Proc Natl Acad Sci U S A 118.

Xie KT, Wang G, Thompson AC, Wucherpfennig JI, Reimchen TE, MacColl AD, Schluter D, Bell MA, Vasquez KM, Kingsley DM. 2019. DNA fragility in the parallel evolution of pelvic reduction in stickleback fish. Science 363: 81-84.

Yang H, He BZ, Ma H, Tsaur SC, Ma C, Wu Y, Ting CT, Zhang YE. 2015. Expression profile and gene age jointly shaped the genome-wide distribution of premature termination codons in a Drosophila melanogaster population. Mol Biol Evol 32: 216-228.

Reviewer #3 (Recommendations for the authors):

– In recent years several gene loss pipelines were already published (e.g. TOGA – https://github.com/hillerlab/TOGA) and it would be highly beneficial for this study to compare their gene loss reports with output obtained from existing pipelines (which also address the false discovery rate issue)

– We did not appreciate Figure 5. We strongly recommend finding a more quantitative approach to visualise these results, since heat maps are misleading and show different x-axis and y-axis ranges (zoom-ins/ zoom-outs)

– All Figures and Supplementary Figures showing violin plots need to report the number of genes that underly these distributions.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Gene fate spectrum as a reflection of local genomic properties" for further consideration by eLife. Your revised article has been evaluated by George Perry (Senior Editor) and a Reviewing Editor. A previous reviewer also read and commented on the revised manuscript.

Apparently, the manuscript has been significantly improved but there are some remaining issues that need to be addressed, as outlined below. In view of these comments, we kindly request that you consider revising the manuscript once more. Our hope is that, through this additional revision, the manuscript will be written more clearly and rigorously. Thank you for your understanding and continued efforts in improving your submission.

Reviewer #2 made the following comment. Please consider it during revision.

"The authors attempt to argue that the elusive status is ancient ("Thus, the heterogeneous genomic features driving gene fates toward loss have been in place since the ancestral vertebrates"). However, in response to my previous suggestion regarding chicken microchromosomes, the authors present mixed results. They observed high GC content, high gene density, and short gene length in chicken, similar to the findings in humans (Figure 3). Yet, the critical functional data between the two species are conflicting: human elusive genes exhibit low expression and fewer ATAC-seq peaks, while their chicken counterparts display the opposite pattern. In other words, chicken elusive genes exhibit higher pleiotropy, which may decrease the likelihood of their loss. Thus, these genes are not elusive, and high GC content, high gene density, and short gene length do not necessarily predict elusiveness. Given that the authors only analyzed the functional data of human and chicken genomes, it is not possible to determine whether the "elusive" status is ancient or derived from the human or mammalian lineage. I suggest that the authors analyze the transcriptome data in shark or spotted gar to provide further phylogenetic context. Otherwise, the authors should significantly tone down their statement."

The reviewing editor also has some comments on the title and abstract.

1. Please consider revising the title to "The Impact of Local Genomic Properties on the Evolutionary Fate of Genes Across Vertebrates".

2. Please consider adding a sentence (or something similar) at the end of the abstract. "This study sheds light on the complex interplay between gene function and local genomic properties in shaping gene evolution across vertebrate lineages."

https://doi.org/10.7554/eLife.82290.sa1

Author response

Essential revisions:

1) The identification of "elusive genes" in the current manuscript requires additional scrutinization, given it is the foundation of the whole study. Please check published pipelines in the identification of gene losses (e.g., TOGA – https://github.com/hillerlab/TOGA) and use additional tools such as BLASTX search to test for known technical artifacts when calling genes (homology detection failure, refer to https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3000862). Please give a reason for the parameters used in the analysis (e.g., CD-Hit clustering) or examine if the conclusions remain supported by various parameters in the computational pipeline. Also, please take a look at the enrichment of elusive genes in human chr19, and use the synteny-based age estimation of the elusive genes (Shao et al. 2019).

Based on this suggestion, we have repeated the analysis, and to more rigorously exclude possible false positives, and as a result, we have obtained a refined set of 813 elusive genes. The validation of gene loss has employed the ortholog annotations implemented in the RefSeq gene prediction pipeline, a synteny-based ortholog clustering platform other than TOGA, as well as similarity sequence search by MMSeqs2 that performs superior to NCBI BLAST.

For the revised manuscript, we have refined the elusive gene set as the reviewer suggested. In the genome assemblies, we have searched for the orthologs of the elusive genes for the species in which they were missing. The search has been conducted by querying amino acid sequences of the elusive genes with tblastn as well as MMSeqs2 that performed superior to tblastn in sensitivity and computational speed. In addition, regarding a comment by Reviewer 3 we have searched for the orthologs by referring to existing ortholog annotations. We used the ortholog annotations implemented in RefSeq instead of those from the TOGA pipeline: both employ synteny conservation. We have coordinated the identified orthologs with our gene loss criteria–absence from all the species used in a particular taxon–and excluded 268 genes from the original elusive gene set. These genes contain those missing in the previous gene annotations used in the original manuscript but present in the latest ones, as well as those falsely missing due to incorrect inference of gene trees. Finally, the refined set of 813 elusive genes were subject to comparisons with the non-elusive genes. Importantly, these comparisons retained the significantly different trends of the particular genomic, transcriptomic, and epigenomic features between them except for very few cases (Author response table 1). This indicates that both initial and revised sets of the elusive genes reflect the nature of the ‘elusiveness,’ though the initial set contained some noises. We have modified the numbers of elusive genes in the corresponding parts of the manuscript including figures and tables. Additionally, we have added the validation procedures in Methods.

Author response table 1
Difference in statistical significances across different elusive gene sets.
FeaturesNon-significant in theinitial gene set(1,081 elusive genes)Non-significant in thecurrent gene set(813)Non-significant in the current gene setexcluding chr19(669)
Gene density in the turkey genome
Gene density in the green anole genome
Gene density in the bamboo shark genome
Gene density in the whale shark genome
KS in avians
KA in avians
KA in sharks
ATAC-seq peak density for GM23338
Lamin B1 ChIP-seq peak density for K562
  1. The other features showed significantly different trends between the elusive and non-elusive genes for all of the elusive gene sets and thus are not included in this table.

cd-hit clustering with 100% sequence identity only clusters those with identical (and sometimes truncated) sequences, and, in the cluster, the sequences other than the representative are discarded. This means that the sequences remain if they are not identical to the other ones. If the similarity threshold is lowered, both identical and highly similar sequences are clustered with each other, and more sequences are discarded. Therefore, our approach that employs clustering with 100% similarity may minimize false positive gene loss.

Please refer to the reply to [Rev2 point 3] for the abundance of the elusive genes on chromosome 19. To examine this point, we have performed the comparisons between the elusive and non-elusive genes excluding the genes on chromosome 19, and the characteristics remained unchanged even when chromosome 19 was excluded.

2) Please also present the features of other genes (the genomic background other than elusive and non-elusive genes). Are these genes show intermediate patterns between those of elusive and non-elusive genes? Please edit the manuscript accordingly if the definition of non-elusiveness is actually equivalent to the genomic background.

Our aim in this study is to extract the characteristics of the genes that differentiate their fates from retention to loss. To achieve this, we compared the genes with clearly different phylogenetic signatures for gene loss, namely elusive and non-elusive genes.

The remainders excluding the elusive and non-elusive genes do not necessarily exhibit intermediate features. Because our definitions of the elusive and non-elusive genes were stringent, the remainders may contain considerable numbers of genes with loss in more restricted taxa than our criterion. In addition, the reminders contain those with other particular phylogenetic signals (e.g., frequently duplicated). These do not necessarily exhibit intermediate features, and at least the former is rather closer to elusive genes.

3) It would be informative to test the links between recombination rate / LD and the genomic locations of elusive genes (compared against randomly sampled genes).

We have retrieved fine-scale recombination rate data of males and females from https://www.decode.com/addendum/ (Suppl. Data of Kong, A et al., Nature, 467:1099–1103, 2010) and have compared them between the gene regions of the elusive and non-elusive genes. Both comparisons show no significant differences: average 0.829 and 0.900 recombinations/kb for the elusive and non-elusive genes, respectively, p=0.898, for males; average 0.836 and 0.846 recombinations/kb for the elusive and non-elusive genes, respectively, p=0.256, for females.

4) Please control confounding factors such as gene expression level and confirm whether the proxy of mutation rate (i.e., Ks) is actually confounded by gene importance.

We thank the reviewer for this important comment. We totally agree that transcriptomic and epigenomic features cannot be easily distinguished from gene dispensability and do not think that these features of the elusive genes can be explained solely by intrinsic properties of the genomes. Our motivation for investigating the expression profiles of the elusive gene is to understand how they lost their functional indispensability (original manuscript L285-286 in Results). We also discussed the possibility that sequence composition and genomic location of elusive genes may be associated with epigenetic features for expression depression, which may result in a decrease of functional constraints (original manuscript L470-474 in Discussion). Nevertheless, we think that the original manuscript may have contained misleading wordings, and thus we have edited them to better convey our view that gene expression and epigenomic features are related to gene function.

(P.2, Introduction) “This evolutionary fate of a gene can also be affected by factors independent of gene dispensability, including the mutability of genomic positions, but such features have not been examined well.”

(P6, Introduction) “These data assisted us to understand how intrinsic genomic features may affect gene fate, leading to gene loss by decreasing the expression level and eventually relaxing the functional importance of ʻelusiveʼ genes.”

(P33, Discussion) “Another factor is the spatiotemporal suppression of gene expression via epigenetic constraints. Previous studies showed that lowly expressed genes reduce their functional dispensability (Cherry, 2010; Gout et al., 2010), and so do the elusive genes.”

Additionally, responding to the advices from Reviewers 1 and 2, we have added a new section Elusive gene orthologs in the chicken microchromosomes in which we describe the relationship between the elusive genes and chicken microchromosomes. In this section, we also argue for the relationship between the genomic feature of the elusive genes and their transcriptomic and epigenomic characteristics. In the chicken genome, elusive genes did not show reduced pleiotropy of gene expression nor the epigenetic features relevant with the reduction, consistently with the moderation of nucleotide substitution rates. This also suggests that the relaxation of the ‘elusiveness’ is associated with the increase of functional indispensability.

(P27, Elusive gene orthologs in the chicken microchromosomes in Results) “Our analyses indicates that the genomic features of the elusive genes such as high GC and high nucleotide substitutions do not always correlate with a reduction in pleiotropy of gene expression that potentially leads to an increase in functional dispensability, although these features have been well conserved across vertebrates. In addition, the avian orthologs of the elusive genes did not show higher KA and KS values than those of the non-elusive genes (Figure 3; Figure 3—figure supplement 1), likely consistent with similar expression levels between them (Figure 5—figure supplement 1) (Cherry, 2010; Zhang and Yang, 2015). With respect to the chicken genome, the sequence features of the elusive genes themselves might have been relaxed during evolution.”

Also, please refer to [Rev1- point2] for the mutability of the elusive genes. To examine this point, we have computed nucleotide sequence differences in introns, namely KI, between the human and chimpanzee genomes. This analysis revealed higher KI values in the elusive genes than in the non-elusive genes, which is in line with our original hypothesis. The results have been added in the revised manuscript.

5) Please consider extending the analyses on fish and birds to other genomic features.

Please refer to the reply to [Rev2 point 4a] for the comparison between the elusive and non-elusive gene orthologs of nonmammalian vertebrates. We analyzed expression profiles and ATAC-seq peak densities of the elusive and non-elusive gene of these animals. We have created a new section Elusive gene orthologs in the chicken microchromosomes in Results and described the results in this section and Discussion.

We analyzed expression profiles and ATAC-seq peak densities of the elusive and non-elusive gene of these animals. We have created a new section Elusive gene orthologs in the chicken microchromosomes in Results and described the results in this section and Discussion.

We appreciate the reviewer for this meaningful suggestion. As a response, we have computed the differences in intron sequences between the human and chimpanzee genomes and compared them between the elusive and non-elusive genes. As expected, we found larger sequence differences in introns for the elusive genes than for the non-elusive genes. In Figure 2c of the revised manuscript, we have included the distribution of KI, sequence differences in introns between the human and chimpanzee genomes for the elusive and non-elusive genes. Additionally, we have added the corresponding texts to Results and the procedure to Methods as shown below.

(P11, Identification of human ‘elusive’ genes in Results) “In addition, we computed nucleotide substitution rates for introns (KI) between human and chimpanzee (Pan troglodytes) orthologs and compared them between the elusive and non-elusive genes.”

(P11, Identification of human ‘elusive’ genes in Results) “Our analysis further illuminated larger KS and KI values for the elusive genes than in the non-elusive genes (Figure 2b, c; Figure 2—figure supplement 1). Importantly, the higher rate of synonymous and intronic nucleotide substitutions, which may not affect changes in amino acid residues, indicates that the elusive genes are also susceptible to genomic characteristics independent of selective constraints on gene functions.”

(P39, Methods) “To compute nucleotide sequence differences of the individual introns, we extracted 473 elusive and 4,626 non-elusive genes that harbored introns aligned with the chimpanzee genome assembly. The nucleotide differences were calculated via the whole genome alignments of hg38 and panTro6 retrieved from the UCSC genome browser.”

6) The authors should reconcile the findings in this study with previous reports about microchromosomes.

Please refer to the reply to [Rev2-point 4b] for the reconciliation of our findings of the elusive genes with the features of microchromosomes. In the human genome, we found a significant overlap between the elusive genes and the genes whose chicken orthologs are located on microchromosomes. Although the chicken microchromosomes shared some sequence features such as high GC content and high KS values in common with the elusive genes in our sense but exhibited opposite transcriptomic and epigenetic trends. Although the result does not change the basis of our study in the human genome, it indicates that, in the course of evolution, genomic features of the elusive genes are not always associated with a reduction of pleiotropy of gene expression. The results have been described in the newly created section Elusive gene orthologs in the chicken microchromosomes in Results.

7) Please consider improving the clarity/presentation of Figures 5 and 7 and examine whether the pattern remains robust using various parameter sets.

Please refer to the replies to [Rev3- point 2] for this point. First, following this suggestion, we have conducted a statistical test to see whether the elusive genes contain more genes with restricted expression profiles (H’ < 1) than the non-elusive genes, and this trend was statistically supported. We have added the gene numbers of the individual categories and the result of this statistical tests to Figure 5. Therefore, we think the classification of the elusive genes based on the threshold (H’ = 1) is reasonable in Figure 7.

In addition, we have conducted statistical tests in a similar way with different thresholds for H’ (2, 3, and 0.5) and found that the pattern remains robust (p = 3.85x10-67, 1.29x10-50, and 9.40x10-46 setting thresholds H’ at 2, 3, and 0.5, respectively, for the GTEx dataset and p = 8.35x10-61, 2.01x1061, and 1.25x10-24 setting thresholds H’ at 2, 3, and 0.5, respectively, for the Descartes dataset).

To use Figure 7 in a new section in Results, we have added an ideogram showing the distribution of the genes that retain the chicken orthologs in microchromosomes. In response to the comment by Reviewer 2 [Rev2- point 4b], we have performed statistical tests and found that the elusive genes were significantly more abundant in orthologs in microchromosomes than the non-elusive genes. Furthermore, the observation that the elusive genes prefer to be located in gene-rich regions was already statistically supported (Figure 2f).

As shown in Figure 5, Shannon’s H' ranged from zero to approximately 4 (exact maximum value is 3.97) and 5 (5.11) for the GTEx and Descartes gene expression datasets, respectively. Although the threshold H'=1 was an arbitrarily set, we think that it is reasonable to classify the genes with high pleiotropy from those with low pleiotropy.

8) Please think of a better term than "elusive gene" to describe the genes that were lost independently in different lineages. Please also clearly define other terms in the manuscript (e.g., functionally indispensable vs. importance, are they the same concept?)

The phrase ‘elusive gene’ was already used in our previous paper to follow the advice of peer-reviewers for our past manuscripts, and for consistency, we would like use this term, while we have modified the sentence introducing this term to more carefully explain it.

9) Please consider presenting the current study in the framework of mutation-selection balance, and better explain the novelty of the study over previous tremendous studies about gene losses.

Please refer to the replies below to convey our novelty. We have cited existing literature in the revised manuscript.

Reviewer #2 (Recommendations for the authors): Overall, I believe that this is an interesting study. However, this version of the manuscript could be significantly improved in terms of logical depth and methodological stringency.

1. Authors actually support the concept of mutation-driven evolution, i.e., the high mutation rate in genomic regions harboring the elusive genes would predispose their fate toward death. To increase the significance of their work, I suggest authors cite (Nei 2013; Xie et al. 2019) and put their work in a bigger context.

We appreciate the reviewer for considering a way to enhance the significance of our study. We now recognize the importance of the study by Xie et al. (2019) reporting the recurrent loss of an enhancer that modulates fin morphogenesis in stickleback. As suggested, in the revised manuscript, we have cited this paper in Introduction.

(P4, Introduction) “In the stickleback genome, a Pitx1 enhancer was independently lost in multiple lineages inhabiting freshwater due to its genomic location in a structurally fragile site, leading to recurrent loss of pelvic fins (Xie et al., 2019). Genes and genomic elements in such particular regions may be prone to loss in a more neutral manner than the relaxation of functional importance or via functional adaptations.”

Additionally, to enhance the broad interests of our study, we have cited the Dr. Nei’s book in Discussion as shown below.

(P31, Discussion) “Close coordination of the studies on gene evolution with germline mutation repertoires and spectra, which can be approximated from the collection of de novo mutations obtained by trio sequencing, may further facilitate our understanding of gene fates driven by heterogeneous genomic features—this would be viewed as ‘mutation-driven’ evolution (Nei, 2013).”

2. Authors mentioned that elusive genes are less important and thus more prone to loss. In my view, pleiotropy is a better term compared to importance. That is, elusive genes are less pleiotropic [e.g. narrowly expressed, Figure 5] and thus their loss are more tolerable or easily compensated by other genes. Actually, narrow expression breadth has been observed to be correlated with gene loss in both humans and flies (MacArthur et al. 2012; Yang et al. 2015).

We thank the reviewer for the thoughtful suggestion. We agree that the word pleiotropy is more suitable in our manuscript. We have modified the manuscript as shown below.

(P20, Transcriptomic natures of elusive genes in Results) “Our findings demonstrate that some elusive genes harbor low-level and spatially-restricted expression profiles, i.e., less pleiotropic states, which are rarely observed in the non-elusive genes.”

(P23, Epigenetic nature of elusive genes in Results) “… we classified the elusive genes into two groups based on the pleiotropy in terms of gene expressions: that is, 589 elusive genes with Shannon’s diversity index H¢ ≥ 1 were ubiquitously expressed, i.e, more pleiotropic, and 150 of those with H¢ < 1 were expressed in only a few or none of the tissues examined, i.e., less pleiotropic (Figure 5).”

(P34, Discussion) “Elusive genes with reduced pleiotropy may have limited opportunities to function, potentially leading to loss of their important roles.”

3. I am generally convinced that authors reliably identified elusive genes by identifying gene loss events in the common ancestor of multiple descendant species (to control for errors induced by assembly or annotation, Figure 1). However, Figure 7 shows the enrichment of elusive genes in human chr19. This chromosome is well known to be enriched with tandemly duplicated Krueppel-associated box C2H2 zinc-finger protein family (KZNF), many of which are primate-specific (Shao et al. 2019). I suspect that tree-based strategy implemented in Figure 1 could not be able to dissect the evolution of this super complex gene family. I am proposing two specific analyses: how many elusive genes encoded by chr19 are KZNFs? how many of them have Ensembl one-to-one orthologs across mammals?

Responding to this comment, we investigated the KZNF genes on chromosome 19. We identified 206 those genes in chromosome 19, of which 75 were found to be elusive. Of these, 30 genes retained one-to-one orthologs of the mouse or dog. Although we excluded the nearly identical paralogs in our pipeline in the original manuscript and have refined the elusive gene set with synteny-based ortholog annotation in the current manuscript, the elusive KZNF genes have remained in this refined gene set. The elusive KZNF genes on chromosome 19 were lost in some mammalian lineages and duplicated in early primates.

Motivated by this comment, we conceived a possibility whether the enrichment of the paralogs of the elusive genes on chromosome 19 overrepresents the features relevant to the ‘elusiveness’. We thus have created another set of elusive genes, those excluding the genes on chromosome 19, and performed comparisons on the genomic, transcriptomic, and epigenomic features of the elusive and non-elusive genes. The results showed that the significant/non-significant differences have been maintained except for a few cases (see Author response table 1), indicating that the characteristics remained unchanged even when chromosome 19 is excluded.

4. With the patterns in Figure 3 and 7, authors argued that features of elusive genes are deeply ancient and could be inherited from the microchromosomes of early vertebrates. This statement has multiple problems.

a) Figure 3 only show genomic level features (e.g., high GC content) conserved in multiple vertebrates including shark and chicken. Epigenetic features analyzed in Figure 5 to 6 were only based on human data. I suggest authors to extend these analyses to shark or chicken. Although some epigenetic data could not be available for these species, transcriptome data analyzed in Figure 5 should be available for at least some species.

In response to this reviewer’s comment, we have investigated the transcriptomic and epigenetic characteristics of the orthologs of the elusive genes in non-mammalian vertebrates as we had done for human in the original manuscript. We have retrieved gene expression profiles of normal tissues and early embryos from bulk RNA-seq data of chicken, tropical clawed frog, coelacanth, and spotted gar, respectively, from the Bgee database (https://bgee.org/). We have then compared expression profiles between the orthologs of the elusive and non-elusive genes for the individual species using the same procedure as those in the manuscript. Only the anole, coelacanth, and gar orthologs of the elusive genes show the enrichment of low H’ values (newly created figure in Figure 5—figure supplement 1). This indicates that the low pleiotropy of expression of the elusive genes is not always observed in the non-mammalian species.

Furthermore, we have compared epigenomic properties between the orthologs of the elusive and non-elusive genes in the chicken genome. We have retrieved ATAC-seq narrow peak data for tissues of chicken embryos from NCBI GEO and compared the density of ATAC peaks between the orthologs of the elusive and non-elusive genes. The result indicates that, in five samples out of the eight, the orthologs of the elusive genes retained more ATAC peaks than those of the non-elusive genes, and that the reminders did not show this difference (newly created figure in Figure 6—figure supplement 5). This observation may remind us of a link between the reduction of the ‘elusiveness’ and the decrease of functional dispensability in gene evolution. However, it should be interpreted carefully, as different sets of tissues were used for the transcriptomic and epigenomic analyses between human and chicken. As described above, the chicken ATAC-seq experiments were mainly performed with developing embryos, while the human ATAC-seq used in our study were performed with cell lines. Nevertheless, the cross-species comparison of the epigenetic features suggests that the sequence features relevant to the elusive genes do not always induce the epigenetic conditions for gene expression depletion.

We have newly created Figure 5—figure supplement 1 and Figure 6—figure supplement 5 as shown below, described the results in the newly created section Elusive gene orthologs in the chicken microchromosomes in Results, and discussed them in Discussion. Also, we described the procedures in Methods.

(P-26-27, Elusive gene orthologs in the chicken microchromosomes in Results)

“Elusive gene orthologs in the chicken microchromosomes

The heterogeneous locations of the elusive genes can also be examined from a chromosome-scale viewpoint (Figure 7; Figure 7—figure supplement 1). The visualization via chromosome ideograms indicated an overlap of the elusive genes with the genomic regions enriched for the genes whose chicken orthologs are on the microchromosomes (chromosomes 11-38 and W), providing a statistical support for this trend (p=0.0175; Figure 7a). Indeed, microchromosomes of the chicken and other vertebrate exhibit genomic features including high GC-content, high gene density, and rapid nucleotide substitutions in comparison with their macrochromosomes (Groenen et al., 2009; International Chicken Genome Sequencing Consortium, 2004; Schield et al., 2019; Waters et al., 2021), which also characterize genomic regions containing elusive genes. On the contrary, previous studies revealed that the chicken microchromosomes are preferentially located in the A compartments of the nucleus (Perry et al., 2020) and are early replicating (McQueen et al., 1998). These characteristics associated with the microchromosomes were opposite characteristics to the human genomic regions preferentially containing the elusive genes.

We further analyzed the ATAC-seq peaks in the chicken genome and found more peaks in the genomic regions including the elusive gene orthologs than in those containing non-elusive gene orthologs in four samples out of eight and no significant differences in the peak density in the four remaining samples (Figure 6; Figure 6—figure supplement 5). These observations indicate that, in an epigenetic manner, the chicken orthologs of the elusive genes are not regulated to reduce their expression level. This idea was further supported by a comparison of the expression profiles between the chicken orthologs of the elusive and non-elusive genes, showing no significant differences between them (Figure 5; Figure 5—figure supplement 1). Our analyses indicate that the genomic features of the elusive genes such as high GC and high nucleotide substitutions do not always correlate with a reduction in pleiotropy of gene expression that potentially leads to an increase in functional dispensability in the course of vertebrate evolution. In addition, avian orthologs of the elusive genes did not show higher KA and KS values than those of the non-elusive genes (Figure 3; Figure 3—figure supplement 1), likely consistent with not significant difference in gene expression levels between them in the species (Figure 5—figure supplement 1) (Cherry, 2010; Zhang and Yang, 2015). With respect to the chicken genome, the sequence features of the elusive genes might have been relaxed during evolution.”

(P32-33, Discussion) “A chromosomal-scale view of the distribution of elusive genes illuminated their significant correlation with the genes whose chicken orthologs are located on microchromosomes (Figure 7a). More importantly, genomic regions rich in elusive genes were traced back to the microchromosomes of the ancestral gnathostomes by reconstructing chromosomes of the ancestral genomes (Figure 7b). This inference of ancestral karyotypes augments our observations that some elusive natures of genomic sequences have been retained for hundreds of millions of years (Figure 3). In other words, the result suggests that the disparity of genomic regions which allows the ‘elusiveness’ for the genes has been retained during vertebrate evolution. On the other hand, comparisons of the expression profiles between the orthologs of the elusive and non-elusive genes for non-mammalian vertebrates indicate that the elusive genes are not always associated with the restricted expression profiles (Figure 5; Figure 5—figure supplement 1). Additionally, in the chicken genome, this trend in gene expression may be correlated with the abundance of ATAC-seq peaks in the elusive genes (Figure 6—figure supplement 5). These findings again suggest that the elusive genes are not always associated with a reduction in pleiotropy of gene expression, which may lead to an increase of functional dispensability during evolution. It should be noted that the choices of tissues used in these analyses were largely different between the human and non-mammalian vertebrates (Tables S3, S4). The chicken ATAC-seq data could be obtained only from developing embryos, while the human ATAC-seq in ENCODE were performed with cell lines. Therefore, the aforementioned interpretation should be treated carefully.”

(P40, Methods) “We also compared expression profiles and ATAC-seq peak densities between the orthologs of the elusive and non-elusive genes in nonmammalian vertebrates in a similar way as we did with the human datasets.

Normalized gene expression profiles from RNA-seq data of normal adult tissues and early embryos for chicken, green anole, Western clawed frog, coelacanth, and spotted gar were obtained from the Bgee version 15 database (Bastian et al., 2021) (Table S4). ATAC-seq narrow peak signals of chicken tissues were retrieved from NCBI GEO (Table S4) followed by coordination of the genome assembly with galGal5 with the UCSC liftOver tool (Hinrichs et al., 2006) as needed.”

b) In Figure 7, authors propose the concept related with microchromosomes. These chromosomes have been extensively studied, especially in birds. Some features of microchromosomes are consistent with that of elusive genes [e.g., high GC, (Bravo et al. 2021)]. However, microchromosomes are conserved in terms of gene order and their genes generally show high protein-level constraints as shown by low Ka/Ks (Waters et al. 2021; Li et al. 2022). Authors need to reconcile their discovery with the previous rich literature.

Responding to this reviewer’s comment, we have first examined whether the chicken orthologs of the elusive genes tend to be located on microchromosomes. We have classified the chicken orthologs of the elusive and non-elusive genes into those on macro- and microchrosomes. Fisher’s exact test has indicated that elusive genes are enriched in the chicken microchromomes (p=0.0175, odds ratio=1.46; Author response table 2). In addition, the elusive genes are enriched in the genomic regions corresponding to the ancestral microchromosomes in early vertebrates (p=9.50x10-24, odds ratio=2.31; Author response table 3).

Author response table 2
Number of chicken orthologs of elusive and non-elusive genes locating in macro- and microchromosomes.
ElusiveNon-elusive
Macrochromosome (chr1-10, Z)934211
Microchromosome (chr11-, W)682078
  1. Genes in non-chromosome scaffolds and mitochondrial genome were excluded Fisher's exact test p=0.0175, odds ratio=1.48

Author response table 3
Number of elusive and non-elusive genes locating in the genomic regions derived from ancestral macro- and microchromosomes.
ElusiveNon-elusive
Ancestral macrochromosome3955950
Ancestral microchromosome2961929
  1. Genes in the genomic regions that did not correspond to the ancestral macro/microchromsomes were excluded. Fisher's exact test p=9.50x10-24, odds ratio=2.31.

As the reviewer noted, the genomic regions including the elusive genes share several genomic characteristics with microchromosomes: high GC content, high gene density, and high KA and KS values (although KA/KS values are lower for newly identified genes in the chicken microchromosomes in Li et al., 2022). On the contrary, previous analyses showed that the chicken microchromosomes tend to be early replicating and located in the A compartment in the nucleus, where the genes are actively transcribed. These observations are concordant with our finding that the chicken orthologs of the elusive genes are rich in ATAC peaks (see in the previous reply). We could not determine which epigenomic states are ancestral, but the genomic features of high GC-content and high nucleotide substitutions are not always correlated to reduced pleiotropy of gene expression, potentially leading to an increase in functional dispensability.

In the revised manuscript, we have created a new section Elusive gene orthologs in the chicken microchromosomes in Results to describe the enrichment of the elusive genes in the chicken microchromsomes and the association of their genomic characters with those of the chicken orthologs, while we cite the literature mentioned by the reviewer. In addition, we have mentioned the relevance of the elusive genes to the ancestral macrochromosomes in early vertebrates in Discussion. The corresponding sentences in the revised manuscript have already been included above in the previous.

c) Line (L) 204, among 982 human elusive genes, only 540~390 are shared by other species (e.g., shark). I suggest taking advantage of genome-level synteny based age data generated in Shao et al. 2019 to examine the age distribution of human elusive genes. If a high proportion of them are dated as being old (e.g., shared by jawed vertebrates), the statement that these genes have an ancient origin could be better supported.

We thank the reviewer for an opportunity to revisit the age of the elusive genes. We have examined this for the revised manuscript. Though the reviewer suggested using GenTree (http://gentree.ioz.ac.cn/) for this purpose, this database consists of euteleostome animals and does not include chondrichthyans and more distantly related animals. Instead, we have used the ortholog groups of jawed vertebrates built in this study and further added older branching dates of their orthologs and paralogs from the Ensembl Gene Tree. Of the 813 elusive genes in the revised dataset, 152 retain only mammalian orthologs, and this proportion is higher than that of the non-elusive genes (65 out of 8,050, p=2.50x10-110). We have then extracted the 517 elusive and 7,900 non-elusive genes whose ancestors were dated at the early evolution of jawed vertebrates or older. This subset of old gene age allowed us to examine if the elusive genes also retain loss-prone nature in non-mammalian vertebrates. On average, 40% of the 517 elusive genes are found to be retained in the nonmammalian vertebrates. On the other hand, more than 90% of the 7,900 non-elusive genes are retained by these species. These results indicate that while relatively young genes are frequently found in the elusive genes, the loss-prone nature of the elusive genes in the non-mammalian vertebrates is recapitulated by using this gene set of old origin. We have replaced Figure 3—figure supplement 1 with a revised version.

(P14, Tracing elusiveness back along the vertebrate evolutionary tree in Results) “We found that 152 out of 813 elusive genes originated in mammalian lineages, and this proportion was larger than those of the elusive genes (65 out of 8050, p = 2.50 ´ 10-110), indicating that the elusive genes are more abundant in recently born genes than non-elusive genes. We then selected 517 elusive and 7,900 non-elusive genes that originated in the common ancestors of jawed vertebrates or earlier. These subsets allowed us to examine the degree of retention of non-mammalian vertebrate orthologs in the elusive and non-elusive genes. On average, approximately 40% of these elusive genes were found to be retained by non-mammalian vertebrates, while this proportion increased up to 90% for the non-elusive genes. (Figure 3—figure supplement 1a). In the coelacanth, gar, and shark, the orthologs of the elusive genes were less frequently retained by all the species than those of the non-elusive ones (Figure 3—figure supplement 1b). The results suggest that the origins of the loss-prone propensity of the elusive genes potentially date back to the period long before the emergence of the Mammalia.”

Reviewer #3 (Recommendations for the authors):

– In recent years several gene loss pipelines were already published (e.g. TOGA – https://github.com/hillerlab/TOGA) and it would be highly beneficial for this study to compare their gene loss reports with output obtained from existing pipelines (which also address the false discovery rate issue)

Please refer to the reply to Essential revisions 1. We understand that TOGA is an ortholog search pipeline, whose output include gene loss information, employing synteny information. TOGA inputs genome alignments between the reference and all the target species, which requires tremendous computational time in our case involving a number of species. Instead, we have incorporated ortholog annotations implemented in RefSeq, another synteny-based ortholog detection platform, into our validation for gene loss.

– We did not appreciate Figure 5. We strongly recommend finding a more quantitative approach to visualise these results, since heat maps are misleading and show different x-axis and y-axis ranges (zoom-ins/ zoom-outs)

For the visibility of Figure 5, we have added the percentage of the genes of H'≥1 and H'<1 in the plots and performed the statistical test to see the fraction of the genes of H'≥1 and H'<1 is equal between the elusive and non-elusive genes. Additionally, we have coordinated x-axis and y-axis ranges between the elusive and non-elusive genes. We have replaced Figure 5 with a revised version containing these modified figure panels.

– All Figures and Supplementary Figures showing violin plots need to report the number of genes that underly these distributions.

We have added the sample numbers of the plots to all figures or figure legends.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Apparently, the manuscript has been significantly improved but there are some remaining issues that need to be addressed, as outlined below. In view of these comments, we kindly request that you consider revising the manuscript once more. Our hope is that, through this additional revision, the manuscript will be written more clearly and rigorously. Thank you for your understanding and continued efforts in improving your submission.

Reviewer #2 made the following comment. Please consider it during revision.

"The authors attempt to argue that the elusive status is ancient ("Thus, the heterogeneous genomic features driving gene fates toward loss have been in place since the ancestral vertebrates"). However, in response to my previous suggestion regarding chicken microchromosomes, the authors present mixed results. They observed high GC content, high gene density, and short gene length in chicken, similar to the findings in humans (Figure 3). Yet, the critical functional data between the two species are conflicting: human elusive genes exhibit low expression and fewer ATAC-seq peaks, while their chicken counterparts display the opposite pattern. In other words, chicken elusive genes exhibit higher pleiotropy, which may decrease the likelihood of their loss. Thus, these genes are not elusive, and high GC content, high gene density, and short gene length do not necessarily predict elusiveness. Given that the authors only analyzed the functional data of human and chicken genomes, it is not possible to determine whether the "elusive" status is ancient or derived from the human or mammalian lineage. I suggest that the authors analyze the transcriptome data in shark or spotted gar to provide further phylogenetic context. Otherwise, the authors should significantly tone down their statement."

We are aware of your concern about the inconsistency of the results between human and chicken, and therefore already included comparisons of the gene expression profiles between the orthologs of the elusive and non-elusive genes with several non-mammalian vertebrates in the previous revision (Figure 5—figure supplement 1). This comparison indicated that the orthologs of the elusive genes were rich in those with lower pleiotropy in anole, coelacanth, and gar, while this tendency was not found in chicken and frog. The result suggests that the orthologs of the elusive genes are likely to have retained the “elusive” features in the common ancestor of bony vertebrates.

In addition, regarding the reviewer’s comment, we should have clearly stated that in non-mammalians, the genes we focused on are not necessarily ‘elusive genes’ but rather ‘the orthologs of the human elusive genes’ in the previous manuscript. These orthologs do not always retain the ‘elusiveness’ as the human genes do. Interestingly, the orthologs of the elusive genes in the chicken genome do not exhibit significant differences in KA and KS values from those of the non-elusive genes, which is potentially associated with increased pleiotropy of gene expression in these genes. These observations suggest that the orthologs of the human elusive genes we identified have increased functional importance in the lineage leading to chicken.

We have revised the following parts of Results and Discussion to describe transcriptomic characteristics of the orthologs of the human elusive genes in non-mammalians more clearly in the current manuscript as included below.

Elusive gene orthologs in the chicken microchromosomes in Results

“We further compared expression profiles between the orthologs of the human elusive and non-elusive genes in several non-mammalian vertebrates and found that the orthologs of the elusive genes tend to exhibit low pleiotropy in green anole, coelacanth, and gar but not in Western clawed frog. The result suggests that the low pleiotropy of the elusive genes has persisted at least since the bony vertebrate ancestors. With respect to the chicken genome, the ‘elusive’ features the genes orthologous to human elusive genes might have been relaxed —functional importance of the orthologs have increased—during evolution leading to chicken.”

Discussion

“On the other hand, comparisons of the expression profiles between the orthologs of the elusive and non-elusive genes for non-mammalian vertebrates suggest that the orthologs of the elusive genes have been associated with a reduction in pleiotropy of gene expression since vertebrate ancestors but acquired the diverse expressions in chicken and frog (Figure 5; Figure 5—figure supplement 1). Additionally, in the chicken genome, the diverse expressions of the chicken orthologs of the human elusive genes may be correlated with the abundance of ATAC-seq peaks (Figure 6—figure supplement 5). These findings again suggest that the chicken orthologs of the human elusive genes have increased pleiotropy of gene expression, which may lead to a lineage-specific acquisition of functional indispensability.”

The reviewing editor also has some comments on the title and abstract.

1. Please consider revising the title to "The Impact of Local Genomic Properties on the Evolutionary Fate of Genes Across Vertebrates".

Thank you for your suggestion. We would like to use the title that you suggested without the last two words (‘The Impact of Local Genomic Properties on the Evolutionary Fate of Genes’). This title is more likely to increase the accessibility of our paper to a broad readership beyond those who are primarily investigating vertebrates, potentially prompting them to investigate genomic properties associated with the fate of genes in diverse organisms such as invertebrates, fungi, and plants.

2. Please consider adding a sentence (or something similar) at the end of the abstract. "This study sheds light on the complex interplay between gene function and local genomic properties in shaping gene evolution across vertebrate lineages."

Thank you for your suggestion that strengthens our focus and conclusion. We have revised the abstract as follows.

Thus, the heterogeneous genomic features driving gene fates toward loss have been in place and may sometimes have relaxed the functional indispensability of such genes. This study sheds light on the complex interplay between gene function and local genomic properties in shaping gene evolution that has persisted since the vertebrate ancestor.

https://doi.org/10.7554/eLife.82290.sa2

Article and author information

Author details

  1. Yuichiro Hara

    Research Center for Genome & Medical Sciences, Tokyo Metropolitan Institute of Medical Science, Tokyo, Japan
    Contribution
    Conceptualization, Data curation, Funding acquisition, Writing – original draft, Writing – review and editing, Investigation, Methodology
    For correspondence
    hara-yi@igakuken.or.jp
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7817-8963
  2. Shigehiro Kuraku

    1. Molecular Life History Laboratory, Department of Genomics and Evolutionary Biology, National Institute of Genetics, Mishima, Japan
    2. Department of Genetics, Sokendai (Graduate University for Advanced Studies), Mishima, Japan
    3. RIKEN Center for Biosystems Dynamics Research, Kobe, Japan
    Contribution
    Conceptualization, Supervision, Writing – review and editing, Funding acquisition
    Competing interests
    Reviewing editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1464-8388

Funding

Japan Society for the Promotion of Science (KAKENHI)

  • Yuichiro Hara

Japan Society for the Promotion of Science (KAKENHI)

  • Shigehiro Kuraku

Mochida Memorial Foundation for Medical and Pharmaceutical Research

  • Yuichiro Hara

Japan Society for the Promotion of Science (21K06132)

  • Yuichiro Hara

Japan Society for the Promotion of Science (20H03269)

  • Shigehiro Kuraku

Japan Agency for Medical Research and Development (JP21wm0325050)

  • Yuichiro Hara

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Yoichiro Nakatani for providing the information of orthologous regions of ancestral chromosomes in the human genome, and Hideya Kawaji and Ichiro Hiratani for insightful comments. We also would like to thank the peer reviewers for their constructive feedback and insightful comments, which have significantly improved this paper. This work was supported by RIKEN to SK, JSPS KAKENHI under Grant Number 20H03269 to SK and 21K06132 to YH, AMED under Grant Number JP21wm0325050 to YH, and Mochida Memorial Foundation for Medical and Pharmaceutical Research to YH. Computations were partially performed on the NIG supercomputer at the ROIS National Institute of Genetics.

Senior Editor

  1. George H Perry, Pennsylvania State University, United States

Reviewing Editor

  1. Wenfeng Qian, Chinese Academy of Sciences, China

Reviewer

  1. Wenfeng Qian, Chinese Academy of Sciences, China

Publication history

  1. Received: July 29, 2022
  2. Preprint posted: August 11, 2022 (view preprint)
  3. Accepted: April 25, 2023
  4. Version of Record published: May 24, 2023 (version 1)

Copyright

© 2023, Hara and Kuraku

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 603
    Page views
  • 70
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Yuichiro Hara
  2. Shigehiro Kuraku
(2023)
The impact of local genomic properties on the evolutionary fate of genes
eLife 12:e82290.
https://doi.org/10.7554/eLife.82290

Further reading

    1. Ecology
    2. Evolutionary Biology
    Jason P Dinh, SN Patek
    Research Article Updated

    Evolutionary theory suggests that individuals should express costly traits at a magnitude that optimizes the trait bearer’s cost-benefit difference. Trait expression varies across a species because costs and benefits vary among individuals. For example, if large individuals pay lower costs than small individuals, then larger individuals should reach optimal cost-benefit differences at greater trait magnitudes. Using the cavitation-shooting weapons found in the big claws of male and female snapping shrimp, we test whether size- and sex-dependent expenditures explain scaling and sex differences in weapon size. We found that males and females from three snapping shrimp species (Alpheus heterochaelis, Alpheus angulosus, and Alpheus estuariensis) show patterns consistent with tradeoffs between weapon and abdomen size. For male A. heterochaelis, the species for which we had the greatest statistical power, smaller individuals showed steeper tradeoffs. Our extensive dataset in A. heterochaelis also included data about pairing, breeding season, and egg clutch size. Therefore, we could test for reproductive tradeoffs and benefits in this species. Female A. heterochaelis exhibited tradeoffs between weapon size and egg count, average egg volume, and total egg mass volume. For average egg volume, smaller females exhibited steeper tradeoffs. Furthermore, in males but not females, large weapons were positively correlated with the probability of being paired and the relative size of their pair mates. In conclusion, we identified size-dependent tradeoffs that could underlie reliable scaling of costly traits. Furthermore, weapons are especially beneficial to males and burdensome to females, which could explain why males have larger weapons than females.

    1. Evolutionary Biology
    Anthony V Signore, Phillip R Morrison ... Kevin L Campbell
    Research Article

    The extinct Steller's sea cow (Hydrodamalis gigas; †1768) was a whale-sized marine mammal that manifested profound morphological specializations to exploit the harsh coastal climate of the North Pacific. Yet despite first-hand accounts of their biology, little is known regarding the physiological adjustments underlying their evolution to this environment. Here, the adult-expressed hemoglobin (Hb; a2β/δ2) of this sirenian is shown to harbor a fixed amino acid replacement at an otherwise invariant position (β/δ82Lys→Asn) that alters multiple aspects of Hb function. First, our functional characterization of recombinant sirenian Hb proteins demonstrate that the Hb-O2 affinity of this sub-Arctic species was less affected by temperature than those of living (sub)tropical sea cows. This phenotype presumably safeguarded O2 delivery to cool peripheral tissues and largely arises from a reduced intrinsic temperature sensitivity of the H. gigas protein. Additional experiments on H. gigas β/δ82Asn→Lys mutant Hb further reveal this exchange renders Steller's sea cow Hb unresponsive to the potent intraerythrocytic allosteric effector 2,3-diphosphoglycerate, a radical modification that is the first documented example of this phenotype among mammals. Notably, β/δ82Lys→Asn moreover underlies the secondary evolution of a reduced blood-O2 affinity phenotype that would have promoted heightened tissue and maternal/fetal O2 delivery. This conclusion is bolstered by analyses of two Steller's sea cow prenatal Hb proteins (Hb Gower I; z2e2 and HbF; a2g2) that suggest an exclusive embryonic stage expression pattern, and reveal uncommon replacements in H. gigas HbF (g38Thr→Ile and g101Glu→Asp) that increased Hb-O2 affinity relative to dugong HbF. Finally, the β/δ82Lys→Asn replacement of the adult/fetal protein is shown to increase protein solubility, which may have elevated red blood cell Hb content within both the adult and fetal circulations and contributed to meeting the elevated metabolic (thermoregulatory) requirements and fetal growth rates associated with this species cold adaptation.