Multiple origins of green blood in New Guinea lizards.

Several species of lizards from the megadiverse island of New Guinea have evolved green blood. An unusually high concentration of the green bile pigment biliverdin in the circulatory system of these lizards makes the blood, muscles, bones, tongue, and mucosal tissues bright green in color, eclipsing the crimson color from their red blood cells. This is a remarkable physiological feature because bile pigments are toxic physiological waste products of red blood cell catabolism and, when chronically elevated, cause jaundice in humans and all other vertebrates. Although these lizards offer a promising system to examine the evolution of extraordinary physiological characteristics, little is known about the phylogenetic relationships of green-blooded lizards or the evolutionary origins of green blood. We present the first extensive phylogeny for green-blooded lizards and closely related Australasian lizards using thousands of genomic regions to examine the evolutionary history of this unusual trait. Maximum likelihood ancestral character state reconstruction supports four independent origins of green blood. Our results lay the phylogenetic foundation necessary to determine the role, if any, of natural selection in shaping this enigmatic physiological trait as well as understanding the genetic, proteomic, and biochemical basis for the lack of jaundice in those species that have independently evolved green blood.

INTRODUCTION

Understanding the evolutionary history of physiological adaptations is of critical importance to biology. In particular, unusual or outlier physiologies may provide important evolutionary insight into the breadth, tempo, and mode of key innovations. One of the most unusual amniote physiological traits is green blood that has evolved in several ecologically diverse species of New Guinea lizards (, ). Whereas most vertebrates have crimson red blood derived from the ferrous form of heme, these lizards have lime green–colored blood, resulting in bright green coloration of their muscles, bones, tongue, and mucosal tissues (). The green coloration is due to exceptionally high concentrations of the bile pigment biliverdin (714 to 1020 μM) in the circulatory system, which eclipses the vivid red color of their nucleated red blood cells (). In other vertebrates, bile pigments (biliverdin and bilirubin) remain generally undetectable in the circulatory system because their accumulation is known to be genotoxic, cytotoxic, neurotoxic, and causes jaundice (–). Bile pigments are normal degradation products of heme catabolism. As red blood cells senesce, degradation of the heme constituent of red blood cell hemoglobin is catalyzed by the enzyme heme oxygenase, resulting in biliverdin and the recycling of iron and carbon monoxide. In addition to immediate cellular damage, sustained biliverdin accumulation is associated with permanent nerve damage and developmental disorders (, ). Furthermore, recent research has shown that babies treated for jaundice are more likely to develop autism (). Bile pigment accumulation is commonly seen in individuals with impaired metabolic and liver function, such as newborn babies or patients with the life-threatening Gilbert’s and Crigler-Najjar syndromes (–). Nevertheless, several New Guinea lizards have the highest concentration of blood biliverdin measured for any animal () without being jaundiced. In contrast, the highest recorded level of biliverdin measured in humans, 50 μM, resulted in death (). We use the term hyperbiliverdinism for this nonpathological biliverdin accumulation. The unusual physiology found in these lizards offers a promising system to examine the underlying genomic and proteomic changes that allow these lizards to remain jaundice-free. However, the lack of a well-resolved phylogeny for these lizards and close relatives has prevented an understanding of the evolutionary history of hyperbiliverdinism.

Five described species of lizards with green blood are currently grouped into the same genus (Prasinohaema) based solely on blood coloration despite the fact that the five species are morphologically and ecologically divergent. We use phylogenomic data from thousands of genome-wide loci to establish a robust phylogeny for Prasinohaema and related Australasian lizards with normal blood to determine the evolutionary history of green blood. Our phylogenetic and ancestral character state reconstruction results produce a highly robust phylogeny that supports four independent origins of green blood, demonstrating surprising evolutionary dynamism of green blood. Our findings lay the foundation for determining the role of natural selection, if any, in shaping this enigmatic physiological feature in these lizards as well as understanding the genetic and biochemical basis for hyperbiliverdinism. Understanding how these lizards avoid jaundice with such high bile pigment levels may translate into nontraditional approaches to addressing specific practical health problems, providing valuable benefits to society.

RESULTS

Phylogenetic analyses

