Discovery of a New TLR Gene and Gene Expansion Event through Improved Desert Tortoise Genome Assembly with Chromosome-Scale Scaffolds.

Toll-like receptors (TLRs) are a complex family of innate immune genes that are well characterized in mammals and birds but less well understood in nonavian sauropsids (reptiles). The advent of highly contiguous draft genomes of nonmodel organisms enables study of such gene families through analysis of synteny and sequence identity. Here, we analyze TLR genes from the genomes of 22 tetrapod species. Findings reveal a TLR8 gene expansion in crocodilians and turtles (TLR8B), and a second duplication (TLR8C) specifically within turtles, followed by pseudogenization of that gene in the nonfreshwater species (desert tortoise and green sea turtle). Additionally, the Mojave desert tortoise (Gopherus agassizii) has a stop codon in TLR8B (TLR8-1) that is polymorphic among conspecifics. Revised orthology further reveals a new TLR homolog, TLR21-like, which is exclusive to lizards, snakes, turtles, and crocodilians. These analyses were made possible by a new draft genome assembly of the desert tortoise (gopAga2.0), which used chromatin-based assembly to yield draft chromosomal scaffolds (L50 = 26 scaffolds, N50 = 28.36 Mb, longest scaffold = 107 Mb) and an enhanced de novo genome annotation with 25,469 genes. Our three-step approach to orthology curation and comparative analysis of TLR genes shows what new insights are possible using genome assemblies with chromosome-scale scaffolds that permit integration of synteny conservation data.

Introduction

Expansion of gene families through duplication creates the opportunity for the evolution of new functionalization, specialized subfunctionalization, or differences in gene expression (Ranz and Parsch 2012). Conversely, loss within a gene family may reflect reduced selective pressure to retain gene function sometime during a lineage’s evolutionary history (Albalat and Cañestro 2016). Within vertebrates, the immune response is a fundamental part of how an animal interacts with its environment. Toll-like receptors (TLRs) are a large family of signaling transmembrane proteins that trigger an innate immune response following pathogen exposure (Akira and Takeda 2004). Although somatic recombination produces genetic diversity in the adaptive immune system, genes in the innate immune system (e.g., TLRs) must evolve rapidly in other ways to keep pace with the evolution of pathogens. Although TLRs are well studied in endothermic mammals and birds (Hopkins and Sriskandan 2005; Velová et al. 2018), there is less known about TLR gene family evolution in nonavian sauropsids (nonavian reptiles; but see Kahn et al. 2019; Liu et al. 2019). The innate immune system of nonavian sauropsids remains poorly characterized overall (Zimmerman et al. 2010) due to the lack of experimental tools for functional immunology and seasonal changes in ectotherm physiology (Zapata et al. 1992; Sandmeier et al. 2016). Studying TLR evolution of nonavian sauropsids can provide insight into the adaptation of the innate immune response in this group, provide a broader amniote context to TLR evolution in birds and mammals, and shed light on gene evolution broadly.

Genome resources yield information about the gene repertoire of nonmodel species through comparative evolutionary analyses. These analyses depend in part on the contiguity of reference genomes. Here, we analyze two TLR subfamilies across 22 tetrapod genomes: the TLR7 subfamily (TLR7, TLR8, TLR9) that recognizes nucleotide-based pathogen-associated molecular patterns, and the more diverse TLR11 subfamily (TLR11, TLR12, TLR13, TLR21) that recognizes pathogen-associated molecular patterns in bacterial proteins, bacterial RNA, and CpG DNA viruses (Heil 2004; Keestra et al. 2010). To aid these analyses, we used new sequencing and scaffolding technology to yield chromosome-scale scaffolds and a de novo annotation for the threatened Mojave desert tortoise (Gopherus agassizii; genome build: gopAga2.0). We then employed a three-step integrative approach to manually curate TLR orthology assignments across 22 tetrapods. This study demonstrates that much can be learned about nonmodel organisms and well-studied gene families through a comparative genomic approach, which is broadly applicable to studies of other complex gene families.

Results

