Tomkins on the Human Vitellogenin Pseudogene: Who Needs Signal When You Have Noise?

Tomkins on the Human Vitellogenin Pseudogene: Who Needs Signal When You Have Noise?

I’ve blogged previously on the subject of pseudogenes and how they can constitute evidence for common ancestry. In short, they can represent the remnants of genes that once served a specific purpose in our distant ancestors, but were lost in a particular lineage for one reason or another. Imagine if there were a “gill gene” in all fish, and we found the broken remains of this gene in all tetrapods, including humans. That would be pretty compelling evidence that the ancestors of all tetrapods, including humans, had a functional copy of the gill gene and therefore had gills, suggesting that these ancestors could reasonably be called “fish”. Ergo, common descent between tetrapods and extant fish. That’s a cartoon example (there’s no such thing as a “gill gene”, embryological development isn’t nearly that simplistic) meant to illustrate the principle, but the real examples are not too dissimilar. For example, my aforementioned blog post discussed toothless mammals and birds possessing pseudogenes for enamel proteins, suggesting their descent from toothed ancestors (Meredith et al., 2013).

Another example is the case of the vitellogenin-1 pseudogene, which will be the centrepiece of this blog post. Vitellogenin-1 (which I’ll call “VTG”) is a gene that encodes a protein by the same name that is a precursor of the nutritious egg yolk proteins produced by nearly all oviparous (egg laying) animals. According to evolution, the ancestors of Therian (placental+marsupial) mammals also laid eggs, as we’re ultimately descended from reptile-like creatures. So our ancestors would have likely had a functioning VTG gene. We no longer produce egg yolk to nourish our offspring, so this gene must have lost its function along the way. One way we can confirm this then, would be to find the remains of the VTG gene in our genomes, as these remains would prove that we’re descended from creatures that produced VTG – ostensibly egg-laying animals.

This confirmation came in 2008, in a paper by Brawand et al. (2008). For more information about this, see the “vitellogenin test” section of this 2012 Biologos article by Dr. Dennis Venema.

However, not everyone agrees with this evidence or the conclusion. Enter Dr. Jeffrey Tomkins, the “Director of Life Sciences” at the Institute for Creation Research (ICR). Tomkins received his PhD in Genetics in 1996 from Clemson University, and worked primarily in plant genetics performing BAC sequencing until he joined ICR in 2009 (incidentally, the same year as Nathaniel Jeanson). He’s been a major contributor to various “creation science” journals, publishing nearly 40 “papers” in the last 8 years. The main topics of Tomkins’ published “creation science” articles are human/chimp similarity, human chromosome 2, and pseudogenes. It’s obviously that last one that’s relevant today.

In late 2015, Tomkins published an article in “Answers Research Journal” (Answers in Genesis’ Journal) entitled “Challenging the BioLogos Claim that a Vitellogenin (Egg-Laying) Pseudogene Exists in the Human Genome“. It was short, just ~3,000 words, but in it, Tomkins claims to have provided definitive evidence that the sequence purported to represent the human vitellogenin pseudogene is in fact not a pseudogene at all, but a functional enhancer element of a gene involved in neurological processes in the brain.

I’d been peripherally aware of this article for a couple of years now, but I never sat down and really read through it with a critical eye. That was until a few weeks ago, and I was astonished by what I found. I know that Tomkins has a track record of incredibly sloppy analyses (for examples see here, here, and here), but I still wasn’t quite prepared for this article. In hindsight, I was far too naive.

In this blog post I’ll expose the many flaws with Tomkins’ data and conclusions. Be warned, this gets quite technical, but I’ll do my best to explain everything. As usual, I’ll summarise the key points at the end.

Sequence similarity

Literally the very first sentence of the “results” section of Tomkins’ article is false:

The sequence identified by Brawand, Wahli, and Kaessmann (2008) as being a vtg pseudogene is only 150 bases long and listed in the UCSC genome browser at the genome coordinates of chr1:79254632-79254781 (version hg19).

It’s really quite hard for me to understand how Tomkins could have got this impression. A mere glance at Figure 2 of Brawand et al. (2008) reveals that there are 5 sequence fragments in the human genome with matches within the chicken VTG sequence, not just 1. This isn’t a very good start, but it gets much worse. The 150 base pair (bp) fragment that Tomkins mentions is just part of one of the sequence fragments: the part that matches to chicken VTG exon 3. Brawand et al. (2008) highlight this particular fragment in their Figure S2 because they show that humans, dogs, and armadillos all possess this part of the fragment, and share 1bp deletion mutations indicating that the VTG gene was already pseudogenised in our common ancestor. Tomkins seems to have just seen this figure (Figure 1 below) and managed to conclude that this sequence is the entirety of the VTG pseudogene in these animals.

screenshot-2019-06-01-at-21.14.50.png
Figure 1 | A 150bp sequence alignment of sequences corresponding to chicken VTG exon 3. Gallus gallus = chicken, Homo sapiens = human, Canis familiaris = dog, Dasypus novemcinctus = armadillo. Figure S2 of Brawand et al., 2008 (CC-BY 4.0)

I’m actually a bit surprised Tomkins didn’t try to align the chicken VTG gene with the homologous genomic region in humans himself, as he’s very keen on doing lots of BLAST searches in his work comparing the human and chimp genomes. In his article, Tomkins mentions that BLAST isn’t really a sensitive enough tool to do the job, so maybe he’s just not very comfortable using any tools other than BLAST.

Whatever the case may be, I did decide to do the alignments myself. I extracted sequences from a range of representative mammals and non-mammalian amniotes from Ensembl representing the genomic regions either containing the VTG gene or the homologous region, depending on the species (Figure 2). As you can see in Figure 2B, only the chicken and softshell turtle have intact VTG genes in this location.

