Leveraging an existing whole-genome resequencing population data set to characterize toll-like receptor gene diversity in a threatened bird.

Species recovery programs are increasingly using genomic data to measure neutral genetic diversity and calculate metrics like relatedness. While these measures can inform conservation management, determining the mechanisms underlying inbreeding depression requires information about functional genes associated with adaptive or maladaptive traits. Toll-like receptors (TLRs) are one family of functional genes, which play a crucial role in recognition of pathogens and activation of the immune system. Previously, these genes have been analysed using species-specific primers and PCR. Here, we leverage an existing short-read reference genome, whole-genome resequencing population data set, and bioinformatic tools to characterize TLR gene diversity in captive and wild tchūriwat'/tūturuatu/shore plover (Thinornis novaeseelandiae), a threatened bird endemic to Aotearoa New Zealand. Our results show that TLR gene diversity in tchūriwat'/tūturuatu is low, and forms two distinct captive and wild genetic clusters. The bioinformatic approach presented here has broad applicability to other threatened species with existing genomic resources in Aotearoa New Zealand and beyond.

INTRODUCTION

Advances in high‐throughput sequencing technologies and bioinformatic tools are enabling scientists to generate and analyse whole‐genome resequencing data (Auwera et al., 2013 ; Ekblom & Wolf, 2014). This progress has largely increased the measurement of neutral genetic diversity in nonmodel organisms often to inform the conservation management of threatened species (Allendorf et al., 2010; Angeloni et al., 2012; Primmer, 2009). However, studies suggest that neutral genetic diversity may be a poor proxy for functional diversity (Grueber et al., 2015; Marsden et al., 2013; Sommer, 2005). Further, whereas neutral genetic diversity can be used to measure inbreeding and assess inbreeding depression, determining the genetic mechanisms that underlie inbreeding depression requires information about specific functional loci associated with adaptive or maladaptive traits (Kohn et al., 2006; Mable, 2019; Ouborg et al., 2010; Won et al., 2021). When whole‐genome resequencing data and a reference genome for a population of interest already exists, bioinformatic tools may be used to analyse these data to characterize functional genetic diversity (Hoelzel et al., 2019; Kohn et al., 2006).

There is a growing interest in characterizing the functional diversity of immune genes within intensively managed species to improve outcomes and reduce the burden of wildlife disease (Morris et al., 2015; Zhu et al., 2020). Low immune gene diversity increases susceptibility to a variety of pathogens (Spielman et al., 2004). Inbreeding leads to an overall decrease in genetic diversity at immune genes, loss of rare and potentially advantageous alleles, and a decreased ability to adapt to novel or rapidly evolving pathogens (Altizer et al., 2003; Spielman et al., 2004). Research from both wild and laboratory populations show that inbreeding contributes to an increased parasite load, a greater susceptibility to pathogens, a higher likelihood that individuals will act as disease reservoirs, and ultimately, higher rates of mortality due to disease (Acevedo‐Whitehouse et al., 2003; Ross‐Gillespie et al., 2007; Whitehorn et al., 2011; Whiteman et al., 2005). Further, increased immune gene diversity may help populations respond more rapidly to pathogens (Bonneaud et al., 2012; Davies et al., 2021).

Toll‐like receptor (TLR) genes are an innate immune gene family involved in pathogen recognition and immune response (Acevedo‐Whitehouse & Cunningham, 2006; Grueber et al., 2012; Vinkler & Albrecht, 2009). TLRs are an ancient part of the immune system, present in nearly all multicellular organisms (Singh et al., 2003). They are located either on the cell membrane or intracellularly on the membrane of lysosomes or endosomes within innate immune or somatic cells (Behzadi et al., 2021; Takeda & Akira, 2005). TLRs recognize conserved patterns of pathogens by binding residues, known as antigens, from pathogens that enter the body (Singh et al., 2003). Pathogen recognition by TLRs is necessary for the proper activation and direction of the adaptive immune response (Chen et al., 2021; Clark & Kupper, 2005; Kawasaki & Kawai, 2014; Pasare & Medzhitov, 2005; Xia et al., 2021). Studies have found that TLR gene diversity tends to be low in small, highly threatened populations compared to large populations of species of least concern (Alcaide & Edwards, 2011; Dalton et al., 2016; Knafler et al., 2017; Morrison et al., 2020). In previous studies, having higher TLR gene diversity has been associated with greater disease resistance, survival, and reproduction (Antonides et al., 2019; Bateson et al., 2016; Davies et al., 2021; Heng et al., 2011; Podlaszczuk et al., 2021; Quéméré et al., 2021). Having a robust TLR gene diversity within a population allows for selection of the most beneficial variant(s) to effectively deal with pathogens. Also, immune gene diversity may not correlate with neutral measures of genome‐wide diversity, so directly characterizing TLR genes will explicitly reveal how diversity at innate immune genes are affected in highly threatened avian populations (Grueber et al., 2015; Hartmann et al., 2014; Marsden et al., 2013).

To date, TLR gene diversity has been mostly characterized using species‐specific TLR primers and PCR protocols to perform SNP genotyping of individuals in a population (Allendorf et al., 2010; Dalton et al., 2016; Grueber et al., 2015; Lara et al., 2020; Morrison et al., 2020). However, these approaches can be expensive and time‐consuming. When whole‐genome resequencing data and a reference genome for the species of interest already exists, researchers can use bioinformatic tools to analyse these existing resources to identify TLR genes and characterize TLR gene diversity (Hoelzel et al., 2019; Wang et al., 2005).