We investigated the phylogenetic relationships among 51 species (119 individuals) of Australasian skinks including 6 species with green blood, 2 of which are undescribed new species (table S1). We targeted 5060 ultraconserved elements (UCEs) for phylogenomic analyses and recovered an average of 2690 loci (100 to 4241) per individual with a mean length of 586 base pairs (bp) (224 to 5494) per locus (Table 1). We subsampled our full data set to examine the effects of number of loci on phylogeny estimation (our “1K” data set), to explore the influence of loci versus taxonomic sampling on phylogeny estimation (our “TWO” data set), and to examine the effects of uneven taxonomic sampling (our “ONE” data set). The phylogeny presented in the main text, from the ONE data set, was trimmed to include only one representative per species, resulting in 4333 loci shared among more than 50% of taxa. Concatenation and coalescent approaches to species tree inference recovered a well-supported topology (Fig. 1). Overall, highly supported phylogenetic relationships were concordant among our data sets, and any incongruencies had low bootstrap support (figs. S1 to S7). Our data were also robust to biases caused by missing data because analysis of multiple data sets with varying levels of matrix completeness and information content recovered highly concordant topologies within each data set (figs. S1 to S7).

Table 1

Summary of UCE alignments.

Ten samples with fewer than 499 (10%) UCE loci or missing voucher information were removed from our “ALL” data set. Our 1K data set removed all individuals with fewer than 1000 UCE loci, the TWO data set trimmed sampling to no more than two individuals per species, and the ONE data set removed all but one individual per species.

Data set	Taxa	Total length (bp)	No. UCE loci in % taxa
			50%	60%	70%	80%
ALL	109	1,545,997	4088	3331	1931	542
1K	98	1,641,106	4354	3952	3195	1624
TWO	63	1,748,774	4333	3946	3137	1606
ONE	43	1,838,540	4368	3986	3220	1910

Fig. 1

Ancestral reconstruction indicates four origins of green blood.

(A) Transitions between red blood and green blood in Australasian lizards were summarized from 1000 stochastic character simulations on a fixed species topology. The phylogeny was estimated from a 70% complete sequence matrix from our ONE data set using an unpartitioned concatenated analysis in RAxML. Bootstrap values for all nodes are 100. Branches are colored continuously according to posterior probability (PP) support for red (PP = 1.0) or green (PP = 0.0) blood under our “mixed” model. Pie charts indicate proportion of character histories that reconstructed red or green blood as the ancestral state for the most recent common ancestor of all Prasinohaema species under each transition rate model: all rates different (ARD), equal rates (ER), irreversible (IR), or mixed. Values for the character rate matrix were either fixed at their maximum likelihood estimates (MLE) or sampled from a posterior distribution of rate matrices (MCMC). Pictures show the following green-blooded species: (B) Prasinohaema flavipes, (C) P. prehensicauda, (D) Prasinohaema semoni, (E) Prasinohaema sp. nov., and (F) Prasinohaema virens.

Monophyly assessment

Surprisingly, our phylogenomic analysis strongly supports the nonmonophyly of green-blooded lizards; in addition, two other genera (Lipinia and Scincella) are not monophyletic. Our sampling included many Australasian taxa; however, the clade containing the most recent common ancestor for all green-blooded lizards is composed only of New Guinea species, indicating a Papuan origin for this radiation. We recover a robust fully resolved phylogeny with short internal branches that are likely the result of a rapid radiation. Given the large number of taxa and loci in our data set and the absence of particular loci in different taxa, Bayesian concordance analyses were not computationally feasible; however, statistical tests for monophyly corroborate the polyphyletic origins of the six Prasinohaema species. Shimodaira-Hasegawa (SH) tests on our ONE and TWO concatenated data sets supported the best unconstrained maximum likelihood (ML) tree over any of the constraint trees that forced either partial or complete Prasinohaema monophyly (). We explored the possibility of hemiplasy by searching for gene tree histories that support fewer than four lineages with green blood and found no support for Prasinohaema monophyly among more than 4000 individually estimated gene trees (tables S2 and S3) (). In addition, modeling for incomplete lineage sorting using a coalescent approach (ASTRAL) yielded topologies that were concordant with our unpartitioned ML approach (RAxML), indicating that our data were robust to problems of incomplete lineage sorting. Despite strong support for four green-blood lineages (Fig. 1), the number of evolutionary transitions between red and green blood is still unresolved: The evolution of green blood in lizards could be explained by four independent origins, a single origin with multiple reversals, or a complex pattern involving multiple gains and multiple losses.

Ancestral state reconstruction