Figure 2.png
Figure 2 | Extracting homologous VTG-containing sequences. A) The genomic coordinates of the sequences from 7 species: 2 non-mammalian amniotes and 5 mammals. B) Diagram depicting the location of the extracted sequences (red dotted box) in the context of upstream and downstream genes. The armadillo sequence is cut short downstream of ADGRL4 as the scaffold is quite short, so it likely doesn’t represent the full extent of the non-coding sequence prior to the nearest gene. An alternative name for ADGRL4 is ELTD1 or LTD1.

I aligned the extracted sequences using mVISTA, a sensitive tool for global pairwise alignment of large sequences. The purpose being to see if there are any mammalian sequences in this region that match up with the sequence of the VTG gene. Figure 3 shows the results of this alignment, representing sequences that match the chicken reference as peaks.

Figure 3.png
Figure 3 | mVISTA alignment of homologous sequences showing the remnants of VTG in mammalian genomes. There are at least 5 major peaks of sequence conservation between the VTG gene sequence of chicken reference sequence and mammalian sequences (green arrows). The turtle track serves as an example of how a fully intact sequence containing functional VTG looks in a reasonably distant relative of chicken. Alignment and display were performed using the default parameters of mVISTA.

As you can see, there are 5 peaks in mammalian genomes that correspond to the chicken VTG gene sequence – this isn’t news, it perfectly matches the results of Brawand et al. (2008). These peaks represent 5 sequences that total approximately 3,100bp in length. A far cry from Tomkins’ claim that the entire alleged human VTG pseudogene is represented by just 150bp. These 3,100bp in humans have a weighted average identity of approximately 70% with the chicken VTG sequence, discounting INDELs, and 62% including INDELs. More than enough to establish clear homology, especially given the synteny (order) of these 5 sequences being conserved, in addition to the wider gene synteny of the extracted sequences.

Let’s get back to Tomkins’ version of events:

He explains that 150bp is a very small sequence compared to the intact chicken VTG gene, and in doing so, makes a very odd statement (my emphasis):

The chicken vtg1 genomic region, not including five prime promoter sequence, is 42,637 bases long (chr8:17537100-17579736) and is a large and complex gene with 35 exons. Thus, a mere 150 base fragment of alleged similarity is but only 0.35% of the actual size of a real functional vtg1 gene in chicken and hardly representative of anything approximating a real pseudogene.

He’s basically saying here that a very small fragment of a functional gene can’t be a “real pseudogene”. What makes him think that? He should know that many pseudogenes are non-functional, so over time they will be eroded by mutations. This means that depending on factors like their age, pseudogenes can be in a range of states of “preservation”. A brand-new pseudogene might look 100% like the functioning sequence except for a single premature stop codon that abolishes protein production, while a very old one could have almost disappeared entirely. Anyway, I demonstrated earlier that we don’t have just 0.35% of the VTG gene, we have at least 7.3% that is clearly identifiable, and likely more that is harder to ascertain with certainly due to being too divergent as a result of mutations. I wonder if that’s enough for Tomkins to consider it a “real pseudogene”.

Synteny

Tomkins briefly mentions the synteny of the VTG locus at the end of his article, but it seems more logical for me to cover it here. Tomkins writes:

One of the supporting arguments for the vtg pseudogene fragment being authentic is that it shares gene neighborhood synteny with chicken. When this was investigated, it was found that gene synteny surrounding the chicken vtg1 gene (~360,000 bases) compared to the region surrounding the alleged vit1 fragment in human, was completely different except for the presence of the LTD1 gene which was about three times the distance (~100,000 bases) from the alleged vit fragment as its homolog is in chicken (~38,000 bases).

This is a paragraph where I really struggle to give Tomkins the benefit of the doubt. The supporting argument regarding synteny, as he says in the first sentence, is that the VTG pseudogene is in the same gene neighbourhood as the chicken VTG gene. However, rather than evaluating the gene neighbourhood, Tomkins elects to narrowly target the ~360,000bp surrounding the VTG locus, analogous to only looking for the immediate flanking neighbours. He found just what I show in the diagram in Figure 2B: the human VTG locus is flanked on the right-hand side by 2 genes unique to placental mammals (with the possible exception of Xenarthrans): IFI44 and IFI44L, which are then followed by PTGFR. This is in contrast to the VTG gene in chicken, which is immediately flanked on the right by PTGFR. Tomkins seizes upon this as damning evidence against the synteny of the region, but neglects to look at the bigger picture of synteny in the region. Other than the addition of IFI44 and IFI44L, the 10 genes upstream and the 10 genes downstream of the VTG locus are identical between humans and chickens (Figure 4). Even outside of this block of ~20 genes with the VTG locus at the centre, there is still very good synteny for long stretches.

Figure 4.png
Figure 4 | The gene synteny surrounding the VTG locus in human and chicken. Each coloured arrow represents a gene (not to scale). Homologous genes are represented by the same colours. Direction of the arrow represents the orientation of the gene. The human VTG pseudogene is represented by the orange “X”. The red brace spans the sequence that Tomkins analysed for synteny.

It’s overwhelmingly obvious that the VTG pseudogene locus is exactly where it should be in the genome. I don’t know how Tomkins can possibly have written that paragraph with a straight face.

Functionality and the elephant in the room

Up until this point Tomkins has been attempting to cast doubt on the status of the human VTG pseudogene by refuting the commonly-cited evidence in favour of it – synteny and sequence similarity. I think I’ve demonstrated that those attempts were entirely fruitless.

However, the main content of Tomkins’ article is actually dedicated to his claim that the VTG pseudogene (just his cherry-picked 150bp part of the pseudogene, mind you), is actually a regulatory element influencing the expression of a nearby functional gene. Before I get into the technical details of this claim, let’s stop and consider a point that should be obvious, but apparently Tomkins hasn’t considered – “functional” and “pseudogene” are not contradictory descriptors! Old sequences can be adapted for new functions (“exapted”), so it’s perfectly conceivable that an old sequence fragment of a pseudogene could have been repurposed into a functional regulatory element. This is a point made by Dennis Venema in his response to Tomkins’ article in early 2016. In fact, Venema’s response covered the same topics I have up until this point in the blog post, albeit not in quite as much detail. However, he ends his article on the note that Tomkins has constructed a false dichotomy between pseudogene and function, which I don’t think goes far enough. I think we should also be critically examining Tomkins’ claim for functionality, for two main reasons.