TLR genes may be more suited for identification using bioinformatic tools, in comparison to other immune genes like MHC, because they are relatively conserved and gene duplications are rare and well‐defined (Grueber et al., 2015). The TLR protein is composed of three protein domains: the extracellular binding domain, transmembrane protein, and intracellular toll/interleukin 1 (TIR) signalling domain (Yilmaz et al., 2005). There is mostly conserved evolution and synonymous substitutions within TLR sequences, and the majority of variation is within the region coding for the extracellular binding protein, which is what comes into contact with antigens (Werling et al., 2009). The binding domain has a pattern of leucine‐rich repeats, and so is also known as the leucine‐rich repeat (LRR) domain (Alcaide & Edwards, 2011). Inserts of lecuines within the LRR domain may impact pathogen recognition (Offord et al., 2010), and single nucleotide polymorphisms (SNPs) within this region can also affect the binding affinity of TLRs (Keestra et al., 2008; Matsushima et al., 2007). In contrast, the TIR signalling domain is largely conserved across the TLR family, and phylogenetically among related species (Beutler & Rehli, 2002; Narayanan & Park, 2015; Yilmaz et al., 2005). The implications of these protein structures are twofold: (1) the pattern of repeats within the LRR domain and conserved TIR domain are ideal for bioinformatic identification due to the similarity of sequences within taxonomically related species, and (2) the sequences within the LRR region are variable and may contain SNPs that are adaptive for the recognition of particular pathogens.

Advances in bioinformatic tools utilizing comparative genomics and sequence similarity comparison and the growth of online repositories and databases like NCBI with extensive genomic data that includes functional genes from a diversity of species, together facilitate the identification of functional genes within the genomes of nonmodel organisms (Feng et al., 2020; Grueber, 2015; Ren et al., 2021). However only a few studies to date have used bioinformatic tools to identify functional genes and characterize functional genetic diversity within a whole‐genome resequencing population data set (Brandies et al., 2020; Zhang et al., 2014). Previous research on avian TLRs have used tools such as NCBI Basic Local Alignment Search Tool (blast) to design primer sequences to amplify these genes (Chávez‐Treviño, 2017; Grueber et al., 2015; Morrison et al., 2020) and to compare sequences among closely‐related species (Raven et al., 2017), so using established methods for primer design may also help in the identification of TLR genes (Yilmaz et al., 2005). A study on immune gene diversity in Tasmanian devils used bioinformatics to characterize diversity of several immune genes previously identified through a combination of transcriptome analysis and comparative genomes (Morris et al., 2015). Research using genomes from the B10K consortium highlights the efficacy of comparative genomics to identify orthologues and conserved regions using bird species in the same taxonomic class (Feng et al., 2020; Zhang et al., 2014). Increasingly, immune gene databases are available to identify and annotate immune genes within specific taxonomic groups (Grueber, 2015; Mueller et al., 2020; Wong et al., 2011). These resources provide an avenue for characterizing immune gene diversity in nonmodel organisms.

Tchūriwat’/tūturuatu/shore plover (Thinornis novaeseelandiae) is a threatened shorebird endemic to Aotearoa New Zealand (Robertson et al., 2021). Tchūriwat’/tūturuatu is thought to have been widespread on the coastal areas of mainland Aotearoa, but populations decreased dramatically during the 19th century because of introduced mammalian predators (Davis, 1994). Tchūriwat’/tūturuatu was eventually extirpated from the mainland and the remaining wild population is now restricted to the Chatham Islands archipelago (Davis, 1994). The current wild population totals 234 birds and is mostly confined to a single predator‐free island, Hokorereoro/Rangatira (Department of Conservation, Shore Plover Recovery Group herein DOC SPRG). A conservation breeding program was established between 1991–1996 with eggs brought from Hokorereoro/Rangatira and was augmented in 2003 by a single adult male from Western Reef (Dowding et al., 2005). At the time of this study, the captive population had not been augmented since 2003 and the number of individuals in the captive population was 41. Based on the pedigree of the captive population, individuals were highly related to one another (average mean kinship = 0.119) and mean inbreeding was relatively high (0.048) (SPRG, unpublished data).

Captive birds are also highly vulnerable to contracting avian pox (DOC SPRG). Avian pox is caused by avipoxvirus (APV), a large dsDNA virus with hundreds of strains that can infect a wide range of bird species (Bolte et al., 1999; Boyle, 2007). Infection with APV causes lesions on the body at the site where it enters the skin and on feather‐free areas of the body. APV is generally transmitted through insect vectors, but if there is a break in epithelial integrity, it can also spread through contact with an infected bird or via contact with shared objects like feeders used by infected birds (Hansen, 1999). In rare instances, APV can also be transmitted through the inhalation of viral particles (Hansen, 1999). Aspects of captive breeding programmes like shared aviaries and resources make them particularly vulnerable for the spread of infection (van Riper & Forrester 2007).

In addition to being highly transmissible within a captive setting, avian pox infections in captive birds are often severe, can last for several months, and may never resolve (DOC SPRG). In contrast, wild tchūriwat’/tūturuatu contract mild avian pox infections which are most often cleared in 1–2 weeks (DOC SPRG). Longer infection periods for captive birds allow pox lesions more time to develop, making it more likely that birds will contract secondary bacterial infections which increase the severity of and rate of mortality from APV (Alley & Gartrell, 2019; Gartrell et al., 2002; Hansen, 1999; Weli & Tryland, 2011). While death from the virus remains low, the severity of infection is high. Also, chicks and juveniles are more likely to contract APV (DOC SPRG), and without a fully developed adaptive immune system, the innate immune system may be especially important for young birds in fighting off poxvirus infections (Fellah et al., 2008; Palacios et al., 2009). Protecting juveniles from infection is important for the timely release of young birds into the wild (Dowding, 2013, DOC SPRG).

Ongoing research in an aligned project based on approximately 50 K single nucleotide polymorphisms (SNPs) generated using a reduced‐representation sequencing approach shows that wild and captive tchūriwat’/tūturuatu populations are genetically distinct (I. Cubrinovska, unpublished data). These SNPs reflect genome‐wide diversity, but may also be indicative of differences at the gene level, including functional genes like TLRs that contribute to the immune response to poxvirus. If the wild population–which appears to be less susceptible to avian pox (DOC SPRG)–has higher immune gene diversity, then genetic rescue of the captive population through augmentation with wild individuals may improve future disease outcomes.