The final assembly length for the gopAga2.0 assembly is 2.34 Gb, with a longest scaffold length of 106.5 Mb, N50 of 28.36 Mb, and L50 of 26 scaffolds (table 1). The number of predicted haploid chromosomes is 26 (Olmo 1981), suggesting the new assembly contains chromosome-scale scaffolds. The final assembly length for gopAga2.0 is 2.34 Gb, with a longest scaffold length of 106.5 Mb, N50 of 28.36 Mb, and L50 of 26 scaffolds (table 1). The number of predicted haploid chromosomes is 26 (Olmo 1981), suggesting the new assembly contains chromosome-scale scaffolds. Bioinformatic chromosome painting reveals strong synteny conservation between chicken and desert tortoise (fig. 1B, top). As expected for more evolutionarily distant species, there is a decreased synteny conservation between desert tortoise and green anole, as fewer linkage blocks are identified and there is greater mixing of synteny blocks (mixed colors in fig. 1B, bottom).

Table 1

Assembly and Gene Statistics across Genome Builds

Assembly Statistics
	Total Length (Mb)	L50 (No. Scaffolds)	N50 (Mb)	Longest Scaffold (Mb)
gopAga1.0	2,399	2,592	0.25	2
gopAga2.0	2,338	26	28.36	107
	Assembly BUSCO Statistics			Annotation Statistics
	Complete	Fragmented	Missing	Total Genes	Mean Gene Length (bp)
gopAga1.0	92.7%	4.6%	2.7%	20,172	19,062
gopAga2.0	93.7%	3.7%	2.6%	25,469	21,364

note.—Assembly BUSCO statistics are from the Tetrapoda gene set in the genome assemblies. Total genes and mean gene length are calculated based on the annotations for gopAga1.0 and gopAga2.0 (see supplementary table S2, Supplementary Material online).

Regarding TLR evolution, we found that the TLR7 and TLR8 genes are tandemly arranged in a highly conserved syntenic region across tetrapod genomes (fig. 2A) and exhibit a large amount of gain/loss. The three-step orthology approach also revealed a new, misidentified TLR homolog exclusive to nonavian sauropsids, which we refer to as TLR21-like (fig. 3). For more details, see the Discussion section.

. 3.

—TLR21-like gene discovered. (A) Syntenically conserved regions of TLR21-like (gold, middle) for 22 tetrapods. TLR21-like does not colocalize near TLR21 within the genome. Anolis carolinensis has LINGO1 (a leucine-rich repeat gene) instead of TLR21-like (yellow outline). Synteny information was not available for Celonoides abingdonii; therefore, only the present/absence of genes are shown (i.e., no black line). (B) Phylogenetic reconstruction of TLR11 gene subfamily including the TLR21-like gene discovered here (gold highlight). Posterior support values for nodes are 1 unless labeled. (C) Cartoon showing how TLR21 protein initiates the innate immune response by recognizing the CpG nucleotide sequence.

Discussion

TLR8 Gene Expansion

All species surveyed have one copy of TLR7, whereas the TLR8 gene varies among tetrapods. Our analyses confirm the duplication of TLR8 in turtles (Liu et al. 2019) and further reveals that TLR8 duplication (TLR8B) is also present in crocodilians. Mammals have only one copy of TLR8, and squamate reptiles and birds lack a TLR8 ortholog (fig. 2). In addition, we find that the third copy of TLR8 discovered in Chinese softshell turtle (Pelodiscus sinensis; Liu et al. 2019) is actually part of a second duplication event specific to testudines (TLR8C). The TLR8 protein recognizes single-stranded RNA viruses, suggesting this could produce subspecialization or expression differences of this function in crocodilians and especially turtles (fig. 2).

. 2.