We conducted an ML ancestral state reconstruction and stochastic character mapping of green blood on the unpartitioned concatenated ML species trees to investigate the most likely character state for blood physiology in the ancestor of all species with green blood. Model fitting for both ML reconstructions and character mappings resulted in similar support for alternative models (table S4). When the transition rate matrix was fixed at its empirical ML value, each model reconstructed red blood as the ancestral state with varying levels of support (27, 99, 100, and 84% probability) for both marginal reconstructions and stochastic mappings. Although the transition rates and mean number of reversals varied considerably between models, all reconstructions recovered an average of 4 (3.9 to 4.2) transitions from red to green blood (Table 2). In contrast, sampling a new transition rate matrix for each character history from a model’s posterior distribution of 1000 simulations using MCMC unequivocally supported (100% probability) red blood as the ancestral state and four independent origins of green blood. Model choice did not significantly change the number or placement of transitions using the MCMC approach but did vary the results of the ML approach, likely due to the use of a single rate matrix for all reconstructions under a given model (Fig. 1 and Table 2). All reconstructions presented in the main text were performed using our “ONE-p70” phylogeny as input, but we recovered similar results using all subsets of our ONE and TWO data sets (table S5). Furthermore, stochastic character mapping across randomly sampled subsets of posterior probability trees and gene trees never yielded less than four transitions from red to green blood for any given character history (table S6). These results suggest that our ancestral reconstruction and character mapping results are robust to hemiplasy caused by phylogenetic uncertainty. Together, our analyses indicate that hyperbiliverdinism is likely the result of repeated evolution.

Table 2

Integrating over multiple models for stochastic character mapping.

Three transition rate models (ARD, ER, and IR) and one mixed model were considered for stochastic character mapping. The number of simulations (Nsim) for each model was normalized using Akaike information criterion weights (wAIC) for that model and the mixed model sampled character histories from all three models. The transition rate matrix for a given model was either fixed at its empirical MLE or sampled using MCMC from its posterior distribution (MCMC). Rates (q), AIC scores, Likelihood (LnL), mean number of forward and reverse transitions, and PP of red blood reconstructed for the most recent common ancestor (MRCA) for all Prasinohaema are shown for the ONE-p70 tree. Models are described in the main text.

Transition rate matrix	Model	Model summary						Mean no. of changes		PP of red-blooded MRCA
		qRG	qGR	AIC	wAIC	ln (L)	Nsim	Red→Green	Red←Green
MLE (fixed)	ARD	11.7	137.2	32.2	0.386	−14.1	386	4.0	9.0	27.1
	ER	11.8	11.8	32.5	0.332	−15.3	332	4.2	0.5	99.9
	IR	10.6	0.0	32.8	0.282	−15.4	282	4.0	0.0	100
	Mixed	—	—	—	—	—	1000	3.9	2.2	81.4
MCMC	ARD	3.5	1.0	36.3	0.156	−16.1	156	3.9	0.1	100
	ER	3.8	3.8	33.5	0.622	−15.7	622	3.9	0.1	100
	IR	3.9	0.0	35.5	0.222	−16.7	222	4.0	0.0	100
	Mixed	—	—	—	—	—	1000	3.9	0.1	100

DISCUSSION

Surprisingly, we find robust support for the nonmonophyly of green-blooded lizards, revealing a complex and dynamic evolutionary history of this unusual physiological trait. This conclusion is supported by extensive phylogenomic analyses of 5060 UCE loci and tests for monophyly. Considering the pathological consequences of elevated concentrations of bile pigments in most animals, green blood was assumed to have originated only once in amniotes. However, we reconstructed red blood as the most likely ancestral state for all Prasinohaema lineages with strong support for four separate origins of green blood. Although these lizards are the only amniotes to have hyperbiliverdinism, several other vertebrates and invertebrates have been shown to have moderately elevated levels of bile pigments. Elevated levels of biliverdin have been documented in the blood and scales of several unrelated species of fishes (–) and insects (, ), and biliverdin is deposited in eggshells of some birds (). Biliverdin has been suggested as the cause of green bones, skin, and blood in several unrelated frogs (, ). The rare but widespread occurrence of moderately elevated levels of bile pigments in animals suggests that there may be some adaptive value of elevated biliverdin.

How cytotoxicity of bile pigments might be selected for is unclear; however, bile pigments have been shown to reduce pathogen activity (–), so elevated levels of bile pigments may preclude or reduce infection. Because lizards are the only amniotes that have this unusual physiology and remain healthy, further investigation will inform us how these lizards are not jaundiced. In fish with green blood plasma, excess biliverdin is tightly bound to albumin transport proteins, and this complex is neither metabolized nor excreted by the liver but remains in the bloodstream, which may obviate damage from biliverdin (). This mechanism for biliverdin retention and tolerance in unrelated fish is the result of convergent changes in the serum albumin protein, lending support to the possibility of green blood being an advantageous trait, shaped by natural selection (). Convergent and parallel evolution of unusual phenotypes may be more frequent than previously thought due to shared constraints in developmental patterns, genomic architecture, and body form ().