Sourcing individuals from genetically diverse or different populations is recognized as a way to improve both genome‐wide diversity and immune gene diversity (Glassock et al., 2021; Grueber et al., 2017; Heber et al., 2013; McLennan et al., 2020). Leveraging existing genomic resources and bioinformatic tools, we identify TLR genes and characterize TLR gene diversity in captive and wild tchūriwat’/tūturuatu. In addition to informing future research to determine whether a recent translocation of tchūriwat’/tūturuatu from the wild to captivity has increased TLR gene diversity and an associated immune response to the avipox infection, the bioinformatic approach presented here is broadly applicable to threatened birds with existing genomic resources.

MATERIALS AND METHODS

Sample extraction and sequencing

Originally generated for an aligned project (see Benefit‐sharing statement), we used short‐read sequence data for 66 tchūriwat’/tūturuatu (39 captive, 27 wild), including a wild male named Maui that was chosen for the tchūriwat’/tūturuatu reference genome.

Blood samples for captive birds were collected during routine health checks by Brett Gartrell and Isaac Conservation and Wildlife Trust staff. Blood samples for wild birds were collected as part of the 2018/2019 breeding season juvenile banding trip to Maung’ Rē/Mangere and Hokorereoro/Rangatira in the Chathams Island archipelago. Blood was taken from juveniles and adults and stored in lysis buffer for shipment. All samples were subsequently stored at −0°C until extraction. High quantity and quality DNA was extracted using a tailored lithium chloride extraction method (Galla, 2019). Extractions were assessed for quality by running 2 μl of DNA on a 2% agarose gel. A Qubit 2.0 fluorometer (Fisher Scientific) was used for DNA quantification. Blood samples and DNA extractions are stored at the University of Canterbury on behalf of the Shore Plover Recovery Programme.

Libraries were prepared with the Nextera DNA Flex Library Prep Kit according to the manufacturer's specifications and sequenced across one lane of an Illumina Novaseq 6000 to achieve an average coverage depth of approximately 10× (Galla et al., 2020), excluding the individual used for the reference genome which was sequenced to an average depth of approximately 100×. Libraries for 21 individuals that had low depth of coverage were subsequently resequenced as above to achieve an average coverage depth of at least 10×.

Reference genome assembly

fastqc version 0.11.8 (FastQC, 2015) was used to evaluate the quality of the raw Illumina data and assess potential sample contamination. Initial read trimming was performed using trimgalore version 0.6.6 (Kreuger, 2021), using pair‐end mode, a minimum length of 54 bp, and with the ‐‐nextseq two‐colour chemistry option. A kmer abundance plot was created using jellyfish 2.3.0 (Marçais & Kingsford, 2011) with a kmer length of 31 prior to assembly to assess heterozygosity and contamination. The genome was assembled in two stages using two assemblers, masurca version 3.3.4 (Zimin et al., 2013), and meraculous version 2.2.6 (Goltsman et al., 2017). While the assemblies with MaSuRCA give the best contiguity, the assembler can have problems with two‐colour chemistry Illumina data when finding runs of low‐quality bases or reads with a high proportion of “G” bases (De‐Kayne et al., 2020). In those cases, the assembler tends to create regions where the consensus is simply a string of “G” bases and/or unknown “N” bases, as a result of both the error‐correction module and the superread creation during the assembly process. Meraculous is much less prone to these artefacts due to their agnostic approach to error correction, which relies on coverage and bubble resolution in de Bruijn graphs, at the expense of contiguity when compared to MaSuRCA. Thus, both assemblers were combined.

In the first stage, MaSuRCA was used to create the main assembly, using the trimmed reads padded with low quality “N” bases to a uniform length of 150 bp. Parameter adjustments include a grid batch size of 500,000,000, a Jellyfish hash size of 8,000,000,000, use of the Celera assembler for the final step, and the inclusion of scaffold gap closing. All other parameters were set to default for nonbacterial Illumina assemblies. For Meraculous, a kmer size of 61 with a minimum depth cutoff of five was selected, and the assembly was run on Diploid mode 1 (to create a single haploid reference). Finally the Meraculous assembly was used to correct the ambiguities and gaps in the MaSuRCA assembly by first aligning both genomes using Last version 980 (Kielbasa et al., 2011), filtering for contigs with low quality regions (due to the presence of poly‐C/poly‐G regions right next to a poly‐N region) using awk, then using the merge_and_replace.pl Perl script available at https://github.com/Lanilen/SemHelpers. The hardware used for the assembly was a 12‐core, 24‐thread workstation with 128 Gb of RAM. Assembly took 5 days for MaSuRCA and 10 h for Meraculous.

TLR identification and characterization overview

Figure 1 provides an overview of the methods used for TLR identification and characterization. We used a virtual machine (VM; 16 vCPUs, 9 TB of memory, 128 GB of RAM) in the host Research Compute Cluster (RCC) at the University of Canterbury.

Figure 1

Workflow schematic to show how bioinformatic tools and genomic resources are used to identify and characterize TLR genes in threatened birds. See text for details. bcf, BIM collaboration format; LRR, leucine rich repeat; ORF, open reading frame; SNP, single nucleotide polymorphism; TLR, toll‐like receptor. All remaining abbreviations refer to names of bioinformatic tools utilized

TLR identification

blast alignment

To identify TLR genes, we used NCBI Basic Local Alignment Search Tool (blast), given its prior use to design primer sequences that amplify TLR genes (Chávez‐Treviño et al., 2017; Yilmaz et al., 2005) and previous comparisons of TLR sequences among closely‐related avian species (Mueller et al., 2020; Raven et al., 2017). We used bird TLR nucleotide (BLASTn) and protein (tBLASTn) sequences from close relatives and/or species with high quality genome assemblies to search for similar nucleotide sequences within the reference genome for tchūriwat’/tūturuatu (Madden, 2013). Unlike other analyses that focused solely on the most variable LRR binding region, we chose to identify the entire sequence for each TLR, because it provided additional confirmation that the whole TLR was captured. For inclusion in our search through the tchūriwat’/tūturuatu reference genome, TLR sequences were either characterized for other species in the laboratory through targeted amplification and sequencing, or identified using the NCBI genome annotation pipeline (Table 1).