First, because left unchallenged, Tomkins’s claim of functionality has spread far and wide, and is now essentially taken as a brute fact by many in creationist circles. ICR (via Tomkins) proclaims “It’s a functional enhancer element”, and evolutionnews (Discovery Institute) suggests that Venema’s silence on this topic is an indication that Tomkins is right:

“Venema appears to concede Tomkins’s central point that the VIT1 gene could be functional”

If you’ve ever talked to a creationist and the subject of vitellogenin came up, odds are that they cited Tomkins in their replies. To many creationists, even if the presence of function doesn’t entirely refute the evolutionary argument in favour of common descent, it at least weakens it.

Second, because it serves as an exemplar of “creation science”. Creationists are very proud of their handful of trained scientists, and often cite in-house creationist journals like Answers Research Journal to make the point “look, we’re engaging in real scientific research too!” They maintain that the only thing keeping these research articles out of “secular journals” is idealogical bias, and that the science is every bit as rigorous in creationist journals. I think it’s worth pointing out that this is clearly false.

I set out to write this blog post because I could find little online actually responding the substance of Tomkins arguments in favour of the “alleged” human pseudogene sequence being functional. It wasn’t until after I’d finished my own analysis that I found a couple of reddit posts (here and here), and the summary of a discussion on uncommondescent (here) that got into some of the details. Still, I feel like this blog post will fill a niche, and hopefully be more accessible and comprehensive than what has come before.

So without further ado, let’s get into it.

Tomkins evidence for function

Tomkins rests his case on numerous lines of evidence, from expression data to epigenomic profiling, so this is where it really gets technical. To briefly summarise, he claims that the 150bp pseudogene sequence resides in the last intron of a long non-coding “GAM” gene that is expressed in the embryonic and adult brain, and that the sequence itself is a regulatory element of this gene because it contains transcription factor-binding sequence motifs and is characterised by open-chromatin marks, a localised reduction of methylation, and is involved in long-range chromatin interactions with RNA polymerase II. All that sounds quite impressive, especially to a layperson. Just look at all that jargon that Tomkins is using – he must be right, right?

Well, no. Tomkins actually goes into very little detail in his technical research article, breezing rapidly through these supposed lines of evidence without spending any time quantifying anything or even properly describing how his figures support his claims. A layperson trying to read the article will be met with a lot of confident assertions, but have little idea where they actually come from. Let’s go through it all.

Genscan and GAM genes

The first thing Tomkins sets out to establish is that the VTG pseudogene sequence resides within a functional gene. He shows that a program called Genscan predicts a 69kb (69,000bp) gene sequence (“chr1_41.66”) that surrounds the 150bp pseudogene sequence (Figure 5). Genscan predicts genes in sequences based on a number of factors, basically looking for particular features that tend to characterise known genes. However, it has an extremely high false positive rate. As you can see in Figure 5, not only does Genscan fail to accurately predict the structure of the three known genes in the region, it predicts “gene sequences” all over the place, and few if any are likely to represent actual genes. This Genscan prediction would be weak evidence on its own, but it’s especially weak evidence of the existence of a long non-coding RNA gene, as Tomkins claims, because the program is designed to predict protein-coding genes, which have entirely different characteristics!

Figure 5.png
Figure 5 | VTG pseudogene in Genscan-predicted sequence. Known genes are shown in blue, gene sequences predicted by Genscan are shown in gold. Vertical bars represent exons, horizontal lines represent introns. The 150bp pseudogene sequence is indicated by the red arrow (also labelled “Brawand2008” in reference to the original paper). Screenshot from UCSC genome browser (hg19).

Tomkins goes on:

“Bolstering the fact that this is a functional sequence, is gene expression data that runs across the entire length of the predicted gene (fig. 2). More specifically, highly specialized gene expression experiments have functionally identified these expressed sequences in this gene as belonging to a group of novel “genomic address messenger” or “GAM” genes (Bentwich 2007).”

This is problematic in so many ways. I’ve recreated the top half of Tomkins’ fig. 2 in my Figure 6 for context.