In contrast to the toxicity aspect of bile pigments, recent data on biosynthetic pathways indicate possible positive functional roles for bile pigments, but at low concentrations (). Excess biliverdin appears to have antioxidant and cytoprotective properties in vitro (, ), predominantly due to their ability to scavenge free radical species (, –). Bile pigments also appear to have a beneficial role in oxidative stress defense in animals () and have been identified as having significant neuroprotective antioxidant and antimutagenic properties (, ).

Our analyses indicate four independent origins of green blood from a normal, red-blooded ancestor with no losses. The primary alternative hypothesis is that there was only a single origin of green blood followed by four or more losses, an evolutionary process presumably driven by natural selection. To fully tease apart these hypotheses, our phylogenetic results lay the groundwork for future functional genomic and proteomic approaches to identify the genes responsible for the origins, and possible loss, of hyperbiliverdinism. We predict that independent gains will be identified by four unique genomic changes whereas a single origin will be identified by a shared genomic change and that we will be able to track the signature of independent reversals to red blood, if any, from genomic data. The exciting future challenges of this nonmodel system are to address the potential convergent and parallel evolution of green blood in these lizards, determine whether bile pigment retention confers any major selective advantages for cells in vivo, and examine potential proteomic or biochemical innovations associated with nonpathological biliverdin accumulation.

METHODS

DNA sampling and sequencing

We took a genome-wide next-generation sequencing approach and targeted 5060 UCE loci using an in-solution sequence capture protocol (). Our sampling included 24 lizards with green blood and 95 closely related red-blooded individuals. Genomic DNA was extracted from tissue preserved in ethanol using a standard DNA salt extraction protocol () and sonicated into sizes between 400 and 600 bp, with a target of 400 bp using an EpiSonic 1000 (EpiGentek). The Kapa Hyper Prep Kit (Kapa Biosystems Inc.) was used to prepare libraries for target enrichment and sequencing, including ligation of custom dual-indexed adapters with unique barcode sequencing for each individual. Prior to sequence capture, we pooled indexed libraries into groups of 8 in equimolar ratios. For sequence capture, we used the Tetrapods-UCE-5Kv1 probe set (myBaits, MYcroarray) targeting 5060 UCE loci following a freely available open-source protocol. This approach has been successfully used in lizards and other vertebrates (, ). Following enrichment and limited-cycle polymerase chain reaction, pools of equimolar reactions were combined and sent to Oklahoma Medical Research Foundation (Oklahoma City, OK) for 150-bp paired-end sequencing on an Illumina HiSeq3000.

Data processing

To process demultiplexed raw reads for phylogenomic analyses, we followed the standard PHYLUCE pipeline (). Illumiprocessor, in combination with Trimmomatic, was used to filter low-quality reads and remove adapters (, ), yielding an average of 8,504,473 reads (302,468 to 119,426,985) and 1,141,548,626 bp (42,543,058 to 17,007,264,440) per sample. Filtered reads were assembled into contigs using Trinity (). To check for contamination, all contigs were aligned to the nonredundant protein database by BLASTP with an e value cutoff of 1 × 10−6 and checked for best hits to squamate reptiles. We also performed an all-versus-all blast search on the assemblies and removed six individuals that had reciprocal best hits to a different species, likely the result of cross contamination. The resulting assemblies were screened for matches to the UCE probes, and nontargets and paralogs were removed in PHYLUCE. The resulting regions were aligned using MAFFT (), and quality-trimmed using a parallelized wrapper around Gblocks (). We created three data sets from our “all taxa” data set. First, to examine the effects of number of loci in phylogeny estimation, we removed 21 individuals with <1000 UCE loci to create a 98-taxon 1K data set. Second, to explore the influence of number of loci versus taxonomic sampling, we created a TWO data set, which restricted sampling to no more than two individuals per species, independent of loci number. Third, to examine the effects uneven taxonomic sampling, we constructed a ONE data set, for which we removed individuals with less than 1000 UCEs and limited taxon sampling to one representative per species. To assess the sensitivity of phylogenetic inference to missing data, we subsampled each data set and generated alignments that allowed up to 20, 30, 40, or 50% missing data (data sets are appended with p80, p70, p60, p50, respectively). Finally, we generated four additional alignments for each data set containing only the top 100, 250, 500, or 750 loci with the most parsimonious informative sites (data sets are appended with i100, i250, i500, or i750, respectively). In summary, each of our data sets was subsampled into eight alignments using two metrics: missing data and information content.

Phylogenetic inference