TABLE 1

Reference species for toll‐like receptor (TLR) identification in tchūriwat’/tūturuatu, listed alongside the taxonomic group shared with tchūriwat’/tūturuatu. For each referece species, TLR sequences used for comparison were either identfied through the use of the NCBI annotation pipeline or through amplification and sequencing of targeted TLRs in the laboratory

Reference TLR species	Shared taxonomic group	Method of TLR identification
Killdeer (Charadrius vociferus)	Family Charadriidae	NCBI genome annotation a
Black‐headed Gull (Chroicocephalus ridibundus)	Order Charadriiformes	Targeted amplification and sequencing b
Chicken (Gallus gallus)	Class Aves	Targeted amplification and sequencing c
Zebra Finch (Taeniopygia guttata)	Class Aves	NCBI genome annotation excluding targeted amplification and sequencing of TLR4 d

The NCBI Eukaryotic Genome Annotation Pipeline, https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/.

Podlaszczuk et al. (2021).

Yilmaz et al. (2005).

Vinkler et al. (2009).

BLAST uses an algorithm to map input query sequences to the most similar region within the tchūriwat’/tūturuatu genome. If BLAST is able to map the input query sequence to a region or regions within the database, it outputs the region(s) as a list. The list is ordered by the expect value (e‐value) of each alignment, which is the number of alignments with a similar score that are expected by chance. The lower the e‐value of the alignment, the less likely that the alignment is due to chance and the more likely it is that the alignment reflects a biological similarity between the sequences (NCBI, 2020). In addition to the e‐value, the percent query cover and percent identity of the alignment help the user to judge the quality of each alignment. The percent query cover is percent of the query sequence that is aligned to a database sequence, out of the whole length of the query sequence, regardless of the identity of the bases within a sequence. The percent identity is how many of the bases within the query sequence alignment match the bases within the database sequence, out of the whole length of the query sequence. The closer both of these percentages are to 100%, the greater the likelihood that the alignment is biologically relevant.

To further ensure that agreement between the query sequence and individual target sequences was not due to chance, two strategies were implemented: (1) a comparative approach seeking agreement among multiple species' TLR sequences aligning to a specific region in the tchūriwat’/tūturuatu genome; and (2) the reciprocal best hit (RBH) test. In the first instance, if the same region in the tchūriwat’/tūturuatu genome was being mapped to with the same type of TLR gene from different bird species, there was greater confidence in that region. In the second, the RBH test is an approach that has been used within comparative genomics to confirm the identity of orthologs in nonmodel species using well‐annotated reference genomes (Kristensen et al., 2011). This limited the test to just chicken and zebra finch genomes, since both are high‐quality genome assemblies and are well‐annotated. This test involves taking the “best hit” region identified within tchūriwat’/tūturuatu using the TLR sequence from a reference species and BLASTing it against the genome of that reference species (Irizarry et al., 2016). If the “best hit” of this BLAST search matches the original reference TLR, it provides greater support that the genes are orthologs of one another (Kristensen et al., 2011). Each preliminary TLR sequence identified in tchūriwat’/tūturuatu passed the RBH test, mapping to the original chicken and zebra finch reference TLR sequences within each genome. This helped to confirm not only that a TLR gene was likely identified, but also that the specific TLR gene of interest had been located (i.e., a BLAST search with preliminary TLR3 sequence brought up only the TLR3 sequence in reference species).

After conducting these tests of gene alignment, each preliminary TLR sequence was BLASTed against the remainder of the tchūriwat’/tūturuatu genome to find whether there was alignment of this sequence to other regions within the genome. This test was conducted to determine whether there were nontarget sequences within the tchūriwat’/tūturuatu genome that may align to these TLR sequences during the process of whole‐genome alignment. There was little cross‐alignment of one type of TLR mapping to other TLR genes (e.g., TLR1A reference aligning to the TLR3 region). When this did happen, it was often within pairs of duplicated TLRs (TLR1A & TLR1B, TLR2A & TLR2B), and it was only a partial alignment (at most 50%), so these cross alignments were easy to distinguish from the true alignment. False positives were also investigated within the genome, to find out whether there were TLR sequences aligning to other, non‐TLR regions. There were no false positive sequences within the tchūriwat’/tūturuatu genome. Based on these tests, the Burrows‐Wheeler Aligner version 0.7.17 (bwa) (Li & Durbin, 2009) with a default maximum mismatch value of 4% for read alignments (allowing only 4% of bases to differ in identity within sequence alignments) would be sufficient to prevent the cross‐alignment of duplicated TLRs.

Protein analysis

