Abstract
Nuclear receptors (NRs), a series of key transcription factors that are mostly activated by endogenous ligands or environmental xenobiotics, are reportedly good phylogenetic markers of animal genome evolution. As the early diverging class of bilaterians, however, a comprehensive view of the NR family in a marine free-living flatworm Macrostomum lignano and comparative information in flatworms are still lacking, which is of significance to address the evolutionary diversification of the NR family and imply the adaptive evolution in the early diverging Bilateria. Herein, a total of 51, 26, and 23 putative NR genes were identified in M. lignano, Sparganum proliferum, and Clonorchis sinensis, respectively, which were classified into eight subfamilies, implying an extensive expansion of the NR family in M. lignano. It is presumed that the extensive expansion was mainly attributed to the M. lignano-specific hidden polyploidy, segmental, and tandem duplication events. The duplicated NR pairs in M. lignano and the NR orthologs in flatworms all experienced the purifying selection. Phylogenetic analyses indicated the presence of NR3-like genes in M. lignano, which is first reported in flatworms. Intron loss and reduced intron size were mainly contributed to the structural divergence of NR genes in flatworms. The combined data provide indispensable information for a better understanding of the complexity and the adaptive evolution of the NR gene family in metazoans.
Introduction
Nuclear receptor (NR) gene superfamily comprises a large group of ligand-regulated transcription factors which are involved in various functions such as reproduction, differentiation, development, homeostasis, metabolism, and metamorphosis (; ; ). The corresponding ligands of NRs contain a variety of endogenous molecules (e.g., estrogen, androgen, and secondary bile acids) and environmental xenobiotics (e.g., pharmaceutical agents, synthetic hormones, biocides, plastics, and personal care products) (Weatherman et al., 1999; ). Consequently, NRs act as a conduit between the internal and external environments, closely linking to the function of endocrine systems (). Although NRs are responsible for a tremendous diversity of functions, they share a common structural modular. A canonical NR possesses five to six functional regions including the N-terminal A/B domain, a highly conserved C domain (DNA binding domain, DBD), D domain (hinge), a moderately conserved E domain (ligand-binding domain, LBD), and F domain (; ). Characterized by the most conserved domains DBD and LBD, the typical NRs are clustered into seven subfamilies, NR1–NR6 and NR8 (; ). However, atypical NRs are grouped into NR0 with either DBD (NR0A) or LBD (NR0B) and NR7 with two DBD and an LBD in some animal lineages (Wu et al., 2007; ). Additionally, all ctenophore NRs were atypical without a DBD domain ().
Nuclear receptors are considered strong phylogenetic markers of animal genome evolution due to their common ancestral origin, distribution in all metazoan genomes, and highly conserved throughout the whole animal taxa (). Therefore, the genome-wide identification of NR genes can provide a path to link evolutionary and functional genomics, and trace the dynamic evolutionary route in metazoans. Recently, the rapid developments in the genome sequencing have facilitated a genome-wide identification of NR family members in many species. The NR gene repertoires have been characterized in different lineages with a clear species-dependent pattern from invertebrates to vertebrates, such as 2 members in sponge, 4 in placozoan, ≥250 in nematodes, ∼21 in arthropods, ∼48 in mammals, and 66–137 in teleost (; ; ; Yang et al., 2020). These data suggest the evolutionary complexity of the NR family, which is likely driven by gene duplication and loss (). The early diverging bilaterians offer an exceptional opportunity to address the evolution and diversification of the NR family in bilaterians.
Flatworms occupy an assumed central position in the evolution of the Bilateria for that they have been widely reported as the early diverging bilaterians (; ). Moreover, the flatworms have a characteristic of a simple bilateral body plan with the development of complex endocrine systems, which is crucial to coordinate the reaction of organisms to the environment (). They are conventionally divided into the free-living turbellarians and parasitic tapeworms and flake, and consequently exposed to different environments, including terrestrial, freshwater, and marine environments for the former and within the host tissues for the latter (). Some parasitic flatworms (such as Schistosoma mansoni and Echinococcus multilocularis) and free-living fresh turbellarians (Schmidtea mediterranea) have been surveyed for NR homologs (Wu et al., 2006; Tharp et al., 2014; Wu and LoVerde, 2019); however, our knowledge regarding the repertoires of NRs in marine groups and the early branching flatworm lineages remains to be filled.
Macrostomum lignano from the high-tide interstitial sand fauna of the coast of the Mediterranean Sea is a marine free-living regenerative flatworm, which is an affiliate to the Macrostomorpha that is an earlier diverging Bilateria than the other often-studied free-living and parasitic flatworms (Wasik et al., 2015). Thus, M. lignano is an exceptional model organism for addressing the diversification of the NR family and implying the evolution of the endocrine system in the early diverging Bilateria. Moreover, a comparative study of NR genes in M. lignano and the other representative flatworms may provide important clues for understanding the adaptative evolution of early bilaterians. Herein, the comprehensive analyses based on the genome-wide identification, genomic distribution, and gene expansions of the NR family were performed in M. lignano. Furthermore, systematic investigations were also conducted to uncover the phylogenetic relationship, exon-intron organization, and selection pressures of this gene family in the representative flatworms.
Materials and Methods
Identification and Retrieval of NR Genes in the Selected Flatworms
All possible NR homologs were identified in the representative flatworms, including M. lignano (PRJNA371498), S. proliferum (PRJEB35374), and C. sinensis (PRJNA386618), by a four-step genome search strategy: (1) the reported NR homologs of parasitic flatworms were searched in the predicted proteomes of the representative flatworms by the BLASTP program with the default settings and similarity thresholds of ≥40% in the NCBI web1; (2) all NR sequences were aligned using the Bioedit software to retrieve the conserved domain DBD and LBD sequences; (3) to avoid omissions in the repertoires of NR genes, the DBD and LBD domain sequences were used as a query to search against the target genome sequences by iterative TBLASTN (E = 2e–5) until no novel sequence was retrieved; and (4) all of the positive hits were then identified and confirmed for the presence of NR genes by comparing them to the NCBI database using BLASTP, and redundancy was removed by sequence alignments and genomic location. Allowing for the incorrect genome assembly caused by the high repeated sequences content in the M. lignano genome, we artificially removed the NR genes with sequences that are the same in the middle of a nucleotide sequence but slightly different at the ends to minimize the number of potential redundant genes. To avoid the interference caused by the allelic variants, an analysis testing on the nucleotide sequence similarity starting from 90% with an increment of 1% was performed to choose an optimal cut-off value. Afterward, a cut-off value of 95% was chosen to eliminate the redundancy. Eventually, 51 potential NR proteins were identified and renamed based on their phylogenetic relationships.
Phylogenetic Analysis of NR Genes
Multiple sequence alignments of NR proteins from M. lignano, Amphimedon queenslandica, Trichoplax adhaeren, Nematostella vectensi, Crassostrea gigas, Bactrocera dorsalis, Branchiostoma floridae, and Homo sapiens were performed using the Bioedit software to obtain the aligned DBD and LBD amino acid sequences. Based on the alignments, Maximum-likelihood (ML) phylogenetic trees were built by the IQTREE v1.6.8 () with the automatic selection of an optimal model for protein substitution and rate heterogeneity. For branch support analysis, the SH-aLRT test and ultrafast bootstrapping were conducted with 1,000 replicates. Bayesian phylogenetic tree was constructed by MrBayes3.2 () with the following parameters: generations = 2,000,000, number of runs = 2, burnin fraction = 0.25, and temp = 0.15. The ML and Bayesian trees were visualized and polished by the iTOL v6 website2.
Genomic Distribution and Gene Structure Analyses of NR Genes in M. lignano
Chromosome size and genomic distribution of NR genes in M. lignano were obtained from the NCBI3 and displayed by the TBtools software with the Amazing Gene Location plugin. The information of NR gene structure was from the genomic annotations downloaded from the NCBI “(see text footnote 3).” The gene structure of NRs combined with the phylogenetic relationship were displayed by the TBtools software with the Amazing Optional Gene Viewer plugin ().
Investigations on the Gene Duplication Events of NRs in M. lignano
Duplicated NR genes were analyzed by e-value cutoff BLAST-searching (e-value < 10–10) against each other, and identified as duplicated genes when the nucleic acid sequence identity of two or more NR genes reaches at least 80% and aligned region between the two genes were covered >80% of the longer genes (). The NR genes that satisfied these two conditions are considered duplicates. To further explore the potential expansion mechanisms of the NR gene family in M. lignano, we performed all-vs.-all comparisons by the local BLASTP program to identify synteny blocks. SD events were identified according to two conditions as follows: (1) the NR genes derived from gene duplication were located on different scaffolds but clustered into the same group; and (2) the duplicated NRs were found in two or more syntenic regions which were defined to be ≥10 kb in size and ≥80% in sequence identity (). These regions were identified as SDs. The genes adjacent to the NR genes were analyzed to determine whether tandem duplication had occurred. Tandem duplication events were characterized when a pair of duplicated NRs is nearby or separated by few genes in a 100-kb region of a chromosome (Wang et al., 2010).
Selection Pressure Analysis
The identification of positive selection and purifying selection at each site was performed by the SELECTON Server website4 with the ratio of non-synonymous (Ka) to synonymous substitutions (Ks), namely the Ka/Ks ratio (). If the Ka/Ks ratio of the sites is significantly greater than 1, it means that these sites is under positive pressure; that is, these sites in a sequence varied rapidly during evolution. If the Ka/Ks ratio of the sites is less than 1, it is interpreted as an evidence for negative selection. If the Ka/Ks ratio is 1, the sites is under natural selection, that means the absence of natural selection (). Based on these, the CDS sequences in each group of flatworms NR genes were submitted together to the SELECTON Server. The evolutionary model is M8, and other parameters are set as the default values.
Results
Genome-Wide Identification and Phylogenetic Analysis of NR Genes in M. lignano
To identify NR members in M. lignano, the NR proteins of parasitic flatworms and the consensus DBD or LBD domain sequences were used as a query to search against the M. lignano genome database. A total of 140 candidate NR genes with at least one DBD or LBD were originally obtained from the M. lignano genome. By removing 21 redundant sequences, 119 NR genes were obtained based on the presence of an apparently complete DBD or LBD domain. Since the remarkably complex genome structure with ∼75% of the sequence containing repeats and transposon sequences (Wasik et al., 2015), it is difficult to confirm the exact NR gene repertoires in M. lignano. To get more accurate repertoires of NR genes, an analysis testing on the nucleotide sequence similarity starting from 90% with an increment of 1% was performed to choose an optimal cut-off value. The total numbers of NR genes in M. lignano were 90, 62, 55, 52, and 51 with a cut off value of 99, 98, 97, 96, and 95%, respectively, and remained unchanged until to 90% (Supplementary Table 1 and Extended Table). Given the incorrect genome assembly caused by the high repeated sequences content and allelic variants in M. lignano, a cut-off value of 95% was chosen to filter the potential redundancy. Eventually, 51 NR genes were identified and renamed according to the phylogenetic relationships. The identified NR genes ranged from 654 to 2,958 bp, and the corresponding proteins ranged from 217 to 985 amino acids in length. The characteristic information of these NRs was listed in Supplementary Table 2, including the accession number, CDS length, and chromosomal locations.
For a comprehensive understanding of the phylogenetic relationships, an unrooted phylogenetic tree containing 198 NR protein sequences from Amphimedon queenslandica, Trichoplax adhaeren, Nematostella vectensis, Crassostrea gigas, Bactrocera dorsalis, Branchiostoma floridae, Homo sapiens, and M. lignano was constructed using the DBD plus LBD amino acid sequences by the IQTREE procedure with the ML method (Figure 1). All NR genes were clearly divided into nine major groups, as previously reported, referred to as subfamilies NR0–NR8. M. lignano possesses NR members belonging to eight of the nine NR subfamilies. There are two atypical NRs (NR0A1 and NR0A2) containing only one DBD from the NR0 subfamily, which is grouped with NR0A proteins from B. dorsalis rather than C. gigas, B. floridae, and H. sapiens. The largest family is the NR1 subfamily with 27 members, including 1 NR1As, 4 NR1Cs, 2 NR1Fs, and 20 NR1Js. The NR2 subfamily is made up of 7 genes. Interestingly, we have identified 10 members in the NR3 subfamily. Additionally, the rest contained 1 NR4A, 1 NR5B, 2 NR7 (2DBD), and 1 NR8 members. We did not identify any NR6 gene in the M. lignano genome.
FIGURE 1
Identification of NR3 Members Present in M. lignano
The above phylogenetic tree suggested that 10 NR3 proteins from M. lignano all clustered into the NR3B group with the other species. To further confirm the classification of these 10 NR genes, ML and Bayesian phylogenetic trees of the NR3 subfamily were constructed using software IQTREE1.6 and MrBayes 3.2, respectively, based on the multiple alignments of DBD plus LBD sequences (Figure 2). Comparative results showed that the topological structures of the three phylogenic trees were generally similar. However, the NR3 members from M. lignano (MlNR3s) were grouped with the vertebrate NR3C genes in the Bayesian phylogenetic trees, which is different from the result of the ML tree. Normally, the Bayesian Inference method shows a higher accuracy than the other tree-building methods () and, therefore, results were considered more credible. Nevertheless, the bootstrap value is moderate (0.532).
FIGURE 2
The analyses on the sequence identity of DBD and LBD in NR3s between M. lignano and the other species were also performed (Supplementary Tables 3, 4). The results suggested that the DBDs of MlNR3s have 39.3–55%, 44.9–64.4%, and 37.2–55% shared identity compared to the NR3A (ERs), NR3Bs (ERRs), and NR3Cs of the other species, respectively. Similarly, the LBDs of MlNR3s showed a low shared identity with 8.3–31.5%, 13.7–37.8%, and 9.5–28.7% when compared to the NR3A (ERs), NR3Bs (ERRs), and NR3Cs of the other species, respectively. This indirectly supported the result of the ML phylogenetic tree. The DBD and LBD sequences of MlNR3s were compared with the ERs, ERRs, and NR3Cs from the other species (Supplementary Figures 1, 2). The DBDs showed a high level of conservation for all sequences, and were more similar to the ERRs than the NR3Cs. For example, the remarkable conservation of four residues (21E, A22, 30T, and 32Q) in the ERRs was also observed in some of the MlNR3s, while two residues (78P and 79A) displayed conservation in Amphioxus SR, human AR, GR, MR, PR, and MlNR3B8-15. The LBD that constitutes the ligand-binding pocket (LBP) was less conserved. Similar to the other NR3 proteins, the LBDs of the MlNR3s also contained 12 α-helices (H1–H12). The detailed comparison of the LBD sequences suggested that the conserved residues involving ligand binding in the LBP of the ERs and ERRs were also found in the H3, H4–5, H7, H10–11, and H12 of the MlNR3s (Supplementary Figure 2), whereas only one conserved residue of human NR3Cs for oxosteroid binding was observed in the MlNR3s.
Genomic Distribution and Gene Duplication Events of NR Genes in M. lignano
Phylogenetic results in this study revealed a massive expansion of the NR gene family in the M. lignano genome, especially in the NR1 and NR3 subfamilies. Generally, gene duplication events, including whole genome duplication (WGD), segmental duplication (SD), and tandem duplication (TD), has been recognized as the main mechanisms to drive the evolution and expansion of gene families. To explore the reason of NRs expansion in M. lignano, the genomic location was mapped. A complete view of the NR genes localization in the M. lignano genome clearly showed that NR genes displayed a highly dispersed distribution of 51 NR members across the 50 scaffolds (Figure 3 and Supplementary Table 5). Only two NR3-like genes located on Sc-73 formed a cluster, namely NR3-like12 and NR3-like13. As tandem duplicated genes represent an array of at least two homologous genes with a high sequence similarity of the encoded proteins separated by a chromosomal region within 100-kb, these 2 NR genes were clearly clustered into one tandem repeat event region. Therefore, they were identified as TD genes. Additionally, the genomic regions containing the NR members which are potential synteny were searched to locate the segmentally duplicated pairs. The searches at the whole-genome scale confirmed the widespread occurrence of SDs. Only one syntenic blocks containing 2 NR genes were detected in the M. lignano genome (Supplementary Figure 3), indicating that they were involved in SD. However, no TD or SD events were found in the NR gene family of the other flatworms.
FIGURE 3

