Evolutionary Plasticity in Detoxification Gene Modules: The Preservation and Loss of the Pregnane X Receptor in Chondrichthyes Lineages

Abstract

To appraise how evolutionary processes, such as gene duplication and loss, influence an organism’s xenobiotic sensitivity is a critical question in toxicology. Of particular importance are gene families involved in the mediation of detoxification responses, such as members of the nuclear receptor subfamily 1 group I (NR1I), the pregnane X receptor (PXR), and the constitutive androstane receptor (CAR). While documented in multiple vertebrate genomes, PXR and CAR display an intriguing gene distribution. PXR is absent in birds and reptiles, while CAR shows a tetrapod-specific occurrence. More elusive is the presence of PXR and CAR gene orthologs in early branching and ecologically-important Chondrichthyes (chimaeras, sharks and rays). Therefore, we investigated various genome projects and use them to provide the first identification and functional characterization of a Chondrichthyan PXR from the chimaera elephant shark (Callorhinchus milii, Holocephali). Additionally, we substantiate the targeted PXR gene loss in Elasmobranchii (sharks and rays). Compared to other vertebrate groups, the chimaera PXR ortholog displays a diverse expression pattern (skin and gills) and a unique activation profile by classical xenobiotic ligands. Our findings provide insights into the molecular landscape of detoxification mechanisms and suggest lineage-specific adaptations in response to xenobiotics in gnathostome evolution.

1. Introduction

Nuclear receptors (NRs) are central constituents of animal endocrine systems. These ligand-dependent sensors act as transcription factors, regulating key physiological processes including metabolism, development, reproduction and nutrient utilization [1]. Importantly, NRs are also directly exploited by xenobiotics, causing numerous examples of physiological impairment (e.g., [2,3]). Two critical components of the vertebrate’s “chemical defensome” are the pregnane X receptor (PXR) and the constitutive androstane receptor (CAR) [4]. These are part of nuclear receptor subfamily 1 group I (NR1I), which also includes the vitamin D receptor (VDR). PXR and CAR were originally identified as xenobiotic sensors, since they regulate genes involved in drug metabolism such as phase I cytochrome P450 (e.g., CYP3A4, CYP2B6 and CYP2C), phase II transferases (e.g., uridine 5′-diphospho(UDP)-glucuronosyl transferase and glutathione-S-transferase), and drug transporters. Moreover, PXR is notoriously involved in other metabolic processes including energy homeostasis, inflammatory responses, cell proliferation, apoptosis, and tumour development [5,6,7]. The taxonomic distribution of VDR/PXR/CAR gene orthologs is remarkably mutable [8,9,10,11]. In vertebrate species, VDR is found in both cyclostomes (lampreys) and gnathostomes [12]; CAR occurs in tetrapods [10,13]; while PXR genes have been described and characterized in amphibians [10] and mammals [10]. On the other hand, teleost genomes, such as zebrafish (Danio rerio), also retain PXR, but this is not an universal condition found throughout teleost lineages [10,11], consistent with the highly derived nature of their genomes [14]. Importantly, synteny supports the hypothesis that the absence of PXR in birds, reptiles, and some teleosts, as well as CAR in ray-finned fish, is due to secondary gene loss [10]. Genome comparisons between human and teleost PXR, CAR, and VDR orthologous genomic regions further implicates whole genome duplications (2R) as the underlying cause of the NR1I gene expansion [10,15]. These observations suggest that VDR, PXR, and CAR first appeared in the ancestors of vertebrates and should be present in early lineage genomes such as Chondrichthyes (cartilaginous fish). Consistently, VDR has been described and functionally characterized in cartilaginous fishes [12]. PXR and CAR orthologs have not been described in Chondrichthyes, although the presence of the former has been suggested [10]. Here we thoroughly investigate the gene repertoire of the NR1I subfamily, central components of the “chemical defensome”, in Chondrichthyes. Cartilaginous fishes are divided into two branches: Holocephali (chimaeras) and Elasmobranchii (sharks, rays and skate). Together, they are a highly diversified group of early branching vertebrates, representing important components of aquatic ecosystems and food webs, and thus are key ecological indicators [16,17].

2. Results

2.1. Identification of nuclear receptor subfamily 1 group I (NR1I) Ortholog Genes in Chondrichthyes