Once the preliminary TLR regions passed these quality control measures, each sequence was entered into the NCBI Open Reading Frame (ORF) finder (https://www.ncbi.nlm.nih.gov/orffinder). This tool searches for reading frames within the query DNA sequence and transcribes it into a protein sequence with each possible reading frame. The result is a graphic showing the protein sequence resulting from each reading frame that is used for transcription (Wheeler et al., 2003). Then we used BLASTp searches with the resultant protein as the input sequence and the BLAST protein database as comparison sequences. When the correct reading frame for the protein was identified, the search would bring up the reference TLR protein sequences that we had originally aligned with.

Two additional web applications were used to examine TLR protein sequences and investigate the protein products of TLR genes. The first was LRRfinder (Offord et al., 2010), which uses a database of toll‐like receptor TLR sequences acquired from NCBI to identify the LRR regions within a TLR (LRRs, LRRNT, LRRCT, transmembrane protein, TIR signalling domain). The search tool focuses on predicting potential LRR regions, because this region is the most variable and may affect the binding specificity and affinity of TLRs (Keestra et al., 2008; Matsushima et al., 2007). All searches revealed multiple LRR regions within the sequences, a transmembrane protein region in almost all sequences, and a TIR domain in all sequences. Also the LRR identifications matched the type of TLR gene they were from (i.e., all LRRs for TLR3 matched TLR3 LRRs in reference species).

To further visualize these protein products, we used the simple modular architecture research tool (smart version 9.0) (Letunic et al., 2021; Schultz et al., 1998). SMART contains a protein database that uses markov modes to identify protein domains within an input protein sequence by calculating the expected value (e‐value) SWise score, which is the output of the established SWise protein search algorithm (Birney et al., 1996), for each alignment between the query sequence and sequences in the SMART database (Letunic et al., 2021; Schultz et al., 1998). The protein domains assigned with low e‐values are less likely to be assigned by chance, and at a predetermined e‐value threshold, the protein domain will appear as a visual block on a 2D schematic of the protein. This visualization revealed whether we had captured all expected protein domains within the TLR protein sequences and had not missed parts of the sequence. This also facilitated a visual comparison between the tchūriwat’/tūturuatu TLR protein schematics to the TLR SMART protein schematics in chicken (Temperley et al., 2008).

TLR characterization

Population resequencing alignment and analysis

Once TLR regions were identified within the genome, whole‐genome alignment of population resequencing data to the tchūriwat’/tūturuatu reference genome was performed to characterize TLR SNP diversity. To assess read quality, fastqc version 0.11.9 (FastQC, 2015) was run to identify low quality regions and baseline quality scores, and then trim galore version 0.6.5 (Kreuger, 2021) was used to trim Illumina paired end 150 bp reads and remove low quality reads. The two‐colour chemistry option was chosen, which supports trimming and removal of low‐quality sequences of non‐G bases and improves the removal of low quality G's by removing them regardless of their quality (Kreuger, 2021).

Trimmed reads were aligned to the reference genome with bwa version 0.7.17 (Li & Durbin, 2009). samtools version 1.10. (Li et al., 2009) was used to sort BAM files prior to SNP discovery. A custom perl script (“split_bamfiles_tasks.pl”) (Moraga, 2018) was used to split the bam files into chunks that could be processed more quickly with bcftools version 1.11 (Li et al., 2009). bcftools mpileup was run with annotations GT, PL, DP, SP, ADF, ADR, AD to allow for downstream filtering. vcftools version 0.1.16 (Danecek et al., 2011) was used to filter the data set as follows: minor allele frequency (maf) ≥0.05, Phred‐score (quality) >20, max‐missingness = 0.90, minimum depth >5, minGQ >10, and maximum depth <200. bcftools was also used to filter the vcf for strand bias using the parameter of strand‐bias adjusted Phred‐score <60. Hardy–Weinberg equilibrium filtering was not applied because two assumptions are likely to be violated in this study (no selection and random mating). In addition, vcftools was used to examine the site and individual depth and missingness and to remove individuals with low depth and high missingness from further analysis (n = 1 captive individual). Final SNPs were analysed in Geneious Prime 2020 (https://www.geneious.com/) to determine whether they made synonymous or nonsynonymous changes to TLR protein products. Nonsynonymous SNPs were further investigated in Geneious to analyse the amino acid change that resulted from each SNP. Then, a preliminary analysis of nonsynonymous SNPs was conducted to determine which physiochemical changes that might occur due to the change in amino acid identity, using the well‐defined physicochemical attributes of charge and polarity.

beagle version 5.2 (Browning et al., 2021) was used with default settings to phase haplotypes for both populations. Each TLR gene was phased separately by phasing along the whole contig that contained that TLR gene. For each individual, genotype sequences were generated using bcftools. To increase quality and decrease errors within haplotypes, sites with low depth and high strand bias were set to missing using bcftools +setgt with options ‐i ‘FORMAT/DP<5’ and ‘FORMAT/SP>60’. Then, the filtered VCF was used to produce TLR consensus sequences for each individual in the population. samtools faidx was used to target each TLR region in the genome, and bcftools consensus was used with parameters ‐M to output any missing genotypes as “N”, and ‐H 1pIu and 2pIu to produce both phased haplotypes for an individual. If individuals had any sites within the TLR sequence where they were missing a genotype, they were not included in the haplotype analysis for that TLR gene. While the accepted standard for SNP analysis is that nonsynonymous SNPs have significant influence on the functionality of a protein, synonymous SNPs may also change protein function (Sauna & Kimchi‐Sarfaty, 2011), and within TLR genes, synonymous SNPs have been associated with disease resistance and immune response (Cho et al., 2013; Junjie et al., 2012). Given this, both synonymous and nonsynonymous SNPs were used in the construction of haplotypes.

The resulting TLR consensus sequence for each individual was output into a population mega file for either the captive or wild population, and the process was repeated for every polymorphic TLR gene. The resulting mega files were imported into dnasp v. 6.12 (Rozas, 2017) and were analysed to calculate haplotype diversity, nucleotide diversity (𝛑), and Tajima's D. Haplotype nexus files for each TLR were also created in DNAsp and imported into popart (Leigh & Bryant, 2015) to construct minimum‐joining haplotype networks. We also visualized differences in TLR haplotypes between populations with a PCA. Adegenet (Jombart, 2008) was used to construct the PCA with scaled mean haplotype frequencies for all TLR genes, using code adapted from the introductory tutorial to adegenet (Jombart, 2015).

RESULTS

Nine TLRs were identified within the tchūriwat’/tūturuatu genome: TLR1A, TLR1B, TLR2A, TLR2B, TLR3, TLR4, TLR5, TLR7, and TLR21. A partial sequence of the TLR15 gene was identified within the genome, but was not included in the analysis. The protein schematic structure for each TLR identified is shown in Figure 2.

FIGURE 2

SMART protein schematics for each tchūriwat’/tūturuatu toll‐like receptor identified (available via licence: CC BY 2.0). Key shows the visual representation for each protein domain. The pink boxes are areas of low compositional complexity

There were 28 SNPs in the TLR genes of the captive tchūriwat’/tūturuatu population (n = 38) (Table 2). The SNPs were unevenly distributed among TLR genes, with two that were monomorphic (TLR2B, TLR 21), three with one SNP (TLR1A, TLR1B, TLR2A), one with two SNPS (TLR4), one with three SNPs (TLR3), one with 7 SNPs (TLR5), and one with 13 SNPs (TLR7). Out of the total SNPs, half (14) were synonymous, and the remaining half (14) were nonsynonymous. The majority of SNPs (22) were within the LRR binding domain. The majority of nonsynonymous SNPs (12) were also located in the LRR region and half of them (7) caused a physicochemical change due to the change in identity of the amino acid (Table 3). A total of 24 haplotypes were observed across all TLR genes in the captive population, with the TLR genes that have higher SNP diversity also having a higher number of inferred haplotypes (Table 2). Tajima's D was nonsignificant for most loci, except for TLR5 and TLR7, which were both significantly positive (p < .01).

TABLE 2

Tchūriwat’/tūturuatu toll‐like receptor gene diversity statistics and comparisons for captive (n = 38) and wild (n = 26) tchūriwat’/tūturuatu populations

Locus	Length bp (aa)	SNPS (syn:Nsyn)			LRR SNPs (syn:Nsyn)			# inferred haplotypes (sample size)			Nucleotide diversity (𝛑)		Haplotype diversity		Tajima's D
		Captive	Wild	Shared	Captive	Wild	Shared	Captive	Wild	Shared	Captive	Wild	Captive	Wild	Captive	Wild
TLR1A	2289 (761)	1 (0:1)	3 (2:1)	1 (0:1)	1 (0:1)	3 (2:1)	1 (0:1)	2 (35)	3 (25)	1	0.00011	0.00042	0.248	0.483	0.24403	0.91554
TLR1B	1857 (592)	1 (1:0)	1 (1:0)	1 (1:0)	1 (1:0)	1 (1:0)	1 (1:0)	2 (35)	2 (25)	2	0.00024	0.00023	0.437	0.429	1.37017	1.23557
TLR2A	2423 (806)	1 (1:0)	0	0	1 (1:0)	0	0	2 (37)	1 (26)	1	0.00018	0	0.425	0	1.4606	0
TLR2B	2156 (717)	0	1 (1:0)	0	0	0	0	1 (38)	2 (26)	1	0	0.0001	0	0.216	0	−0.0466
TLR3	2367 (783)	3 (3:0)	2 (2:0)	2 (2:0)	1 (1:0)	1 (1:0)	1 (1:0)	4 (38)	3 (26)	2	0.00045	0.00029	0.478	0.521	0.79912	0.46712
TLR4	2259 (752)	2 (1:1)	1 (1:0)	1 (0:1)	1 (0:1)	0	0	3 (38)	2 (24)	2	0.00027	0.00013	0.537	0.284	0.77496	0.35171
TLR5	2585 (860)	7 (4:3)	8 (4:4)	7 (4:3)	4 (2:2)	4 (2:2)	4 (2:2)	3 (35)	4 (23)	3	0.00126	0.00146	0.551	0.673	3.08884*	2.89296*
TLR7	3141 (1045)	13 (4:9)	13 (4:9)	13 (4:9)	13 (4:9)	13 (4:9)	13 (4:9)	6 (38)	6 (26)	6	0.000206	0.00206	0.665	0.681	4.05164*	3.69096*
TLR21	1649 (548)	0	0	0	0	0	0	1 (38)	1 (26)	1	0	0	0	0	0	0
Total Private		3 (2:1)	4 (3:1)	–	2 (1:1)	2 (2:0)	–	5	5	–
Total		28 (14:14)	29 (15:14)	25 (11:14)	22 (9:13)	22 (10:12)	20 (8:12)	24	24	19

Note: Asterisks denote significance at p < .01.

Abbreviations: bp(aa), base pairs (amino acids); syn:nsyn, ratio of synonymous to nonsynonymous SNPs.

TABLE 3

Amino acid analysis of nonsynonymous SNPs within TLR genes of captive and wild tchūriwat’/tūturuatu populations

TLR	Binding specificity a	AA site (codon)	TLR region	AA change	Physicochemical change	Population
TLR 1A	Bacteria, liproproteins	503	LRR	Asp/His	Charge	Captive, Wild
TLR 4	LPS, gram negative bacteria	372	LRR	Asn/Asp	Charge	Captive
TLR 5	Flagellin	9	LRR	Phe/Leu	None	Captive, Wild
		326	LRR	Val/Ile	Charge	Wild
		664	TIR	His/Arg	None	Captive, Wild
		843	TIR	Lys/Glu	Charge	Captive, Wild
TLR7	ssRNA, virus	4	LRR	Ala/Pro	None	Captive, Wild
		79	LRR	Thr/Ile	Polarity	Captive, Wild
		100	LRR	Met/Leu	None	Captive, Wild
		411	LRR	Leu/Phe	None	Captive, Wild
		459	LRR	Gln/His	Charge	Captive, Wild
		465	LRR	Ala/Val	None	Captive, Wild
		467	LRR	Glu/Gly	Charge, Polarity	Captive, Wild
		469	LRR	Asn/Ser	None	Captive, Wild
		755	LRR	Gln/Lys	Charge	Captive, Wild

Abbreviations: AA, amino acid, standard abbreviations for amino acids used; LRR, leucine rich repeat binding domain; SNPs, single nucleotide polymorphisms; TIR, TIR signalling domain; TLR‐ toll‐like receptor.

Alcaide and Edwards (2011).

A total of 29 SNPs were observed in the TLR genes of the wild tchūriwat’/tūturuatu population (n = 26) (Table 2). There were two monomorphic genes (TLR2B, TLR 21), three with one SNP each (TLR1B, TLR2A, TLR4), one with two SNPS (TLR3), one with three SNPs (TLR1A), one with eight SNPs (TLR5), and one with 13 SNPs (TLR7). Out of the total SNPs, slightly more than half (15) were synonymous, and the remaining 14 were nonsynonymous. The majority of SNPs (22) were within the LRR binding domain. The majority of nonsynonymous SNPs (12) were also located in the LRR region and half of them (7) caused a physicochemical change due to the change in identity of the amino acid (Table 3). A total of 24 haplotypes were observed across all TLR genes in the wild population (Table 2). Tajima's D was nonsignificant for most loci, except for TLR5 and TLR7, which were both significantly positive (p < .01).

The wild and captive populations of tchūriwat’/tūturuatu both had low overall TLR SNP and haplotype diversity. Each had a similar ratio of nonsynonymous to synonymous SNPs. Across all loci, the populations shared 25 SNPs and 19 haplotypes (Table 2). The captive population had three private SNPs and five private haplotypes in four genes (TLR1A, TLR2A, TLR3, TLR4). The wild population had four private SNPs and five private haplotypes in four genes (TLR1A, TLR2B, TLR4, TLR5). Haplotype networks reveal there are closely related haplotypes for each TLR gene, and that many of the TLR haplotypes shared between captive and wild tchūriwat’/tūturuatu occur at different frequencies with each population (Figure S1). The PCA reveals two relatively distinct genetic clusters corresponding to wild and captive populations (Figure 3).

FIGURE 3

Principal component analysis genetic clustering of captive (n = 38) and wild (n = 26) tchūriwat’/tūturuatu. Produced using scaled mean allele frequencies for each individual (see text for details)

DISCUSSION

The bioinformatic approach developed here was successful for identifying TLR genes and characterizing TLR gene diversity in tchūriwat’/tūturuatu. This research provides a critical first step towards using TLR gene diversity to inform conservation action for a threatened Aotearoa New Zealand endemic bird. We were able to identify and characterize the same number of or more TLR genes as compared to similar studies done using amplicon sequencing of TLRs within threatened species (Dalton et al., 2016; Grueber et al., 2015; Morrison et al., 2020). Further, as the generation of whole‐genome resequencing data for threatened species becomes routine (Allendorf et al., 2010; Forcina & Leonard, 2020), identifying and characterizing TLRs using existing genomic resources will be more cost effective and more efficient than traditional amplicon sequencing.

All complete TLRs identified display the structure of TLR proteins in other species: a region of leucine rich repeats (the binding domain), then a transmembrane protein (carboxyl‐terminal tail) where the TLR sits within the cell or lysosome membrane, and then the TIR (toll‐like/interleukin receptor) signalling domain (Kannaki et al., 2010; Yilmaz et al., 2005). While LRRs vary in number between TLR sequences, the overall structure and approximate amount of LRRs within each type of TLR is similar to what is seen in chicken TLRs (Temperley et al., 2008). A partial sequence of TLR 15 was identified, but this sequence did not include the TIR signalling domain. The sequence ended prematurely at the end of a contig, so the scaffolding of the short‐read genome may have prevented identification of the full gene. Morris et al. (2015) encountered a similar problem identifying immune gene sequences when they were fragmented or split between contigs in the genome assembly.

In nonpasserine avian species, there are most often 10 avian TLRs, eight of which are orthologous to other vertebrate TLRs (TLR1A/B, TLR2A/B, TLR3, TLR4, TLR5, TLR7), one that is orthologous to bony fish and Xenopus (TLR21), and one that is unique to reptiles and birds (TLR15) (Alcaide & Edwards, 2011; Grueber et al., 2014). Studies of avian TLR evolution show a pattern of both gene loss and duplication in these regions (Kannaki et al., 2010; Temperley et al., 2008; Velová et al., 2018). Recent research suggests there is also a duplication of TLR7 in some avian taxa. The duplicated TLR7 is thought to have a similar function, though with slight difference, and this is an area of ongoing research (Raven et al., 2017). To date, it has been found in Charadriiformes, Cuculiformes, Mesiornithiformes, and some Passeriiformes (Velová et al., 2018). Tūturuatu are within the order Charadriiformes, and this species is likely to have a duplication of TLR7 and thus have 11 TLRs total (Raven et al., 2017; Velová et al., 2018).

However, we did not find evidence of the duplication of TLR7 within tchūriwat’/tūturuatu. After identifying the first gene coding for TLR7, this region was removed from the genome and the BLAST search with reference TLR7 sequences was repeated with the remainder of the genome. This second search was done to ensure that the reference TLR sequence would not align to the previously identified TLR7 region. There was no additional alignment to other sequences in the genome. It is possible that the duplication of TLR7 within the reference genome is not fully resolved, meaning the duplicated genes are collapsed into one region. Alternatively, the duplicated region may be incomplete with gaps in the sequence or it may be of low quality, so that it is either partially or fully cut out of the assembly. In the case of a low quality or partial assembly, the duplicated region would not be easily identified through the use of BLAST. In either case, only one TLR7 gene was identified. If the duplication does exist and is of low quality within the tchūriwat’/tūturuatu genome, this means the population resequencing data may not properly align to the low‐quality TLR7 gene and may instead align to the other TLR7 gene. In this case, SNPs found within TLR7 may not be true SNPs, but instead may be artefacts of the duplicated TLR7 reads misaligned to this region.

Although it was not possible to make direct comparisons to all birds with TLR data, our findings are consistent with the lower TLR gene diversity found within small, isolated populations of other threatened species (African penguins (Spheniscus demersus); Dalton et al., 2016; Mohua (Mohoua ochrocephala) Grueber et al., 2015; Orange‐bellied parrot (Neophema chrysogaster); Morrison et al., 2020).

Comparing captive and wild tchūriwat’/tūturuatu reveals that both populations have low TLR gene diversity. These findings are consistent with Dalton et al. (2016) which show that both captive and wild populations of African penguins (Spheniscus demersus) have low levels of TLR SNP diversity. Unlike tchūriwat’/tūturuatu, captive African penguin populations have a lower number of nonsynonymous SNPs compared to wild populations (Dalton et al., 2016). However, it is also important to note that there is a difference in sample size between the captive and wild tchūriwat’/tūturuatu populations, because the captive population (n = 38) has more samples than the wild population (n = 26). The inclusion of more wild individuals may reveal additional TLR diversity, especially if there are rare TLR SNPs or haplotypes in the wild population. Both captive and wild populations have nonsynonymous SNPs that mostly fall within the LRR binding region, which may have an effect on binding affinity and TLR response (Keestra et al., 2008; Matsushima et al., 2007). Additionally, some of the nonsynonymous SNPs do cause physiochemical changes in the amino acid sequence, but further research is needed to determine how these SNPs may influence functionality of TLR proteins.

In captive and wild populations, genetic diversity is unevenly distributed across TLR genes. Most TLR genes have a low number of SNPs and haplotypes. SNP diversity within TLR7 was the highest of all TLR genes. However, this result may not reflect the true diversity of TLR7 within each population, if there is a TLR7 duplication that is not well resolved within the reference genome. Both captive and wild tchūriwat’/tūturuatu populations also had two TLRs (TLR5 and TLR7) with a value of Tajima's D that was significant and positive, indicating that rare alleles are scarce. The positive, significant values may reflect a rapid decrease in population size and/or balancing selection (Tajima, 1989). Since this result was the same for both wild and captive populations, it may reflect the overall history of species decline, given that the size and distribution of the species decreased dramatically in a short amount of time after mammalian predators were introduced (Davis, 1994) and has remained small and restricted for more than 100 years (DOC SPRG). Nevertheless, it is also possible that this result reflects balancing selection on one or both TLRs (but see below).

The finding that captive and wild tchūriwat/tūturuatu fall into two relatively distinct genetic clusters based on TLR haplotype diversity is consistent with a PCA based on 50 K SNPs (I. Cubrinovska, unpublished data). While most alleles are shared between captive and wild populations, the frequencies of these alleles differ between each population, and there are a small number of private alleles in both populations. Given that there has been little to no genetic exchange between the populations for decades, these differences can most likely be attributed to genetic drift. Previous research in other bird species has also found genetic differentiation in both neutral (microsatellites) and TLR loci, which generally indicates that genetic drift may be stronger than selection (Gonzalez‐Quevedo et al., 2015; Grueber et al., 2013; Knafler et al., 2016). However, it is possible that selection at TLR loci may contribute to the genetic differentiation between populations (Grueber et al., 2013; Knafler et al., 2016, 2017). Further analysis is necessary to determine the impacts of genetic drift and selection on TLR gene diversity in tchūriwat/tūturuatu. In the meantime, low overall TLR gene diversity combined with the genetic differentiation of captive and wild tchūriwat/tūturuatu lends support to the recent decision to augment the captive population with birds from the wild. Future research will investigate how this augmentation affects TLR gene diversity in the captive population and whether it leads to increased immune response in captive individuals.

CONCLUDING REMARKS

Below we highlight a few technical considerations that will influence the broad applicability of our bioinformatic approach to other bird species. If the species reference genome is relatively incomplete, it may not be sufficient for identifying all TLR genes. Undertaking strategies to improve completeness of the genome may help, including alignment of multiple assemblies to fill in sequencing gaps. It may also be necessary to undertake deliberate sequencing strategies to maximize depth and minimize missingness in resequencing data. For example, if preliminary screening of resequencing data shows low depth and high missingness, it may be necessary to invest in resequencing of relevant libraries. Also, if there are individuals who are particularly important, it may be beneficial to (re)sequence them to greater depth. Further, we acknowledge that our approach was used to identify TLR genes and characterize TLR gene diversity in a threatened bird species, so it is unknown whether it may be useful in other species with more complex genomes or heterozygous populations. Nevertheless, advancing technologies such as long‐read sequencing may provide opportunities to apply our approach to species with different levels of heterozygosity and genome complexity (Bayer et al., 2020; Pollard et al., 2018; Tettelin et al., 2005). Despite these caveats, we reiterate that if a whole‐genome resequencing population data set is readily available for a species of interest, these genomic resources can be used to provide a first look at TLR genes and TLR gene diversity. Further, even if TLR gene identification and characterization is incomplete or wholly unsuccessful, a bioinformatic approach is likely to produce useful information about TLR regions that can inform targeted sequencing of TLR genes (e.g., even partial information about TLR sequences can inform primer design).

CONFLICT OF INTEREST

Author declare no conflicts of interest.

BENEFIT‐SHARING STATEMENT

Tchūriwat’/tūturuatu/shore plover is an endemic shorebird classified as Nationally Critical according to the New Zealand Threat Classification System criteria (Robertson et al., 2021). This cocreated research is part of a long‐standing research partnership between the Conservation Systematics and Evolution Research Team (ConSERT, including TES, MM, IC) and Aotearoa New Zealand's Department of Conservation Shore Plover Recovery Group, which is a team of conservation scientists and practitioners (including DH, BG, TES) that provide guidance regarding the development and execution of the Shore Plover Recovery Plan to ensure the recovery of tchūriwat’/tūturuatu in the wild. This research was also conducted in collaboration with Hokotehi Moriori Trust, Ngāti Mūtunga o Wharekauri Iwi Trust, and the local community on Rēkohu/Wharekauri in the Chatham Island archipelago. For example, following engagement with both trusts regarding the selection of an individual bird for the tchūriwat’/tūturuatu reference genome, IC and Stephanie Galla visited the two schools on Rēkohu/Wharekauri: They first visited Kaingaroa School (one class, Years 1–8) and asked students to name both candidate birds for the reference genome. Then they visited Te One School (three classes, Years 1–2, Years 3–5, Years 6–8), and asked students to vote on whether “Kina” or “Maui” should be used for the reference genome, and Maui won. MM also developed an immune system outreach module that highlights our ongoing research with tchūriwat’/tūturuatu that TES subsequently shared with imi and iwi, and the principals of both schools during a recent visit to Rēkohu/Wharekauri.

Figure S1

Click here for additional data file.