Figure 6.png
Figure 6 | A reproduction of Tomkins’ fig. 2a. The “GAM” sequences (JD######) are shown in relation to the Genscan predicted gene sequence (gold), and the 150bp VTG pseudogene fragment (Brawand2008). Screenshot from UCSC genome browser (hg19).

In gold along the top you can see the predicted gene sequence. This is what Tomkins claims is the long non-coding RNA gene. Now he’s trying to support the claim that this gene is functional by citing the “GAM” genes from Bentwich (2007). The RNAs from these “GAM” (genomic address messenger) genes are shown below the gold predicted gene sequence, with names like “JD421963”. It’s important that I point out that those sequences haven’t been shown to be expressed from this particular predicted gene sequence: that’s not what the UCSC browser track is suggesting. Tomkins’ wording leaves this extremely ambiguous, but these sequences are those that were identified by Bentwich (2007), and this track on the genome browser is saying “those DNA sequences are also found here”. In other words, while Tomkins cites the existence of these sequences to be evidence for gene expression at this location in the genome, or at the very least heavily implies it, they don’t support that at all. Once again: those sequences have not been found to be transcribed from this “gene”.

At this point you might be wondering what these “GAM” genes are. I certainly was. The source Tomkins cites for them is Bentwich (2007), which isn’t a scientific paper, but a patent. It’s a bit beyond me, written in legalese rather than providing a clear scientific explanation, but what is apparent is that the patent includes a list of hundreds of thousands of “genes” which the patent holder claims are types of micro RNAs (miRNAs). miRNAs are very short RNAs (these ones are ~22bp long) that can influence gene expression by binding to complementary sequence in mRNAs and causing them to be broken up by enzymes. In his fig. 2 (my Figure 6) Tomkins makes it seem as though the predicted gene sequence, which he now calls a “GAM gene”, is pretty unique, possessing all these “GAM” sequences that he claims are functional. Just look how full the predicted sequence is of these “GAM” sequences, it looks quite compelling, doesn’t it?

Nope. As it turns out these “GAM” sequences are found literally all over the genome. Figures 7, 8, and 9 shows screenshots the UCSC genome browser of a few locations in the genome, and you can see how ubiquitous these sequences are. That’s what you get when you make a list of hundreds of thousands of extremely short sequences – it turns out that some of them will be present everywhere you look.

Figure 7.png
Figure 7 | “GAM” genes around the VTG locus. A zoomed out (~10x) view of the same location Figure 6, showing “GAM” sequences present everywhere. The red block highlights the sequence shown in Tomkins’ fig.2. The bottom of the screenshot cuts off many more “GAM” sequences. Known genes in blue, Genscan-predicted genes in gold. Screenshot from UCSC genome browser (hg19).
Figure 8.png
Figure 8 | “GAM” genes in a random location on chromosome 6. This screenshot is on the same scale as Figure 6, but shows “GAM” sequences present in a random location on a chromosome I chose at random, to illustrate that these sequences really are everywhere. Genscan-predicted genes in gold. Screenshot from UCSC genome browser (hg19).
Figure 9.png
Figure 9 | Showing a larger extent of “GAM” sequences from another random location on chromosome 6. This screenshot is zoomed out to a similar extent to Figure 7, but this time captures more of the extent of “GAM” sequences that are ubiquitous throughout the genome. Known genes in blue, Genscan-predicted genes in gold. Screenshot from UCSC genome browser (hg19).

As these “GAM” sequences are present absolutely everywhere, finding them in a particular location does nothing to support the claim of functionality. Far from demonstrating the expression of functional sequences in this predicted gene, all Tomkins has done is demonstrate that some ubiquitous sequence motifs are present in the vicinity – an extremely mundane point. Once again, it seems impossible to imagine that Tomkins is unaware that these sequences can be found everywhere in the genome, since all it would have taken for him to realise it would be to click a single “zoom out” button from the view where he took the screenshot for his figure 2. The only conclusion is that Tomkins knows these sequences are ubiquitous, but withheld this information from his “technical paper” in service of his narrative.

RNA expression data

What about the second half of Tomkins’ figure 2? This is where he finally discusses some genuine expression data from RNA sequencing (RNA-seq) and microarray experiments. Figure 10 shows a reconstruction of his entire figure 2 for context. We already covered Figure 10A in the last section, now we’re looking at 10B.

Figure 10.png
Figure 10 | A modified reproduction of Tomkins’ figure 2. A) see Figure 6. B) the yellow track represents raw smart-tagged RNA-seq signal, and the green track represents raw RNA-seq signal, both from Maunakea et al. (2010). The green track was absent from Tomkin’s figure. The blue track represents microarray-based human brain expression data generated by the Sestan lab. Screenshot from UCSC genome browser (hg19).

SMART-tagged RNA seq

Tomkins dedicates literally a single sentence to describing the RNA-seq data:

In addition to the GAM gene study, another research project identified expressed sequences from this gene associated with the human brain (the frontal cortex region) based on RNA extracted from a 57 year old male (Maunakea et al. 2010).

That’s it. Notice how he doesn’t even attempt to quantify the amount of expression? There are some sequences from this gene expressed, and apparently that’s good enough for Tomkins, he’s ready to move right along. Well, unfortunately, simply saying the word “expressed” isn’t good enough. Look at Figure 10B – what do those 3 blue bars in the yellow track mean? What’s “RNA-seq Smart-Tagged Raw Signal”? Tomkins doesn’t explain, so I’ll have to do his job for him.

RNA-sequencing is like it sounds – RNA produced in cells is sequenced, and then mapped onto the genome to show where they came from. So-called “SMART-tagged” RNA-seq is a modification of this method, which basically tags the starting end of the string of transcribed RNA with a recognisable tag that can be selected for later. Long story short, sequencing reads with this special tag represent the 5′ (starting) end of the RNA transcript, so can be used to infer transcriptional start sites (TSSs). Of course, nothing in biology is perfectly clean-cut, so these reads can also be found outside of TSSs. In their methods section, Maunakea et al. (2010) explain how they used the tagged and non-tagged RNA-seq reads to infer TSSs:

We tallied number of SMART and NOSMART RNA-seq reads overlapping with each island, and defined TSS activity as (1) having at least five SMART-tagged reads, and (2) at least 70% of total RNA-seq reads are SMART tagged reads.

So, their very first criteria for being a TSS is that the site must be represented by at least 5 SMART-tagged sequencing reads. How many reads do those 3 blue bars in Figure 10B (yellow track) represent? Just one. Now look at the green track, representing the raw signal of the standard RNA-seq experiments, and be aware that this track is actually conspicuously absent from Tomkins’ figure. It shows absolutely no expression data for the entire “GAM gene” whatsoever. None of the “GAM” sequences are represented by any sequencing reads – even the 3 separate SMART-tagged reads don’t overlap with any of the “GAM” sequences. How Tomkins thinks this is supposed to support his claim of these sequences being expressed and functional is beyond me.

I’m not done yet. Figure 11 displays another zoomed-out view of the region, so you can see how the data looks for genuine genes. The red box highlights Tomkins’ claimed “GAM gene”, with the 3 blue bars from Figure 10B that Tomkins seems to think is indicative of functional gene expression barely visible.