We performed unpartitioned concatenated ML analysis of each alignment and data set in RAxML v8.2 (), using 20 searches with GTR GAMMA model of nucleotide substitution for the best ML tree. Nodal support for the best ML topology was assessed using nonparametric bootstrapping with the autoMRE option in RAxML, and we reconciled the best ML tree with the bootstrap replicates in RAxML. Species-tree methods have been shown to be more accurate than concatenated approaches for multilocus data under a variety of conditions (). Because of the large number of gene trees that are typically produced from phylogenomic data sets, species trees were estimated using a summary coalescent approach. Individual gene trees were first estimated in an ML-like framework implemented in RAxML v8.2. Next, multilocus bootstraps were generated by sampling loci with replacement to create 100 resampled data sets, and we performed nonparametric bootstrapping by sites of each locus in each of the resampled data sets using RAxML. Then, bootstrapped gene trees were sorted to ensure that each resampled data set was resampled by loci and by sites (). The 100 resampled data sets were used to infer species trees in a summary coalescent framework, implemented in ASTRAL-II 4.7.8 () that allows unrooted, multifurcating trees as input and has been shown to be robust to effects of both incorrectly estimated gene trees and gene tree heterogeneity (). Extended majority-rule consensus (MRE) trees were estimated using sumtrees in DendroPy v3.2 (). After visually confirming species monophyly and topological concordance among data sets, we elected to proceed with the ONE and TWO data sets for all subsequent analyses because they had even taxonomic sampling, which resulted in a higher number of shared loci between species. For the main text, we presented the results from our 70% complete matrix which was trimmed to one representative per species (ONE-p70), because it was the most complete matrix without a significant reduction in the number of represented species or loci per individual (Table 1).

To assess support for the nonmonophyly of lizards with green blood, we took two approaches. First, we conducted SH tests () on each ML estimate, or the “best tree,” for each concatenated data set in RAxML. This was accomplished by running 20 ML searches for the best “constrained” tree by passing RAxML a constraint file with the -g option, then comparing likelihood scores using the -f h option, and passing RAxML the best constrained tree with -z and the best tree with -t. Second, we searched for both partial and full green blood monophyly in the individually estimated gene trees, excluding the Prasinohaema flavipes + Prasinohaema prehensicauda clade. For this, we used the is.monophyletic option in the R package ape and iterated over more than 4000 gene trees in our ONE and TWO data sets. Our results are summarized in tables S2 and S3.

We used ML ancestral state reconstruction and stochastic character mapping to assess whether hyperbiliverdinism can best be explained by multiple independent origins or a single origin with multiple reversals. Because branch length estimates are important for ancestral state reconstruction, all reconstructions were performed on species trees inferred from the unpartitioned concatenated alignments or on individual ML gene trees for each data set. We reconstructed ancestral states for green blood in an ML framework using the ace command in the R package APE v3.4 (). Three transition rate matrices were considered: (i) a two-parameter model that had different rates for forward and reverse transitions (ARD), (ii) a two-parameter model that assumed once green blood is gained there were no reversals to red blood (IR) model, and (iii) a single-parameter model that assumed forward and reverse transition rates were equal (ER). To determine the model fit, we examined the AIC for each model using the fitDiscrete command from the package geiger v2.0.6 (). Because the three models gave comparable support (table S4), we applied a fourth mixed model for which we integrated stochastic character maps across all three transition models. For the mixed model, we calculated normalized Akaike weights from the log likelihood of the other three models and sampled a number of character histories in proportion to the weight of evidence in support for each model. For example, the Akaike weights for the ONE-p70 data set were 0.386, 0.332, and 0.282, for the ARD, ER, and IR models, respectively. Thus, we simulated 386 character histories under the ARD model, 332 character histories under the ER model, and 282 character histories under the IR model. For all four models, we used stochastic character mapping (), implemented in phytools (), to summarize the results of multiple trait mappings onto our phylogeny using a continuous-time Markov process. The transition rate matrix estimated for each model was applied to these simulations using two approaches. First, we fixed the transition matrix at its empirical ML value estimated for a given model and generated multiple character histories. Second, we sampled a new transition matrix for each character history from the posterior distribution of a given model using MCMC. For the mixed-model MCMC approach, a new transition rate matrix was sampled from the posterior distribution of each model in proportion to the calculated Akaike weights. We generated 1000 character histories for each model using the make.simmap command and summarized the results using describe.simmap to determine the posterior probabilities that each internal node is at each state, the proportion of time spent in each state, and the number of transitions across the phylogeny. Results for the ONE-p70 topology are shown in Fig. 1 and Table 2. Results for other input trees are summarized in table S5.