To determine the gene complement of VDR/PXR/CAR-like genes in Chondrichthyes species, we examined five genome datasets from two subclasses: Holocephali (chimaeras) and Elasmobranchii (sharks and rays) [18,19,20,21,22]. Our search identified one or two genes with similarities to VDR/PXR/CAR genes in the elephant shark (Callorrhyncus milii), little skate (Leucoraja erinacea), cloudy catshark (Scyliorhinus torazame), brownbanded bamboo shark (Chiloscyllium punctatum), and whale shark (Rhyncodon typus). To establish the orthology of the retrieved sequences, we performed phylogenetic analysis (Figure 1A), including the described VDR/PXR/CAR gene sequences from mammals, birds, reptiles, amphibians, teleosts, lepisosteiformes, cyclostomes, and tunicates (Table S1). Three statistically-supported sequence clades were retrieved, consistent with the separation into VDR, PXR and CAR genes (Figure 1A). We then compared the two critical functional domains of the CmiPXR: the DNA-binding domain (DBD) and the ligand binding domain (LBD). CmiPXR shares around 70% of its sequence identity to other PXR-DBDs (Figure 1B), with the LBD identity values displaying significantly lower values (Figure 1B).

2.2. Synteny Analysis of NR1I Ortholog Genes

To further verify the orthology of these novel gene sequences and to discriminate between true gene loss or absence of sequencing data, we next verified the genomic location of VDR, PXR, and CAR within the syntenic locations in the genomes of the elephant shark, cloudy catshark, brownbanded bamboo shark, and whale shark, using the human and zebrafish gene loci composition as reference, as shown in Figure 2. The scattered assembly and small contiguous size of the current little skate genome (LER_WGS_1—GCA_000238235.1) impeded a consistent comparative analysis at this stage. The PXR gene from the elephant shark is flanked by MAATS1 and GSK3B genes (scaffold NW_006890095.1 3.95Mb), and the overall locus composition is similar to that of other vertebrate species (Figure 2). In the Elasmobranchii species analysed here, and despite the global synteny conservation, no intervening PXR-like sequence was found between MAATS1 and GSK3B (brownbanded bamboo shark—scaffold scf_chipu00000056; whale shark—scaffold scf_rhity00002454; cloudy catshark—scaffold scf_scyto00010339) (Figure 2, Tables S2–S4). In the case of the VDR locus, we searched the scaffolds containing the human VDR flanking genes TMEM106C and HDAC7, but no HDAC7 ortholog was found in any of the Chondrichthyes genomes. In the elephant shark, VDR was found in the same scaffold as TMEM106C (scaffold NW_006890370.1 105.7kb), contrary to Elasmobranchii VDRs (brownbanded bamboo shark—scaffold scf_chipu00001415 28.2kb; whale shark—scaffold NW_018047310.1 2.9kb; cloudy catshark—scaffolds scf_scyto00007144 and scf_scyto00012969), which were found on different scaffolds than the TMEM106C orthologs (brownbanded bamboo shark—scaffold scf_chipu00001599; whale shark—scaffold NW_018032445.1; cloudy catshark—scaffold scf_scyto00007676). This was probably due to missing sequencing data for the intervening genomic region (Figure 2, Tables S2–S4). In Chondrichthyes, the CAR locus is dispersed in comparison to humans (Figure 2, Tables S2–S4).

2.3. Gene Expression Analysis of the Elephant Shark Pregnane X Receptor (PXR)

Next, we investigated the gene expression profile of the CmiPXR in a tissue panel. Our analysis indicated that the elephant shark gene ortholog displayed a unique pattern, with restricted expression in the skin and gills (Figure 3).

2.4. Transactivation Assays of CmiPXR

Given the differences in the LBD sequence of CmiPXR, we explored the capacity to transactivate gene expression in the presence of classical PXR ligands from different chemical categories: the natural and the synthetic steroid hormone 17β-estradiol (E2) and 17α-ethinylestradiol (EE2) respectively, and the environmental contaminants trans-nonachlor (TNC) and bisphenol A (BPA). A mammalian cell-based activation assay and the zebrafish PXR (DrePXR) were used as control. Both E2 and EE2 significantly activated (p < 0.05) DrePXR and CmiPXR at high concentrations (Figure 4). Regarding the effect of the two environmental pollutants tested, both DrePXR and CmiPXR were significantly activated (p < 0.05) when exposed to the highest tested concentration of TNC and the two highest concentrations of BPA (Figure 4).

3. Discussion