Figure 11.png
Figure 11 | Zoomed out view of VTG region showing RNA-seq data. The red box highlights the 3 blue columns in Figure 10B (yellow track). Compare the peaks, or lack thereof in the region within Tomkins’ “GAM gene” (chr1_41.66) to the peaks within the 2 real flanking genes: IFI44 and ADGRL4. Screenshot from UCSC genome browser (hg19).

This data is from a single study, the lack of expression data doesn’t prove that this sequence is non-functional, but it certainly doesn’t provide positive evidence in favour of function either. I really don’t know how Tomkins could have thought this complete lack of expression data would be compelling positive evidence of expression.

Sestan microarray data

The second half of Tomkins expression evidence is based on microarrays. Briefly, microarrays are glass chips with lots of short DNA sequences (probes) stuck to their surface. These probes are designed in such a way to be complementary to transcripts of known genes. RNA from cells is isolated, converted into DNA, and washed over the chip. DNA fragments with complementary sequences to the probes on the chip stick to them, and can then be detected using fluorescence. The more fluorescence, the more of that transcript is present, so this is a fairly primitive but effective way of quantifying relative gene expression. More on that in a moment, but first here’s what Tomkins says:

“In confirmation of the RNA sequencing studies for this GAM gene in human brain tissue, the HBT (Human Brain Transcriptome) project at the Department of Neurobiology Yale University School of Medicine also deposited data in the UCSC genome browser from extensive microarray studies showing that this particular GAM gene is not only expressed in the brain tissues of mature humans, but also in the developing fetus (Kang et al. 2011; Pletikos et al. 2014).”

Once again, Tomkins doesn’t attempt to make any quantification of this “expression”, but unlike last time, he’s partly justified in doing so. Unless specifically designed in a particular way, microarray experiments are only capable of quantifying relative gene expression levels. In other words, they take some kind of control microarray, then compare all the fluorescence levels from the tissue sample microarrays to this control, and look for increased or decreased gene expression in the samples as compared to the control.

The relevant data in microarray track (Figure 10B – blue track) is actually the colour of the bars. The bars themselves just represent the sequences covered by the probes on the microarray chip, and then the colour represents the increased or decreased fluorescence signal representing bound transcripts. The typical colour scheme is that red represents increased expression and green represents decreased expression. The more red the bar is, the more signal was found in the tissue sample relative to the control, and the more green, the less signal relative to the control. It’s a little hard to see in Figure 10B (and in Tomkins’ fig.2), but all of the 13 sets of bars (representing results from 13 brain tissue samples) are deep black. In other words, they detected no significant differences between the expression level of those sequences in those 13 brain tissues to the expression level in the control. Just as an aside, in this case, the control wasn’t a physical chip, but an “idealised reference chip”, which basically represents the mean values of all the probes on the chip.

So, the impressive-looking bars Tomkins shows are completely meaningless on their own, and the fact that they are all black just means that no significant differences in expression could be found between the 13 brain tissues samples. There is literally zero information about RNA expression levels here, because the sequences could either be highly expressed in all 13 tissues, or not expressed at all in all 13 tissues, or anywhere in-between. The track Tomkins shows doesn’t support his claim that this “GAM gene” is expressed in the brain tissues of the developing fetus in the slightest.

Yet again, I find myself struggling to believe that Tomkins is unaware of all this. Can we really expect an experienced geneticist to be completely unaware of how to interpret microarray expression data, something that is taught to first-year undergraduates (at least these days)? Can we expect Tomkins to have not done his due diligence when doing his background reading to interpret this data? Is it really plausible that he saw the black bars in the UCSC genome browser track and immediately concluded that they represented expressed sequences, or is it more likely that he just wants his creationist audience to do that instead?

Transcription factor binding sites and epigenomics

Just for context, I’ve reproduced Tomkins’ figure 3 in its entirety in Figure 12 below. I’ll go through each part of it one-by-one, showing the relevant parts of the figure again as necessary.

Figure 12.png
Figure 12 | A reproduction of Tomkins’ figure 3. A) The 150bp VTG pseudogene fragment identified by Brawand et al., 2008 and a MAFK binding site, B) Conserved transcription factor binding sites between human, mouse, rate (HMR), C) Methylation and open chromatin tracks, D) ChIA-PET track representing long-range chromatin interactions in association with RNA Polymerase II.

Transcription factor binding sites

Tomkins says:

“The first line of evidence that the alleged vtg fragment contains functional features is that five different highly conserved mammalian transcription factor motifs are identified in the middle of the 150 base sequence for the following factors: FREAC3, FREAC2, FOX04, FOX FOX03 [sic], FOX01, and SRY”

He lists 6 different transcription factors (despite saying “five different”) and I think that most lay readers will get the impression that this is 6x better evidence than finding a single motif, for example. There’s 6 different independent indications of function! But what you may notice from Tomkins’ own list there, is that 3 of the transcription factors (TFs) have very similar names: FOX04, FOX03, and FOX01. First, a minor quibble: Tomkins names these transcription factors FOX zero 4, zero 3, and zero 1, but their actual names, even as listed on Tomkins’ source, is FOX O (the letter) 4, etc. The fact that he made this mistake leads me to believe that he didn’t do any background reading about any of these transcription factors that are shown in the UCSC genome browser track.

If he had done, he might have noticed that in addition to those 3 named FOX TFs, FREAC3 and FREAC2 are in fact also FOX-family TFs! Their more common names are FOXC1 and FOXF2, respectively. So, we have 5 FOX-family TF motifs out of a grand total of 6. This shouldn’t be particularly surprising, given that Tomkins’ own figure (my Figure 12) shows all 6 motifs being present in a single location, all overlapping. It also isn’t surprising to learn that all 6 TF-binding motifs are highly similar, as you can see in Figure 13.