—TLR8 gene expansion in turtles. (A) Syntenically conserved region for TLR7 and TLR8 (gold, middle) for 21 tetrapods (Chelonoidis abingdonii is excluded because its TLR8 orthology cannot be determined from available data). Genes are labeled at the top and cladogram is shown on left. Black horizontal lines represent contiguous sequence. TLR8 is expanded once in testudines and crocodilian archosaurs (TLR8B) and a second time in testudines (TLR8C) followed by pseudogenization (gray genes). (B) Phylogenetic reconstruction of TLR7 gene subfamily; reptiles are bolded, gene expansions are highlighted in gold. Posterior support values for nodes are 1 unless labeled. (C) Cartoon summarizing how TLR7 and TLR8 proteins initiate the innate immune response by recognizing single-stranded RNA (ssRNA) viruses. (D) A summary of inferred gene duplication based on results in this article.

Within chelonians, the Pinta Island tortoise has two predicted TLR8 genes, but one consists of only 122 amino acids and contains no protein motifs. Gopherus agassizii has three TLR8 genes, one with a conserved open reading frame (TLR8-1/TLR8B), one that is pseudogenized based on accumulation of stop codons (TLR8-2/TLR8C, also pseudogenized in Chelonia mydas), and TLR8-3/TLR8A that is truncated to include only the N-terminal half but retains nine leucine-rich repeats (LRRs; the basis for TLR specificity and function), suggesting protein function is retained (supplementary fig. S2, Supplementary Material online). Predicted truncation of TLR8-3 occurs at the site of a derived TGA stop codon that—surprisingly—is polymorphic among individuals within G. agassizii (CGA/TGA; supplementary fig. S9, Supplementary Material online). The stop codon in the gopAga2.0 type specimen was validated by mapping available RNA-seq and genomic resequencing data of unrelated individuals (see Supplementary Material online). The CGA to TGA mutation can readily occur through deamination of cytosine to uracil and correction of uracil to thymine. A truncated TLR paralog has been observed in several fish species where there is both membrane-bound TLR5M gene and soluble TLR5S gene that lacks a TIR and transmembrane domain but retains 21 LRRs (reviewed in Rebl et al. 2010). These proteins were shown not only to exhibit different functions but also different tissue specificities. Tissue-specific expression patterns were also observed among the TLR8 paralogs in the Chinese softshell turtle (Liu et al. 2019).

Using the above results, we outline a new naming scheme based on the inferred evolutionary duplication history of TLR8. We suggest that TLR8A (TLR8-2, TLR8-3 sensu Liu et al. 2019) originated from TLR7 first because it is the paralog most shared across taxa (fig. 2). TLR8B then duplicated from TLR8A within Archelosauria (testudines, crocodilians, birds) and was secondarily lost in birds. TLR8C then duplicated from TLR8B specifically within testudines, leading to pseudogenization within the two nonfreshwater species (green sea turtle and Mojave desert tortoise). The LRR motifs for TLR8-1 and TLR8-2 in the outgroup (western clawed frog) are nearly identical (supplementary fig. S7, Supplementary Material online) and differ from the other TLR8 motifs (supplementary fig. S2, Supplementary Material online), suggesting this duplication was separate. The duplication of TLR8 genes in nonavian sauropsids and the truncation in tortoises may reflect diversification of the innate immune response that warrants further study.

New TLR21-like Gene in Nonavian Sauropsids (Reptiles)

Many species have two genes labeled as TLR13. The first TLR13 is syntenically conserved across tetrapods (excluding birds and some crocodilians; supplementary fig. S5, Supplementary Material online). The second TLR13 gene is phylogenetically sister to, and reciprocally monophyletic with, TLR21 (fig. 3) and occurs in a separate, highly conserved syntenic region (fig. 3 vs. supplementary fig. S5, Supplementary Material online vs. supplementary fig. S6, Supplementary Material online). Here, we refer to this gene as TLR21-like. The absence of a TLR21-like ortholog in mammals and birds, which include most traditional model systems, may explain why the sequence was not previously identified and reveals a potentially important component of the immune gene repertoire of nonavian sauropsids. The TLR21-like gene is exclusive to nonavian sauropsids (but absent in gecko) and has a highly conserved protein motif pattern (supplementary fig. S2, Supplementary Material online), suggesting the duplication event could have occurred recently (e.g., by gene conversion) or that it arose on the sauropsid stem lineage, was strongly conserved, and was secondarily lost in birds and gecko. This second interpretation is likely given documented loss of TLR8, 9, and 13 in birds (Kahn et al. 2019). In Anolis carolinensis, the LINGO1 (LRR and lg domain-containing 1) gene is in the syntenic region of TLR21-like and contains five LRRs and could potentially be mislabeled or derived from TLR21-like. The Chinese softshell turtle has a tandem duplication of TLR21-like with an identical nucleotide sequence but separate gene model, suggesting a recent lineage-specific duplication (supplementary fig. S10, Supplementary Material online).