By performing a comprehensive search of the NR1I gene repertoire in early branching gnathostome genomes, we unfolded the evolutionary history of this fundamental component of detoxification response. Notably, we were able to deduce that PXR, while present in the elephant shark, a Holocephali, has been most likely lost in the investigated Elasmobranchii species. The two newly-identified genes in the elephant shark fall into the PXR and VDR classes, while the single gene identified in little skate, cloudy catshark, brownbanded bamboo shark, and whale shark are bona fide VDR orthologs. The orthology of the new gene sequences found in this study was further confirmed by the syntenic analysis of VDR, PXR, and CAR locations in the genomes of the elephant shark, cloudy catshark, brownbanded bamboo shark, and whale shark. Regarding the CAR locus, its dispersed composition in Chondrichthyes compared to humans impedes a formal conclusion on the loss of CAR in these species, although this is the likeliest scenario. Additionally, gene orthologs of both PXR and CAR were not found in the two currently-available genomes of cyclostomes (sea lamprey, Pmar_germline 1.0 and Japanese lamprey, LetJap1.0). Furthermore, our analysis of amino acid identity between DBD and LBD of human, mouse, zebrafish, and elephant shark PXRs is consistent with previous studies for other species [23,24]. As expected, despite the substantial variation in sequence identity among vertebrates, we observed that the PXR–DBD is more conserved than LBD between species. This suggests that different PXRs should recognize similar response elements in the promotor of target genes, but their activation might be triggered upon binding to different ligands. The large and flexible ligand-binding pocket of PXR allows this receptor to accommodate a huge and diverse ligand range, such as endogenous ligands (5β-pregnane, progesterone, testosterone, lithocholic acids, and 17β-estradiol), antibiotics, drugs, carcinogens, and an array of environmental pollutants [25,26,27,28,29].

The occurrence of different gene complements of the NR1I subfamily in Chondrichthyes parallels similar findings in other vertebrate lineages [11,30,31]. Recently, an extensive investigation into 76 fish genomes suggests that approximately half of these species have lost PXR [11], in line with the description made here in cartilaginous fishes. Moreover, xenobiotic exposure experiments with classical xenobiotics, PXR ligands in cod (PXR-absent) did not show a clear transcription activation of gene coding for P450 cytochrome enzymes (CYP3A), as observed in mice or zebrafish [11]. In addition, promotor analysis raised the interesting possibility the PXR-absent species might have their CYP3A and CYP1A genes regulated by an unrelated xenobiotic-sensing transcription factor, the Aryl Hydrocarbon Receptor (AHR) [11]. Whether this is the case in Chondrichthyes that have lost PXR remains to be investigated. Together, these results suggest that the transcriptional regulation of detoxification gene modules is exceptionally plastic and has been rewired during vertebrate evolution.

The expression of the PXR gene in the elephant shark is confined to skin and gills, in clear contrast with what is observed in other species. For example, PXR gene transcripts in zebrafish are found in the liver, eye, intestine, brain, heart, and kidney, while in mammals, PXR is expressed in the liver, gastrointestinal tract, brain, and retina [23].

The ligand binding profile of the chimaera PXR ortholog is also puzzling, with CmiPXR exhibiting a smaller activation for the selected environmental pollutants, in contrast to PXR from zebrafish [32]. However, distinct sensitivities towards xenobiotics have been reported across species [32,33,34,35]. Thus, we cannot fully discard the existence of a distinct set of potential PXR-interacting xenobiotics, or other unidentified endogenous ligands, for chimaeras. Nonetheless, our transactivation results, together with the unique expression profile, raise the interesting possibility that CmiPXR could act as a specialized steroid-like sensor in the skin. In fact, previous studies have suggested putative effects of skin-expressed PXR in humans and rodents, including the induction of keratinocyte proliferation, immune hyper-responsiveness, modulation of DNA repair mechanisms, and overall skin barrier functions [36,37]. On the other hand, it is known that estrogens participate in skin homeostasis, by modulating collagen deposition, wound healing and scarring, and maintaining skin hydration and elasticity [38]. Furthermore, both PXR and estrogen have been directly or indirectly linked to fibrous connective tissue equilibrium [37,39]. Thus, we hypothesize that PXR could play a role in estrogen-dependent skin maintenance in chimaeras, contributing to the peculiar appearance of their smooth, rubbery and scale-less skin.

4. Material and Methods

4.1. Sequence Retrieval and Phylogenetic Analysis