Figure 13.png
Figure 13 | Binding motifs of the six transcription factors matching a single motif in the human VTG pseudogene fragment. The height of each letter corresponds to the enrichment of this base in the given position across known binding motifs for the TFObtained from the Jaspar Core 2018 database. Reverse complemented.

Simply put, where you find a binding motif for a particular FOX TF, you are bound to also find motifs matching a number of other FOX TF and probably the SRY TF as well, because their motifs are all very similar. For this reason, Tomkins’ evidence here boils down to finding a single motif of approximately 8bp (TTGTTTAC) that is well conserved (between humans and rodents).

I realise it may seem extremely pedantic of me to spent so much time making the point that Tomkins implication of multiple binding sites supplying multiple independent lines of evidence for function is wrong, but perhaps this next quote, from a few sentences later in the article, provides some justification:

“Not only is this alleged vtg regio n [sic] replete with transcription factor binding sites”

He actually says the sequence is “replete” with binding sites, when in fact his data indicates a single binding site capable of binding several transcription factors (not at the same time!) with similar motifs, which is an important distinction. Functional elements are often represented by clusters of conserved motifs, rarely by single motifs.

Such single conserved motifs are ubiquitous throughout the genome, as you can see in Figure 14. Some of these represent real, functional binding sites in promoters and other regulatory elements like enhancers, but many do not.

Figure 14.png
Figure 14 | Conserved TF-binding motifs are ubiquitous in the genome. A zoomed-out view of the VTG region with the 6 TFs highlighted by Tomkins in the VTG pseudogene fragment indicated by the red brace.

This short motif isn’t unique to humans and rodents, either. Remember this motif is present in the 150bp of sequence that aligned to part of the 3rd exon of the functional chicken VTG gene. As you can see in Figure 15, the human motif (highlighted) aligns very well to the chicken sequence, the main difference being a single 1bp insertion or deletion. In other words, only a single base needed to change in the common ancestor of mammals that lost the function of VTG in order to create this motif that Tomkins is claiming is a functional transcription factor binding site. So, even if Tomkins is right about all his claims of functionality, this would actually be further evidence that new functions are quite readily evolvable.

Figure 15.png
Figure 15 | The FOX transcription factor binding site in the VTG pseudogene fragment. Highlighted in green is the TFBS in humans and corresponding sequences in 3 other species, including chicken. Gallus gallus = chicken, Homo sapiens = human, Canis familiaris = dog, Dasypus novemcinctus = armadillo. Modified from Figure S2 of Brawand et al., 2008 (CC-BY 4.0)

A second piece of TFBS evidence Tomkins mentions is the presence of a MAFK transcription factor binding site that overlaps the VTG fragment. This time it represents an actual, tangible, binding event, since the data comes from ENCODE ChIP-seq assays, meaning that they actually found the MAFK protein in physical contact with this sequence of DNA. In particular, they found MAFK proteins bound in an area, with the binding peak being about 100bp upstream of the VTG fragment in 2 cell lines: HEPG2 (liver cancer cell) and IMR90 (lung fibroblast cell). How strong was the signal? Not very. About 35% of the mean cluster score plus one standard deviation, if you want to be specific. The same set of ENCODE ChIP-seq assays looked for sites bound with any of 161 different transcription factors (many of them general factors, such as EP300 and CTCF), in dozens of different cell lines. The fact that all of that produced a single mediocre bound TF overlapping the VTG fragment is practically evidence against function, rather than for it. Such single and/or weak bound TF signals are a dime a dozen throughout the genome, and I doubt that even Tomkins believes they’re all functional.

A “Dip” in DNA Methylation

Tomkins claims that there is a reduction of methylation levels in the 150bp VTG fragment, and that this is indicative of function:

“In fact, the entire 150 base sequence of the alleged vtg fragment is characterized by a central dip (lowering) in cytosine methylation directly within the predicted binding domains mentioned above”

It’s true that active enhancers display reduced DNA methylation, but is Tomkins accurately describing the data here? In his figure 3C, he displays the “MeDIP-seq Raw signal” track,  and sure enough, there seems to be a distinct “dip” in the track where the fragment is (Figure 12C).

What he doesn’t display is another track that is part of the same track set called “MeDIP CpG Score”. Tomkins would have had to manually turn this track off to make his figure 3, so he’s perfectly aware of what I’m about to say.

The “Score” track is the one that reports significant levels of DNA methylation, rather than displaying the raw signal. There are no significant peaks in the raw signal near the 150bp VTG fragment. The sequences flanking the “dip” have a raw signal value of 3-4. The average raw signal value of the genome as a whole is about 5. In other words, the sequences within and around the VTG region don’t stand out against the background noise (Figure 16). If the flanking regions don’t have a significant level of methylation, then the sequence between them can’t have a significant lowering of methylation that might be indicative of an active enhancer.

There’s always going to be peaks and troughs in background noise, and the vast majority aren’t biologically meaningful. Tomkins should have just said “the fragment isn’t strongly methylated” – it would have been more accurate, and technically still left the door open to his enhancer idea. A lack of strong methylation is consistent with active enhancers, but it’s also consistent with a large proportion of non-functional sequences too (Figure 16).

Figure 16.png
Figure 16 | The “dip” in DNA methylation. A zoomed-out reproduction of part of Tomkins’ figure 3. The red arrow indicates the “dip” in methylation relative to the methylation levels immediately left and right of the arrow.

Open Chromatin – FAIRE and DNase1

The next signature of function Tomkins tries to rally is open chromatin, in this case as measured by FAIRE and DNaseI hypersensitive sites. Active regulatory elements must reside in open chromatin in order to be bound by activator proteins. FAIRE-seq is a technique used to detect regions of the genome that are unoccupied by nucleosomes (open), and DNase-seq searches for sites in the genome that are sensitive to being cleaved by the DNaseI enzyme, indicating that the site is “open” to the enzyme, and therefore in an open conformation, free of nucleosomes. These 2 techniques are used with the same goal in mind, but have different pros and cons that aren’t relevant to the current discussion.