The TLR21 protein function is not well characterized, but involves recognizing CpG motifs of DNA viruses (fig. 3). In zebrafish, TLR21 is expressed in the early embryo and within the adult in organs of the innate immune system including intestine, spleen, and kidney (Yeh et al. 2013). Although TLR21 and TLR9 proteins have distinct CpG-DNA recognition sequences, they cooperatively mediate the pathogen-associated response. The duplication of a TLR21-like gene in nonavian sauropsids may have allowed for paralog specialization of tissue-specific expression and CpG-DNA recognition. Further characterization of this TLR21-like gene will be informative.

Conclusion

Comparative analyses of 22 tetrapod genomes revealed active TLR gene gain/loss across nonavian sauropsids, particularly the expansion of TLR8 in turtles, tortoises, and crocodilians and a novel TLR homolog exclusive to nonavian sauropsids, TLR21-like. The results from our three-step approach underscore the importance of integrating synteny, sequence similarity, and semantic searching to study gene family evolution. The syntenic analysis requires contiguous reference genomes and promises new insight into other biologically important gene families. Our findings were enabled by a revised draft of the Mojave Desert tortoise genome (gopAga2.0) that greatly improved assembly contiguity and the de novo annotation of 25,469 genes.

Materials and Methods

Genome Assembly and Annotation

For the gopAga2.0 assembly, the gopAga1.1 assembly (Webster, Dolby, et al. 2018) and skeletal muscle from the original individual were sent to Dovetail Genomics for chromatin-based long-range proximity ligation sequencing with the Chicago® method. These data were reanalyzed together with the gopAga1.0 paired-end and mate-pair Illumina HiSeq reads (SRX2367114–5, SRX2367341, SRX2367455, SRX2367692–3; Tollis et al. 2017) for rescaffolding using the HiRise software, which corrects for misjoins (Putnam et al. 2016; see Supplementary Material online).

We generated a de novo annotation for gopAga2.0 using genome-guided, tissue-specific transcriptomes generated in Trinity v2.5.1 (Grabherr et al. 2011), predicted proteins from the western painted turtle (NCBI PRJNA210179), and protein evidence from UniProtKB/Swiss-Prot database (fig. 1). These data were mapped to the assembly using MAKER v3 (Campbell et al. 2015) followed by three rounds of ab initio gene model training with SNAP v2017-3-1 (Korf 2004) and AUGUSTUS v3.3.2 (Keller et al. 2011) within MAKER.

. 1.

—Genome annotation and assembly. (A) Summary of methodology. (B) Results from chromosome painting of conserved synteny blocks shown between the Mojave Desert tortoise genome (gopAga2.0) and chicken (top right), and green anole lizard (bottom right). The following pairs of scaffolds (left) and chromosomes (right) between Mojave Desert tortoise and chicken are particularly conserved: B-27/28, D-15, E-14, G-11, I-17, O-21, W-22. Black bars show 50-Mb scale. Microchromosomes are labeled “LG”; there are no conserved synteny blocks for LGa, LGd, and LGg. Only homologous linkage blocks are colored.

Bioinformatic Chromosome Painting