Amino acid sequences were retrieved through blast searches in the publicly-available genome databases, using as reference annotated human VDR, PXR, and CAR sequences. Sequence sampling included major vertebrate lineages: mammals (Homo sapiens, Mus musculus, Rattus norvegicus, Oryctolagus cuniculus, Sus scrofa, Bos taurus), birds (Gallus gallus, Anas platyrhynchos), reptiles (Anolis carolinensis), amphibians (Xenopus tropicalis), euteleostei (Danio rerio, Cyprinus carpio, Oryzias latipes, Oreochromis niloticus), osteoglossomorpha (Scleropages formosus), holostei (Lepisosteus oculatus), chondrichthyans (Elasmobranchii: Leucoraja erinacea, Chiloscyllium punctatum, Scyliorhinus torazame, Rhyncodon typus; Chimaera: Callorhinchus milii), cyclostomes (Petromyzon marinus) and invertebrates (Ciona intestinalis). The retrieved sequences and corresponding protein accession numbers are listed in the Table S1. A multiple alignment of the retrieved sequences was obtained with MAFFT (multiple alignment using fast Fourier transform) alignment software [40] using default parameters. The final sequence alignment, containing 52 sequences and 659 positions, was used to construct a phylogenetic tree with MrBayes v 3.2.3 (CIPRES, San Diego, CA, USA; http://www.phylo.org/index.php/) sited in the CIPRES Science Gateway V3.3 [41]. The Bayesian analyses were performed under a mixed substitution model with two independent runs of four chains (one cooled and 3 heated) for 1 × 10⁶ generations, and the trees were sampled every 500 generations with a burnin set to 0.25 until the average standard deviation of the split frequencies remained <0.01. The statistical support for each branch is indicated at the nodes and expressed as Bayesian posterior probabilities [42]. FigTree v1.3.1 was used to visualize the tree. Geneious^® v7.1.7 (Biomatters Ltd., Auckland, New Zealand) was used to calculate the amino acid identity between human, mouse, zebrafish and elephant shark PXRs.

4.2. Synteny Analysis

The genomic region containing the human PXR, CAR, and VDR genes were localized at chromosomes 3 (119.78 Mb), 1 (161.22 Mb), and 12 (47.84 Mb), respectively. The human neighbouring genes, as well as the respective loci (GRCh38.p7), were collected from the GenBank database and used as references to assemble the synteny maps of zebrafish, elephant shark, cloudy catshark, brownbanded bamboo shark, and whale shark. To find the ortholog genes in the genomes of zebrafish (GRCz10), elephant shark (Callorhinchus_milii-6.1.3), and whale shark (ASM164234v2), we performed a BLAST of the human neighbour genes. Four flanking genes from both sides of each target gene were considered against the above-mentioned genomes. In the case of the cloudy catshark (Storazame_v1.0), brownbanded bamboo shark (Cpunctatum_v1.0), and the new version of whale shark (Rtypus_kobe_v1.0), the flanking genes found in elephant shark, as well the previous reference genes of human, were blasted (blast-n: -word_size 10, -outfmt 6, -num_threads 50) against the three recently-built elasmobranchs genomes (https://figshare.com/authors/Phyloinformatics_Lab_in_RIKEN_Kobe/4815111). Importantly, we used the new version of the whale shark genome (Rtypus_kobe_v1.0) to complete the previous syntenic map of the whale shark (ASM164234v2). Next, we manually inspected the BLAST-n results and, using the qstart, qend, sstart, send, and bit score options of outfmt6 format of BLAST software, reconstructed the structure for each gene. To confirm the neighbors’ homology in non-annotated genomes (C. punctatum, S. torazame, and R. typus new version), we used the following strategy: (1) .fasta and .gff files of each genome were used to extract the predicted coding region of each homolog candidate gene (if the blast approach detected the gene fragmentated in different scaffolds, we considered the biggest); (2) we performed reciprocal blast-n (with megablast and dc-megablast algorithm) searches of all candidates genes in the nucleotide database of NCBI (NT-NCBI); (3) if each candidate gene matched against the expected genes (references mentioned above, or ortholog genes in other species), it was kept and used to build the synteny maps.

4.3. Construction of Plasmid Vectors

The PXR hinge region and ligand binding domain (LBD) were isolated from zebrafish using a polymerase chain reaction (PCR) approach with the specific primers (restriction sites are underlined) F: 5′-atttCTAGAATGAAGAGAGAGCTGATCATGTC-3′ and R: 5′- aattGGTACCCTTTGTGAGGACTTAGGTGTC-3′, and the Phusion Flash master mix (Thermo Fisher Scientific, Waltham, MA, USA), according to the protocol from the supplier. The hinge and LBD of the elephant shark PXR was synthesized by IDT—Integrated DNA Technologies, Inc. (https://eu.idtdna.com/pages/home). Both PXRs were digested with XbaI and KpnI restriction enzymes (NZYtech) and ligated to pBIND (AF264722; Promega) with T4 ligase (Promega, Madison, WI, USA) to produce a GAL4-LBD “chimeric” receptor. The chimeric receptor produces a hybrid protein that contains the Gal4 DNA binding domain (DBD) and acts on an upstream activation sequence (UAS) response element. Plasmid sequences were confirmed using Sanger sequencing (Eurofins GATC, Constance, Germany).

4.4. Gene Expression

Total RNA was extracted from the following tissues of elephant shark using Trizol reagent (Gibco BRL): brain, gill, heart, intestine, kidney, liver, muscle, ovary, skin, and testis. One microgram of total RNA from each tissue was reverse-transcribed into cDNA with Superscript II (Invitrogen, Carlsbad, CA, USA) and used as a template for RT-PCR. The following primer pair which spans an intron was used to amplify elephant shark PXR: PXR_F, 5′-TGGAAGATCTCCTGGAAGCACATC-3′ and PXR_R, 5′-GAAGTTACGCTGGAGCTTGTAGTC-3′. Actin was amplified as an internal control to verify the integrity of cDNA using the primers: Actin_F, 5′-GGTATTGTCACCAACTGGGAC-3′ and Actin_R, 5′-AGATGGGCACAGTGTGGGTG-3. The PCR cycles comprised an initial denature step of 95 °C for 30 s, followed by 35 cycles of 95 °C for 10 s, 60 °C for 30 s, and 72 °C for 30 s, and a final elongation step of 72 °C for 5 min.

4.5. Transfection and Transactivation Assays

Cell culture and transactivation assays were performed as described in Fonseca et al. 2017 [43]. All ligands used (E2, EE2, TNC and BPA) were purchased from Sigma Aldrich (Sintra, Portugal). All compounds were resuspended in dimethyl sulfoxide (DMSO) at a final concentration of 0.1, 1, and 10 μM; and 5, 50, and 100 μM for BPA. Briefly, Cos-1 cells (Sigma, Sintra, Portugal) were maintained in DMEM (PAN-Biotech, Aidenbach, Bayern, Germany) supplemented with 10% fetal bovine serum (PAN-Biotech, Aidenbach, Bayern, Germany) and 1% penicillin/streptomycin (PAN-Biotech, Aidenbach, Bayern, Germany) at 37 °C with a humidified atmosphere and 5% CO₂. Cells were seeded in 24-well culture plates, and after 24 h, cells were transfected with 0.5 μg of pBIND constructs (pBIND-CmiPXRLBD or pBIND-DrePXRLBD) and 1 μg of pGL4.31[luc2P/GAL4UAS/Hygro] luciferase reporter vector (DQ487213; Promega, Madison, WI, USA) containing five UAS elements upstream of the firefly luciferase reporter gene, using lipofectamine 2000 reagent (Invitrogen, Carlsbad, CA, USA), in Opti-MEM (Gibco, Carlsbad, CA, USA), according to the manufacturer’s instructions. After 5 h of incubation, the transfection media was replaced with a medium containing the test compounds (E2, EE2; TNC—1, 1, and 10μM; and BPA—5, 50, and 100 μM) dissolved in DMSO (0.1%). Cells were lysed 24 h after transfection, and assayed for Firefly luciferase (reporter pGL4.31) and Renilla luciferase (pBIND) activities with Dual-Luciferase Reporter Assay System (Promega, Madison, WI, USA), according to the manufacturer’s instructions. All transfections were performed with two technical replicates per condition in three independent assays. The results were expressed as a fold-induction resulting from the ratio between luciferase (reporter pGL4.31) and Renilla (internal control for transfection efficiency luminescent activity), and then normalized by the DMSO control. Transactivation data was presented as means of the normalized values (n = 3) and bars with standard error of the mean (SEM) from the three separate experiments. The means of the technical replicates were used for statistical analysis with one-way analysis of variance (ANOVA), followed by the Holm–Sidak method in SigmaPlot 11.0 software (Systat Software Inc., San Jose, CA, USA). The level of significance was set to 0.05.

5. Conclusions

In this study, we deciphered the early evolution of the central components of the vertebrate “chemical defensome”. Our findings indicate that PXR gene orthologs are present in Holocephali but have been probably lost in Elasmobranchii. Moreover, the chimaera PXR gene displays a unique pattern of gene expression. Future studies will be required to dissect the molecular wiring of detoxification gene modules in Chondrichthyes.