As usual, Tomkins barely even bothers to describe the results of these profiling efforts, this time with a mere sentence fragment:

“…combined with open-active chromatin profiles associated with nucleosome depletion…”

So, he’s claiming that the 150bp VTG pseudogene fragment bears signatures of open chromatin, and therefore of functionality, again without giving any kind of quantification. The only glimpse we get is in his figure 3C, the relevant part of which I have reproduced in Figure 17.

Figure 17.png
Figure 17 | Reproduction of Tomkins’ figure 3C depicting open-chromatin profiling data. A) FAIRE-seq tracks, B) DNase-seq tracks.

In Figure 17A you see the FAIRE profiles – 6 tracks representing “DS” and “OS” tracks for 3 different human cell lines: GM12678, H1-hESC, and K562 (lymphocyte, embryonic stem cell, and myelogenous leukemia cell, if you were wondering). DS and OS refer to “Density Signal” and “Overlap Signal”, but the subtle distinction between them isn’t important to us here. The results of the DNase-seq in 3 different cell human cell lines are shown in Figure 17B. K562 again, HUVEC (umbilical vein endothelial cell), and HepG2 (liver cancer cell). The tracks are named “Sig”, referring to per-base signal.

What do all the grey boxes mean? Tomkins doesn’t explain, and I don’t think I’m stepping out of line to suggest that most of the readers of his article don’t know. Tomkins seems to just present these shaded grey boxes as if to say “look – there’s sciencey-looking stuff in the diagram, so trust me when I say that it means the chromatin is open”.

The darker the bar, the more signal is present in the sequence covered by the bar. The amount of FAIRE and DNaseI signal there is actually no higher than the background level of the rest of the genome. Nothing about the grey bars indicates some kind of “peak” of signal that might be associated with an active regulatory element, as Tomkins tries to imply. I feel like I’m becoming a broken record at this point, but the worst of it is that Tomkins knows the signal he’s showing in his figure isn’t significant. I know he knows, because in order to produce his figure 3C (my Figure 17), he had to deliberately turn off or move out of sight a set of tracks that are associated with the FAIRE and DNaseI profiling tracks by default. As with the “MeDIP CpG Score” track in the previous section, the open chromatin profile tracks also have “peak” and “hotspot” tracks that show where the statistically significant levels of signal are present above the background, but Tomkins decided not to show them in favour of the “signal” tracks which superficially made it look like something was there. He didn’t show these tracks because they are completely empty – no significant FAIRE or DNaseI signal was detected in any of the 5 cell lines that Tomkins chose (Figure 18).

Figure 18.png
Figure 18 | No significant open chromatin signals. A reproduction of Tomkins’ figure 3C, this time including the empty “peak’ (Pk) and “hotspot” (Hot) tracks that Tomkins’ omitted (highlighted in yellow). A) FAIRE-seq tracks, B) DNase-seq tracks.

In other words, the chromatin in and around the 150bp fragment that Tomkins is claiming is an active enhancer is no more open (active) than the background in these datasets. To use this data to try and suggest functionality just doesn’t make any sense at all.

Long-range chromatin interactions with RNA polymerase II

The final signature of function Tomkins mentions is long-range chromatin interactions of the 150bp VTG fragment associated with the enzyme RNA polymerase II, which synthesises RNA from DNA in the process of transcription:

“chromatin interaction analysis of paired-end tags (ChIA-PET) tracks showing that this region is involved in long range chromatin interactions associated with RNA polymerase2 and gene transcription”

Without going into the details of ChIA-PET, it’s a technique designed to detect chromatin interactions (which sequence of the DNA is associated with another sequence in 3D space), and simultaneously see which proteins are present at this interaction. When a regulatory DNA sequence is active, it has proteins bound to it that cause it to loop around to the promotor of the gene it’s regulating and influence transcription. For this reason, when 2 DNA sequences are found by ChIA-PET to be in contact with one another in the presence of the primary transcription enzyme, it can be an indicator that one of those 2 DNA sequences is a regulatory element. This is what Tomkins is claiming is shown in his figure 3D. I’ve reproduced this part of his figure in Figure 19, and as you can see, it’s more grey bars.

Figure 19.png
Figure 19 | A reproduction of Tomkins’ figure 3D.

Exactly as before, the darker the grey in the bar, the more signal is present, in this case the data is from 2 replicates of another human cell line: MCF7 (breast cancer cell). Exactly as before, Tomkins doesn’t attempt to explain the meaning of the grey bars to his readers, he seems to hope that they’ll see a couple of shades of grey and be satisfied that it means something profound.

In the UCSC genome browser, when you turn on this ChIA-PET track for this cell line, the default tracks shown aren’t Sig 1 (signal from replicate 1) and Sig 2, they’re Sig 1 and Int 1. Int 1 is the track that, and stop me if you’ve heard this one before, reports the locations of interacting DNA sites determined by possessing a significant signal in replicate 1. Tomkins decided to turn the Int 1 track off, but Sig 2 on. He also decided not to turn on Sigs 3 and 4, never mind Ints 2-4. Figure 20 shows how his figure would have looked if he had decided to turn all of these tracks on.

Figure 20.png
Figure 20 | Full representation of ChIA-PET tracks for the VTG fragment. A) MCF7 Pol2 Signal tracks for replicates 1-4, B) MCF7 Pol2 interactions tracks for replicates 1-4. The 6 tracks highlighted in yellow were not shown in Tomkins figure 3D.

It’s easy to see why Tomkins didn’t include them in his figure – all 6 tracks are completely empty! By deliberately not showing the Int tracks, Tomkins hid the fact that the signal shown by the grey boxes is insignificant relative to the rest of the genome, and by not showing the Sig 3 and 4 tracks, he hid the fact that any signal at all was only found in 2 out of 4 replicates studied.