To assess synteny conservation, we used bioinformatic chromosome “painting.” We downloaded EMBL genome files from Ensembl for green anole and chicken (Collins et al. 2004; Zhang et al. 2014). To process and reformat the gopAga2.0 genome file, we used a custom Perl script (Gfftofasta.pl modified from https://github.com/ISUgenomics/common_scripts) and subset proteins from the 26 most gene-rich scaffolds in the gopAga2.0 annotation (supplementary table S3, Supplementary Material online) using Biopython v1.73 (Cock et al. 2009). With conversion scripts in the CHROnicle package (v2015), we converted the three species files into SynChro-specific formatting. SynChro computed conserved synteny blocks with delta = 4, which is slightly conservative as it requires four consecutive genes to match across species to be considered a synteny block (the default is 2; Drillon et al. 2014).

TLR Evolution

To study TLR evolution, we used a three-step approach to curate orthology by: 1) searching by sequence similarity, 2) searching by semantic labels, and 3) assessing synteny conservation (supplementary fig. S1, Supplementary Material online). First, using human TLR protein as query sequences, or mouse/chicken sequences when human was unavailable (supplementary table S4, Supplementary Material online), we identified the top ten gene “hits” using BLASTp from the BLAST+ v2.9.0 (Camacho et al. 2009) for each of the seven TLR genes in 22 species: 14 nonavian sauropsids, 4 mammals, 3 birds, and the western clawed frog as outgroup (supplementary table S1, Supplementary Material online and Appendix). From the top ten RefSeq matches, we chose the top hit with an NP tag; if no hit had an NP tag then we used the top predicted hit (XP tag).

Second, we queried the NIH RefSeq database for the text string “Toll-like receptor” with the name of each study species and added to the list any unique gene models that had nonunique gene annotations (e.g., if two separate gene models were both labeled TLR21). Third, and most importantly, we used synteny conservation to verify orthology assignments of these candidate homologs. We visualized the syntenic region of each gene/gene model in every species using the Integrated Genomics Viewer (IGV, v2.3.47; Robinson et al. 2011), manually recorded gene IDs, gene order, and gene direction for five neighboring genes when available, and verified that each gene localized to the expected syntenic region. We excluded long noncoding RNAs, microRNAs, and unidentified loci. For sequences that fell outside the expected syntenic region, if that new syntenic region was the same in another species, it was retained (e.g., TLR21-like), if it fell in a syntenic region that was not recapitulated in other taxa, it was discarded (see Supplementary Material online). This three-step approach provides orthology assignments using complementary sources of evidence and allows for the discovery of novel genes that would be overlooked through sequence similarity alone.

We verified functional protein domains for some TLR sequences of interest using the Simple Modular Architecture Research Tool (SMART; Letunic and Bork 2018). Curated amino acid sequences were aligned using the Gonnet scoring matrix in ClustalW v2.1 (Larkin et al. 2007) with default settings (gap extension penalty = 0.20, gap opening penalty = 10). We visually inspected the amino acid alignments and removed short, highly variable regions that could not be unambiguously aligned. Original and manually curated alignments are available (see Supplementary Material online).

For tree reconstructions, we fit amino acid substitution models for the two subfamilies with ProtTest v3.2 (Darriba et al. 2011) and generated maximum likelihood and Bayesian inference trees on the CIPRES Science Gateway v3.3 (Miller et al. 2010) using the best-scoring models that were implemented in MrBayes (WAG, TLR7 family; VT, TLR11 family; supplementary table S5, Supplementary Material online). We generated maximum likelihood trees using 100 bootstrap iterations in RAxML-HPC v8.2.10 (Stamatakis 2014). In MrBayes v3.2.6 (Ronquist et al. 2012), we generated two runs of 10 million generations each, four chains per run and sampled every 10,000 generations. Runs reached convergence according to the potential scale reduction factor converging on 1; runs were aggregated using sump and sumt commands with the default 25% burn-in discarded. We visually compared topologies and found consistency across ML and BI trees. Markdown tutorials and scripts for all methods are provided (see Supplementary Material online).

Statement of Ethics Approval

Tortoise handling, transport, care, and tissue collection were approved by the Institutional Animal Care and Use Committee at Arizona State University (Protocol No. 13-1319 R). Acquisition and transport of the tortoise was also approved by the United States Fish and Wildlife Service (Recovery Subpermit No. FWSDTRO-1) and the Nevada Department of Wildlife (Export Permit No. S37016).

Click here for additional data file.