Chromosomal locations and gene duplication events of NR genes in M. lignano. Their detailed scaffold IDs are recorded in Supplementary Table 5. The NR genes from SD are linked by blue dashed lines, and from TD are marked by green star. The detailed map of the TD and SD is provided in Supplementary Figure S3.
Identification of NR Genes in Other Representative Flatworms
Given that the comparative analyses of gene family can provide insights into interpreting the biological variation among flatworms in the light of evolution, we also identified the NR genes in other representative flatworms with 26 NR members in Sparganum proliferum (tapeworm) and 23 in C. sinensis (fluke) using the same methods (Supplementary Table 6). By combining the previously reported data of S. mansoni and S. mediterranea NR genes, comparison analyses on the scaffold number, genome size, coverage, total gene number, and identified NR gene number of six representative flatworms, including two free-living species, two tapeworms, and two flukes, were conducted (Table 1). The results revealed that the number of identified NR genes varies with the genome size of parasitic flatworm species in the range of 17--26 with similar genome sequencing quality and bioinformatic methods, but the free-living species. Among these species, the M. lignano genome contains a maximum of 51 NR genes, while E. granulosus has 17 NRs. The number of M. lignano NRs is 0.96--2 times more than that of the other flatworms, while the genome size was approximate. By searching the Genome database of NCBI5, we found that the total gene number of M. lignano is 49,027, which approximately is 1, 3.3, 1.2, 2.3, and 3.2 times larger than that of S. mediterranea (23,657), E. granulosus (11,325), S. proliferum (22,739), C. sinensis (14,936), and S. mansoni (11,713), respectively.
TABLE 1
| Species name | Genome size (Mb) | Coverage | Total gene number | Identified NR genes | References |
| M. lignano | 764.41 | 99.79% | 49,027 | 51 | This study |
| S. mediterranea | 773.94 | 99.99% | 23,657 | 23 | Wu et al., 2006 |
| E. granulosus | 110.84 | 99.34% | 11,325 | 17 | |
| S. proliferum | 653.39 | 92.19% | 22,739 | 26 | This study |
| C. sinensis | 562.77 | 96.84% | 14,936 | 23 | This study |
| S. mansoni | 409.58 | 97.72% | 11,713 | 21 | Wu et al., 2019 |
Genomic information and identified NR gene numbers in representative species of flatworms.
The data were obtained from NCBI (https://www.ncbi.nlm.nih.gov).
The numbers of NR genes in certain groups were distinct among these species (Figure 4). We found that the number of genes in NR1 and NR3 were larger as compared to the other subfamilies in M. lignano. For instance, the NR1 subfamily was composed of 27 genes from M. lignano, but only 5 genes from S. mansoni, and 6 from S. proliferum, C. sinensis, and S. mediterranea. Additionally, 10 NR genes from M. lignano were classified under the NR3 subfamily, while no member in the NR3 subfamily was observed in the other flatworms. These results suggested that these groups have undergone an extensive expansion in M. lignano. Additionally, the largest clade is NR1, followed by NR3, NR2, NR0, NR7, NR5, and NR8 in M. lignano, while the NR2 subfamily is the largest group followed by the NR1, NR7, NR5, NR0, and NR4 subfamilies in the other flatworms. As compared to the other flatworms, NR1B and NR2F homologs were absent in the M. lignano genome. Conversely, the NR8 subfamily is unique to M. lignano. In addition, no homolog to the NR6 subfamily was identified in the flatworms.
FIGURE 4

Different distribution patterns of NR genes in flatworms. The x-axis represents the species, and the y-axis represents the gene numbers in different subfamilies.
Comparison of NR Gene Structure in Flatworms
The structural stability of a gene is a prerequisite to maintain their biological functions, while divergence in the structures of a gene is essential for studying evolutionary divergence within members of gene families. Therefore, the evolutionary dynamics of an exon-intron organization are indicators for the evolutionary history of a gene family (
FIGURE 5

Phylogenetic relationships and exon-intron organizations of NR genes in M. lignano, S. proliferum, and S. mansoni. (A) Phylogenetic relationships. The phylogenic tree was constructed using full length NR protein sequences from the three species. The eight major groups were marked with different color backgrounds. (B) Exon-intron organizations of NR genes. Box represent exons, of which, green boxes signify coding sequence (CDS), and yellow boxes represent untranslated region introns (UTR). Black lines represent introns.
FIGURE 6

Differences of introns in M. lignano, S. proliferum, and S. mansoni. (A) Distribution of intron numbers in flatworm species. (B) Distribution of intron length in flatworm species.
Selective Pressure Analysis of NR Genes
Performing selection pressure analysis at the genetic level helps to understand not only the evolutionary history of organisms, but also the structural and functional variations of genes. According to the phylogenetic analysis, we found that the linage-specific gene duplication events happened in the NR family in M. lignano. For the sake of investigations on the selection pattern of these NR paralogs after duplication events, which might also contribute to the overall divergence of the NR gene family, the non-synonymous substitutions (Ka), synonymous substitutions (Ks), and their ratio Ka/Ks values of two pairs of paralogous NRs were obtained, respectively. As shown in Supplementary Table 8, the Ks value of NR genes derived from TD (1.7553) is far greater than that from SD (0.0449). The Ka/Ks ratios of the duplicated NR pairs in M. lignano were all lower than 1, implying that the duplicated NR genes has evolved under the purifying selection. To further assess which selection patterns drove the evolution of the NR gene family in flatworms, we also calculated the Ka/Ks ratio of orthologous NR genes among flatworms (Supplementary Table 9). The Ka/Ks ratios for all groups of NR genes in flatworms were <1, ranging from 0.2215 to 0.6785 with an average of 0.4405. These observations suggested that the evolution of the NR gene family in flatworms is under an intensive purifying selection pressure.
Discussion
The in-depth investigation of well-chosen gene families appears as the most promising path for linking evolution to functional adaptations (
Among the extant animals, the flatworms are close to one of the earlier groups of bilaterians, occupying the important evolutionary position in the tree of life. Importantly, the marine free-living M. lignano is the early diverging sub-taxon of the Platyhelminthes–Rhabditophora (
Whole-Genome Identification and Phylogenetic Analysis of NR Genes in Flatworms
A total of 119 NR genes were observed in M. lignano. After removing the potential redundancy caused by incorrected assembly and allelic variants, 51 NRs were eventually identified with a cut-off value of 95% sequence identity. Most of the removed NR genes with a high sequence similarity (>97%) were located in the duplicated segments as reported by Wasik et al. (2015). It has been reported that the segmentally duplicated genes, especially those that occurred recently, have a high sequence similarity with a range of 90–99.5% (Samonte and Eichler, 2002). Even the segmental duplication has been defined as “long stretches of duplicated sequences that can span between 1–200 kb and that share a sequence identity higher than 95%” by
Based on the phylogenetic grouping relationships, a total of 198 NR proteins were organized into nine major groups, from NR0–NR8 in selected species (Figure 1), which largely coincided with the previous reports (
The NR gene repertoire is significantly larger than those in S. mediterranea (21), E. granulosus (17), S. proliferum (26), C. sinensis (23), and S. mansoni (21) with similar genome sequencing quality and bioinformatic methods (Table 1 and Supplementary Table 4), clearly showing the size multiplicity of this gene family among different flatworms. The genome size varied in different flatworms, such as M. lignano (764.41 Mb), E. granulosus (110.84 Mb), S. proliferum (653.39 Mb), C. sinensis (562.77 Mb), and S. mansoni (409.58 Mb) (Table 1), suggesting that the number of NR family members has an absolute correlation with the genome size in parasitic flatworms. However, the genome size of free-living species S. mediterranea (773.94 Mb) is much larger than that of the other parasitic flatworms, while it bears less NR genes, implying that the different evolutionary histories for each species is another reason except for the genome size. This was assumedly attributed to the relatively similar environment surrounding the parasitic flatworms, which is one of the evolutionary driving forces that contribute to the final shape of the organism genome (
Comparative analyses suggested that 12 of 17 subfamilies were shared by M. lignano and the other flatworms; this means that 12 subfamilies of M. lignano contained NR orthologs to the other flatworms (Supplementary Table 6). Remarkably, the NR3 and NR8 subgroups were exclusive to M. lignano. As the second-oldest subfamily, a member of the NR3 subfamily (ERR) has been identified in the placozoan T. adhaerens (
Extensive Expansion of NR Family in M. lignano
The number of NR genes identified in M. lignano was much more than the other flatworms, which is nearly three times, showing an extensive expansion in the M. lignano genome. The independent species-duplications were also detected in the homeobox (Hox) gene family in M. lignano, leaving multiple copies of Hox genes (Wasik et al., 2015). Generally, gene duplication was considered a major mechanism for gene expansion, including whole WGD, SD, and TD (
The extensive expansion was mainly attributed to the impressive number of NR1J, NR3, and NR1C genes with 20, 10, and 4 duplicates in M. lignano, respectively. Seemingly, species- or lineage-specific gene duplications appeared to be relatively frequent in invertebrates. For instance, lineage-specific gene duplication of 10 NR1H has been observed in Branchiostoma floridae (
Evolutionary Conservation and Divergence of NR Genes in Flatworms
The structure of a specific gene regarding the intron numbers and length is typically conserved in the closely related species (
During the evolutionary process, the conservation or divergence must occur in the NR gene family of flatworms, especially in M. lignano, for extensive duplication. In general, duplicated genes reportedly undergo different evolutionary fates after duplication events, including functional conservation, neo-functionalization, sub-functionalization, and pseudogenization (
Conclusion
The present study displays a genome-wide comprehensive identification and comparative analyses of NR family in three flatworms and provides insight for understanding the mode of evolution of the NR family in flatworms. Our analyses suggest an extensive expansion of the NR gene family with 51 members in M. lignano as compared to the other flatworm species (17–26 members). Remarkably, the NR3-like genes were identified in M. lignano, which is not presented in the other flatworms. Based on the genomic location and syntenic analysis, the expansion of NR family seems to be attributed to the M. lignano-specific hidden tetraploidy, SDs and TD. Furthermore, the evolutionary conservation and diversity of NR genes in flatworms is denoted by the comparative analyses on the exon-intron organizations and selective pressure. Intron loss and a reduction in intron size have been observed in M. lignano, which may mainly attribute to the massive gene duplications. The purifying selection is the primary force driving the evolution of the NR gene family in flatworms. This study provides important cues for addressing the evolutionary diversification of NR family, and understanding the adaptative evolution in the early diverging Bilateria.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
YYC contributed to the conceptualization, analyses, interpretation, writing – original draft, and writing – review and editing. JLC contributed to the software and data analyses. IM contributed to the writing – review and editing. JMC contributed to the project administration and funding acquisition. All authors have read and agreed to the submitted version of the manuscript.
Funding
This research was funded by China Ocean Mineral Resources Research and Development Association Program (Grant No. DY135-B2-13), Program for Innovative Research Team in Science and Technology of Education Department of Fujian Province, Project of Fujian Provincial Department of Science and Technology (Grant No. 2020J01868), Project of Minjiang University (Grant No. MYK20004), and National Natural Science Foundation of China (Grant No. 31902352).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.653447/full#supplementary-material
References
1
AlviteG.RieraX.CancelaS.PaulinoM.EstevesA. (2019). Bioinformatic analysis of a novel Echinococcus granulosus nuclear receptor with two DNA binding domains.PLoS One14:e0224703. 10.1371/journal.pone.0224703
2
AuwerxJ.BaulieuE.BeatoM.Becker-AndreM.BurbachP. H.CamerinoG.et al (1999). A unified nomenclature system for the nuclear receptor superfamily.Cell97161–163. 10.1016/s0092-8674(00)80726-6
3
BainD. L.HeneghanA. F.Connaghan-JonesK. D.MiuraM. T. (2007). Nuclear receptor structure: implications for function.Annu. Rev. Physiol.69201–220. 10.1146/annurev.physiol.69.031905.160308
4
BakerM. E. (2019). Steroid receptors and vertebrate evolution.Mol. Cell Endocrinol.496:110526. 10.1016/j.mce.2019.110526
5
BertrandS.BrunetF. G.EscrivaH.ParmentierG.LaudetV.Robinson-RechaviM. (2004). Evolutionary genomics of nuclear receptors: from twenty-five ancestral genes to derived endocrine systems.Mol. Biol. Evol.21, 1923–1937. 10.1093/molbev/msh200
6
BodofskyS.KoitzF.WightmanB. (2017). Conserved and exapted functions of nuclear receptors in animal development.Nucl. Receptor. Res.4:101305. 10.11131/2017/101305
7
BridghamJ. T.EickG. N.LarrouxC.DeshpandeK.HarmsM. J.GauthierM. E. A.et al (2010). Protein evolution by molecular tinkering: diversification of the nuclear receptor superfamily from a ligand-dependent ancestor.PLoS Biol.8:e1000497. 10.1371/journal.pbio.1000497
8
CaoJ.HuangJ.YangY.HuX. (2011). Analyses of the oligopeptide transporter gene family in poplar and grape.BMC Genomics12:465. 10.1186/1471-2164-12-465
9
CasesI.de LorenzoV.OuzounisC. A. (2003). Transcription regulation and environmental adaptation in bacteria.Trends Microbiol.11248–253. 10.1016/S0966-842X(03)00103-3
10
ChenC.ChenH.ZhangY.ThomasH. R.FrankM. H.HeY.et al (2020). TBtools: an integrative toolkit developed for interactive analyses of big biological data.Mol. Plant.13, 1194–1202. 10.1016/j.molp.2020.06.009
11
ChengX.ZhangD.ChengZ.KellerB.LingH. Q. (2009). A new family of ty1-copia-like retrotransposons originated in the tomato genome by a recent horizontal transfer event.Genetics181:1183. 10.1534/genetics.108.099150
12
ChengY. Y.TaoW. J.ChenJ. L.SunL. N.ZhouL. Y.SongQ.et al (2015). Genome-wide identification, evolution and expression analysis of nuclear receptor superfamily in Nile tilapia, Oreochromis niloticus.Gene569141–152. 10.1016/j.gene.2015.05.057
13
CollinsJ. J. (2017). Platyhelminthes.Curr. Biol.27252–256. 10.1016/j.cub.2017.02.016
14
CruzeiroC.Lopes-MarquesM.RuivoR.Rodrigues-OliveiraN.SantosM. M.RochaM. J.et al (2016). A mollusk VDR/PXR/CAR-like (NR1J) nuclear receptor provides insight into ancient detoxification mechanisms.Aquat. Toxicol.17461–69. 10.1016/j.aquatox.2016.02.007
15
Doron-FaigenboimA.SternA.MayroseI.BacharachE.PupkoT. (2005). Selecton: a server for detecting evolutionary forces at a single amino-acid site.Bioinformatics21, 2101–2103. 10.1093/bioinformatics/bti259
16
EggerB.LaprazF.TomiczekB.MüllerS.DessimozC.GirstmairJ.et al (2015). A transcriptomic-phylogenomic analysis of the evolutionary relationships of flatworms.Curr. Biol.251347–1353. 10.1016/j.cub.2015.03.034
17
EickG. N.ThorntonJ. W. (2011). Evolution of steroid receptors from an estrogen-sensitive ancestral receptor.Mol. Cell Endocrinol.33431–38. 10.1016/j.mce.2010.09.003
18
EscrivaH.DelaunayF.LaudetV. (2000). Ligand binding and nuclear receptor evolution.Bioessays22717–727. 10.1002/1521-1878(200008)22:8<717::aid-bies5>3.0.co;2-i
19
EscrivaH. G.LaudetV.Robinson-RechaviM. (2003). Nuclear receptors are markers of animal genome evolution.J. Struct. Funct. Genomics3177–184. 10.1007/978-94-010-0263-9_17
20
FernándezR.GabaldónT. (2020). Gene gain and loss across the metazoan tree of life.Nat. Ecol. Evol.4524–533. 10.1038/s41559-019-1069-x
21
FonsecaE.MachadoA. M.Vilas-ArrondoN.Gomes-dos-SantosA.VeríssimoA.EstevesP.et al (2020). Cartilaginous fishes offer unique insights into the evolution of the nuclear receptor gene repertoire in gnathostomes.Gen. Comp. Endocr.295:113527. 10.1016/j.ygcen.2020.113527
22
GarcíaH. E.LaudetV.Robinson-RechaviM. (2003). Nuclear receptors are markers of animal genome evolution.J. Struct. Funct. Genomics3, 177–184.
23
HajiahmadiZ.AbediA.WeiH.SunW.RuanH.ZhugeQ.et al (2020). Identification, evolution, expression, and docking studies of fatty acid desaturase genes in wheat (Triticum aestivum L.).BMC Genomics21:778. 10.1186/s12864-020-07199-1
24
HallB. G. (2005). Comparison of the accuracies of several phylogenetic methods using protein and DNA sequences.Mol. Biol. Evol.22:792. 10.1093/molbev/msi066
25
HerathV.VerchotJ. (2021). Insight into the bZIP gene family in solanum tuberosum: genome and transcriptome analysis to understand the roles of gene diversification in spatiotemporal gene expression and function.Int. J. Mol. Sci.22:253. 10.3390/ijms22010253
26
HuangW.XuF.LiJ.LiL.QueH.ZhangG. (2015). Evolution of a novel nuclear receptor subfamily with emphasis on the member from the Pacific oyster Crassostrea gigas.Gene567164–172. 10.1016/j.gene.2015.04.082
27
HwangD. S.LeeB. Y.KimH. S.LeeM. C.KyungD. H.OmA. S.et al (2014). Genome-wide identification of nuclear receptor (NR) superfamily genes in the copepod Tigriopus japonicus.BMC Genomics15:993. 10.1186/1471-2164-15-993
28
HydeC. J.FitzgibbonQ. P.ElizurA.SmithG. G.VenturaT. (2019). Transcriptional profiling of spiny lobster metamorphosis reveals three new additions to the nuclear receptor superfamily.BMC Genomics20:531. 10.1186/s12864-019-5925-5
29
JinK. M.ZhuoR. Y.XuD.WangY. J.FanH. J.HuangB. Y.et al (2020). Genome-wide identification of the expansin gene family and its potential association with drought stress in Moso bamboo.Int. J. Mol. Sci.21:9491. 10.3390/ijms21249491
30
KaurS.JoblingS.JonesC. S.NobleL. R.RoutledgeE. J.LockyerA. E. (2015). The nuclear receptors of biomphalaria glabrata and lottia gigantea: implications for developing new model organisms.PLoS One10:e0121259. 10.1371/journal.pone.0121259
31
KimD. H.KimH. S.HwangD. S.KimH. J.HagiwaraA.LeeJ. S.et al (2017). Genome-wide identification of nuclear receptor (NR) genes and the evolutionary significance of the NR1O subfamily in the monogonont rotifer Brachionus spp.Gen. Comp. Endocr.252219–225. 10.1016/j.ygcen.2017.06.030
32
KongX.LvW.JiangS.ZhangD.CaiG.PanJ.et al (2013). Genome-wide identification and expression analysis of calcium-dependent protein kinase in maize.BMC Genomics14:433. 10.1186/1471-2164-14-433
33
LallemandT.LeducM.LandèsC.RizzonC.LeratE. (2020). An overview of duplicated gene detection methods: why the duplication mechanism has to be accounted for in their choice.Genes11:1046. 10.3390/genes11091046
34
LazarM. A. (2017). Maturing of the nuclear receptor family.J. Clin. Invest.1271123–1125. 10.1172/JCI92949
35
LecroiseyC.LaudetV.SchubertM. (2012). The cephalochordate amphioxus: a key to reveal the secrets of nuclear receptor evolution.Brief. Funct. Genomics11156–166. 10.1093/bfgp/els008
36
LequimeS.FontaineA.Ar GouilhM.Moltini-ConcloisI.LambrechtsL. (2016). Genetic drift, purifying selection and vector genotype shape dengue virus intra-host genetic diversity in mosquitoes.Plos Genet.12:e1006111. 10.1371/journal.pgen.1006111
37
Martín-DuránJ. M.EggerB. (2012). Developmental diversity in free-living flatworms.Evodevo3:7. 10.1186/2041-9139-3-7
38
MazairaG. I.ZgajnarN. R.LotufoC. M.Daneri-BecerraC.SivilsJ. C.SotoO. B.et al (2018). The nuclear receptor field: a historical overview and future challenges.Nucl. Receptor. Res.5:101320.
39
MiglioliA.CanesiL.GomesI. D. L.SchubertM.DumollardR. (2021). Nuclear receptors and development of marine invertebrates.Genes12:83. 10.3390/genes12010083
40
MooreR. C.PuruggananM. D. (2005). The evolutionary dynamics of plant duplicate genes.Curr. Opin. Plant Biol.8122–128. 10.1016/j.pbi.2004.12.001
41
NguyenL. T.SchmidtH. A.von HaeselerA.MinhB. Q. (2015). IQ-TREE: a fast and effective stochastic algorithm for estimating Maximum-Likelihood phylogenies.Mol. Biol. Evol.32268–274. 10.1093/molbev/msu300
42
NoreñaC.DamboreneaC.BrusaF. (2015). “Chapter 10 – phylum platyhelminthes,” in Thorp and Covich’s Freshwater Invertebrates, 4th Edn, edsThorpJ. H.RogersD. C. (Boston: Academic Press), 181–203.
43
NovotnýJ. P.ChughtaiA. A.KostrouchováM.KostrouchováV.KostrouchD.KaššákF.et al (2017). Trichoplax adhaerens reveals a network of nuclear receptors sensitive to 9-cis-retinoic acid at the base of metazoan evolution.PeerJ.5:e3789. 10.7717/peerj.3789
44
OlefskyJ. M. (2001). Nuclear receptor minireview series.J. Biol. Chem.27636863–36864. 10.1074/jbc.r100047200
45
ParvathaneniR. K.DeLeoV. L.SpiekermanJ. J.ChakrabortyD.DevosK. M. (2017). Parallel loss of introns in the ABCB1 gene in angiosperms.BMC Evol. Biol.17:238. 10.1186/s12862-017-1077-x
46
ReitzelA. M.PangK.RyanJ. F.MullikinJ. C.MartindaleM. Q.BaxevanisA. D.et al (2011). Nuclear receptors from the ctenophore Mnemiopsis leidyi lack a zinc-finger DNA-binding domain: lineage-specific loss or ancestral condition in the emergence of the nuclear receptor superfamily?Evodevo2:3. 10.1186/2041-9139-2-3
47
Robinson-RechaviM.CarpentierA.-S.DuffraisseM.LaudetV. (2001). How many nuclear hormone receptors are there in the human genome?Trends Genet.17, 554–556. 10.1016/s0168-9525(01)02417-9
48
RonquistF.TeslenkoM.van der MarkP.AyresD. L.DarlingA.HöhnaS.et al (2012). MrBayes 3.2: efficient bayesian phylogenetic inference and model choice across a large model space.Syst. Biol.61539–542. 10.1093/sysbio/sys029
49
SamonteR. V.EichlerE. E. (2002). Segmental duplications and the evolution of the primate genome.Nat. Rev. Genet.365–72. 10.1038/nrg705
50
ThomsonS. A.BaldwinW. S.WangY. H.KwonG.LeBlancG. A. (2009). Annotation, phylogenetics, and expression of the nuclear receptors in Daphnia pulex.BMC Genomics10:500. 10.1186/1471-2164-10-500
51
TharpM. E.CollinsJ. J.NewmarkP. A. (2014). A lophotrochozoan-specific nuclear hormone receptor is required for reproductive system development in the planarian.Dev. Biol.396, 150–157. 10.1016/j.ydbio.2014.09.024
52
VogelerS.GallowayT. S.IsupovM.BeanT. P. (2017). Cloning retinoid and peroxisome proliferator-activated nuclear receptors of the Pacific oyster and in silico binding to environmental chemicals.PLoS One12:e0176024. 10.1371/journal.pone.0176024
53
VogelerS.GallowayT. S.LyonsB. P.BeanT. P. (2014). The nuclear receptor gene family in the Pacific oyster, Crassostrea gigas, contains a novel subfamily group.BMC Genomics15:369. 10.1186/1471-2164-15-369
54
WangL.GuoK.LiY.TuY.HuH.WangB.et al (2010). Expression profiling and integrative analysis of the CESA/CSL superfamily in rice.BMC Plant Biol.10:282. 10.1186/1471-2229-10-282
55
WangS.ZhuX. qCaiX. (2017). Gene duplication analysis reveals no ancient whole genome duplication but extensive small-scale duplications during genome evolution and adaptation of Schistosoma mansoni.Front. Cell Infect. Microbiol.7:412. 10.3389/fcimb.2017.00412
56
WangY.XuL.ThilmonyR.YouF. M.GuY. Q.Coleman-DerrD. (2017). PIECE 2.0: an update for the plant gene structure comparison and evolution database.Nucleic Acids Res.45(D1)D1015–D1020. 10.1093/nar/gkw935
57
WasikK.GurtowskiJ.ZhouX.RamosO. M.DelásM. J.BattistoniG.et al (2015). Genome and transcriptome of the regeneration-competent flatworm, Macrostomum lignano.Proc. Natl. Acad. Sci. U.S.A.11212462–12467. 10.1073/pnas.1516718112
58
WeathermanR. V.FletterickR. J.ScanlanT. S. (1999). Nuclear-receptor ligands and ligand-binding domains.Annu. Rev. Biochem.68, 559–581. 10.1146/annurev.biochem.68.1.559
59
WuA.HaoP.WeiH.SunH.ChengS.ChenP.et al (2019). Genome-wide identification and characterization of glycosyltransferase family 47 in cotton.Front. Genet.10:824. 10.3389/fgene.2019.00824
60
WuW.LoVerdeP. T. (2019). Nuclear hormone receptors in parasitic Platyhelminths.Mol. Biochem. Parasit.233:111218. 10.1016/j.molbiopara.2019.111218
61
WuW.NilesE. G.El-SayedN.BerrimanM.LoVerdeP. T. (2006). Schistosoma mansoni (Platyhelminthes, Trematoda) nuclear receptors: sixteen new members and a novel subfamily.Gene366303–315. 10.1016/j.gene.2005.09.013
62
WuW.NilesE. G.HiraiH.LoVerdeP. T. (2007). Evolution of a novel subfamily of nuclear receptors with members that each contain two DNA binding domains.BMC Evol. Biol.7:27.
63
WudarskiJ.EggerB.RammS. A.SchärerL.LadurnerP.ZadesenetsK. S.et al (2020). The free-living flatworm Macrostomum lignano.Evodevo11:5. 10.1186/s13227-020-00150-1
64
XuG.GuoC.ShanH.KongH. (2012). Divergence of duplicate genes in exon–intron structure.Proc. Natl. Acad. Sci. U.S.A.1091187–1192. 10.1073/pnas.1109047109
65
YangP. J.ChenE. H.SongZ. H.HeW.LiuS. H.DouW.et al (2020). Molecular characterization and expression profiling of nuclear receptor gene families in oriental fruit fly, Bactrocera Dorsalis (Hendel).Insects11:126. 10.3390/insects11020126
66
ZadesenetsK. S.JetybayevI. Y.SchärerL.RubtsovN. B. (2020). Genome and karyotype reorganization after whole genome duplication in free-living flatworms of the genus Macrostomum.Int. J. Mol. Sci.21:680. 10.3390/ijms21020680
67
ZadesenetsK. S.SchärerL.RubtsovN. B. (2017). New insights into the karyotype evolution of the free-living flatworm Macrostomum lignano (Platyhelminthes, Turbellaria).Sci. Rep.76066.
68
ZhangY. M.ZhangS. Z.ZhengC. C. (2014). Genomewide analysis of LATERAL ORGAN BOUNDARIES Domain gene family in Zea mays.J. Genet.9379–91. 10.1007/s12041-014-0342-7
69
ZhangZ.BurchP. E.CooneyA. J.LanzR. B.PereiraF. A.WuJ.et al (2004). Genomic analysis of the nuclear receptor family: new insights into structure, regulation, and evolution from the rat genome.Genome Res.14, 580–590. 10.1101/gr.2160004
70
ZhengS.LongJ.LiuZ.TaoW.WangD. (2018). Identification and evolution of TGF-β signaling pathway members in twenty-four animal species and expression in Tilapia.Int. J. Mol. Sci.19:1154. 10.3390/ijms19041154
Summary
Keywords
nuclear receptor, Macrostomum lignano, flatworms, gene family expansion, intron loss, reduced intron size, evolution
Citation
Cheng Y, Chen J, Mukhtar I and Chen J (2021) Genome-Wide Characterization of the Nuclear Receptor Gene Family in Macrostomum lignano Imply Its Evolutionary Diversification. Front. Mar. Sci. 8:653447. doi: 10.3389/fmars.2021.653447
Received
14 January 2021
Accepted
07 May 2021
Published
14 June 2021
Volume
8 - 2021
Edited by
Oleg Simakov, University of Vienna, Austria
Reviewed by
Adam Michael Reitzel, University of North Carolina at Charlotte, United States; Vincent Laudet, Okinawa Institute of Science and Technology Graduate University, Japan
Updates

Check for updates
Copyright
© 2021 Cheng, Chen, Mukhtar and Chen.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jianming Chen, chenjm@mju.edu.cn
This article was submitted to Marine Molecular Biology and Ecology, a section of the journal Frontiers in Marine Science
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.