Just to give you an idea what these tracks look like on a larger scale, Figure 21 shows a view that is zoomed out 10,000x compared to Figure 20, covering 2.25Mb. At this scale, we can see what “real” peaks in signal look like, and that these peaks are represented in the Int tracks. In fact, we can see that Sigs 3 and 4 seem to have a much better signal:noise ratio than Sigs 1 and 2 in this region – they still recover strong peaks (highlighted in yellow) while having much less background signal outside of these peaks. Highlighted in red you can see the VTG fragment that Tomkins is focused on, and that the “signal” Tomkins trumpets from Sig 1 and 2 is really just part of the background noise.

Figure 21.png
Figure 21 | ChIA-PET data for 2.25Mb around the VTG fragment. Some peaks in signal are highlighted in yellow, while the VTG fragment is indicated in red.

If I was feeling like a broken record before, surely now I must be shattered. This is yet another example of Tomkins misrepresenting the significance of his data in service of his narrative of functionality. There is absolutely no way to interpret the ChIA-PET data he shows as any kind of good indicator of function. If there had been peaks in the VTG fragment similar to those highlighted in yellow in Figure 21, then maybe, but he’s literally just showing us background noise.

Summary

We finally made it to the end. Or you skipped ahead to the summary because it all got a bit too technical and repetitive, and to be honest, I don’t blame you. I didn’t expect this blog post to get this long and detailed when I first began work on it, but the more I looked into Tomkins’ data, the more the full extent of the problems became clear to me and I decided it all needed to be spelled out in excruciating detail.

First, I showed that Tomkins’ paper is based on a false premise – that a 150bp sequence is the full detectable extent of the VTG pseudogene in humans. In fact, it’s at least on the order of 3000bp in total. These sequence fragments reside precisely where they would be expected to be if they really were remnants of the ancestral VTG gene: between the ADGRL4 and PTGFR genes in quite a large region of synteny with the Chicken genome.

Next, I addressed a second false premise in Tomkins’ paper – that a sequence doesn’t represent a pseudogene if it can be shown to currently be functional. To the contrary, it’s perfectly possible for a functional element to evolve from a pseudogene sequence, and in many cases still be clearly recognisable as such. There’s no a priori reason for anyone to desperately want the sequence fragment to be non-functional to preserve a piece of evidence in favour of evolution.

Despite this, I decided to evaluate Tomkins’ claimed evidence in favour of functionality, and I really was blown away by what I found. Broadly speaking, Tomkins described seven lines of evidence, and here’s how they turned out:

  1. Genescan predictions – These are ubiquitous throughout the genome.
  2. “GAM gene” sequences – These are ubiquitous throughout the genome.
  3. Gene expression – No significant signal.
  4. TFBSs – Single 8bp motif conserved, and ubiquitous throughout the genome.
  5. Hypomethylation – No significant signal/ubiquitous.
  6. Open Chromatin – No significant signal.
  7. Chromatin interactions + RNA Pol II – No significant signal.

You might think that the strength of each individual piece of evidence is weak on its own, but that together it might add up to good evidence. Certainly not in this case. Nada multiplied by seven really does equal diddly squat.

As an exemplar of “creation science”, Tomkins’ article is a perfect one, because it is just so bad. No scientific journal worth its salt would accept such a flawed analysis, which is why Tomkins would’ve had no choice but to publish it in Answers Research Journal. Tomkins is presumably capable of writing genuine scientific articles – consider his real publication record – but for some reason he throws all his rigour out the window when he puts his “creation scientist” hat on. We can’t be sure whether or not Tomkins does this consciously, but as I mentioned several times in the article, there’s a very good case to be made that he knows exactly what he’s doing, which is carefully misrepresenting scientific data to fit his idealogical biases.

I’ll leave you with one final gaffe in the article that caught my eye, from the end of Tomkins’ conclusion. Pay attention to the bolded text:

“However, the real story is that the alleged 150 base vtg sequence is not a pseudogene remnant at all, but a functional enhancer element in the fifth intron of a “genomic address messenger” (GAM) gene. This particular GAM gene produces long noncoding RNAs that have been experimentally shown to selectively inhibit the translation of known target genes, a majority of which have been implicated in a variety of human diseases. Messenger RNAs from this particular gene are also known to be expressed in a variety of human brain tissues in both fetal and mature subjects in three separate studies.”

Messenger RNAs are RNAs transcribed from protein-coding genes that go on to be translated by ribosomes. It’s quite a specific technical term. Tomkins seems to have forgotten, not for the first time, that he’s claiming that the so-called “GAM gene” is a long non-coding (lnc) RNA gene, not a protein-coding gene!

 

Comments and queries are welcome.

-RM

4 thoughts on “Tomkins on the Human Vitellogenin Pseudogene: Who Needs Signal When You Have Noise?

  1. [Disclosure: I skipped the detail] It’s significant that Tomkins explicitly addresses Biologos, rather than the purely scientific literature. It is clear that his main enemy is the Christian who has the temerity to accept the evidence for evolution. One can only guess what Tompkins’ opinion is of the Pope.

    Like

  2. Excellent work, many thanks for diving into the nitty gritty. An anecdote related to the FOX binding site thing – a few years ago, after the infamous ENCODE papers came out, a colleague in the math department made me some artificial chromosomes of 10 million bps in length in Mat lab, some just purely random sequences, some premised on base frequencies, etc. I then ran these through some online software designed to locate tf binding sites – all the ones known at that time. When I got the results, I was a bit shocked – I had more tf binding sites than I did sequence (due to overlap, negative strand stuff, etc.). Point is, if known tf binding sites show up literally all over the place in random assemblages of bases, the fact that they show up in unexpected places in real genomic data does not mean what it might seem to, at least to me. These are, after all, pretty short stretches of bases, shouldn’t be too surprised that they show up, even just by chance, all over the place.

    Liked by 1 person

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s