The identification of immune genes in the milk transcriptome of the Tasmanian devil (Sarcophilus harrisii)

Ethics statement

The milk collection in this study was approved by The University of Sydney Animal Ethics Committee (Animal Ethics no. 6039).

Sampling

Approximately 10 mL of milk was collected from a Tasmanian devil at approximately four months (∾120 days) post-parturition. The animal was held at the Australian Reptile Park, Somersby, NSW, Australia.

RNA isolation and sequencing

The milk sample was kept on ice during transport and RNA extraction was carried out within 2–3 h after collection. Approximately 7 ml of milk was centrifuged at 2,000 × g at 4 °C for 10 min. The top fat layer was removed and the bottom layer was washed once using 10 ml PBS (pH 7.2, with 0.5 mM EDTA) followed by centrifugation at 2,000 × g and 4 °C for 10 min. The cell pellet (which would have contained neutrophils, macrophages and lymphocytes, with a lower abundance of epithelial cells and granulocytes (Young et al., 1997; Young & Deane, 2001) was recovered and RNA was extracted using 1 ml TRIzol^® Reagent (Life Technologies, Carlsbad, CA, USA), following the RNA extraction protocol of the manufacturer, including removal of fat, proteins and other material, and an RNase-Free DNase Set (Qiagen, Hilden, Germany) was used to remove DNA contamination from the extracted RNA. A second round of purification was performed using an RNeasy Mini Kit (Qiagen, Hilden, Germany) and quality-checked on a Bioanalyser (Agilent Technologies, Santa Clara, CA USA). The final yield of total RNA was approximately 320 ng and RNA integrity number was 7.3. Library construction and sequencing were performed by The Ramaciotti Centre (UNSW) with TruSeq chemistry on a HiSeq2000 (Illumina, San Diego, CA, USA). Approximately 22.5 million paired 100 bp reads were obtained, totalling over 44.9 Gbp of data. Reads were submitted to the NCBI Sequence Read Archive under the BioProject accession PRJNA274196and BioSample accession SUB812082. The assembled transcriptome was submitted to the Transcriptome Shotgun Assembly Sequence Database under accession GEDN00000000.

Transcriptome assembly and annotation

RNAseq reads were assembled with the Trinity pipeline (released 10 November 2013) (Haas et al., 2013) using the default parameters. This assembly resulted in 200,829 contigs, with an N50 contig size of 2,720 bp, a mean contig length of 1,285 bp, and a transcript sum of 258.1 Mb. Functional annotation of the devil milk transcriptome was performed using the Trinotate pipeline (released 10 November 2013 (default settings): http://trinotate.github.io). In brief, BLASTp was performed using devil milk predicted ORFs as the query and the SwissProt non-redundant database (accessed 29th July 2013) as the target and the de novo transcripts aligned against the same using BLASTx (Altschul et al., 1990). HMMER v3.1b1 and Pfam v27 databases (Finn, Clements & Eddy, 2011) were used to predict protein domains, SignalP v4.1 (Petersen et al., 2011) to predict the presence of signal peptides, RNAmmer v1.2 (Lagesen et al., 2007) to predict rDNA, and TMHMM v2.0 (Krogh et al., 2001) to predict transmembrane helices within the predicted ORFs from the milk transcriptome. These transcriptome annotations were loaded into an SQLite database, and abundance estimation was performed using the RSEM v1.2.1 (default settings) (Li & Dewey, 2011) method. GO terms linked to the SwissProt entry of the best BLAST hit were used for ontology annotation. GO functional classifications and plotting was performed by WEGO (http://wego.genomics.org.cn) (Jia et al., 2006).

Top 200 most highly expressed transcripts

The most highly expressed transcript variants were selected based on FPKM (fragments per kilobase of exon per million fragments mapped) gene expression estimation. Twenty-four transcripts of the top 200 had no BLAST hits through the Trinotate pipeline. In order to annotate these transcripts, they were further investigated using BLAST searches against the Tasmanian devil genome on ENSEMBL (release 75) (http://www.ensembl.org/index.html), or tBLASTx to GenBank nucleotide and EST collections. Transcripts that had poor BLAST hits (E-value >1 × 10⁻¹⁰) to SwissProt sequences were also verified using these methods.

For sequences that could not be identified using the above methods, HMMER v3.1b1 (Finn, Clements & Eddy, 2011) and SignalP v4.1 (Petersen et al., 2011) searches were used. Rfam 12.0 (Griffiths-Jones et al., 2003) and Pfam v27 (Finn et al., 2014) were used to identify conserved RNA and protein domains respectively. Finally, genes were identified with the aid of conserved flanking genes in the tammar wallaby and gray short-tailed opossum (Monodelphis domestica) genomes. Genes flanking the unidentified genes were identified in the devil genome. Syntenic regions in the opossum and wallaby genomes were then searched for genes with homology to the devil sequence using FGENESH+ (Solovyey, 2007). Using this process, the top 200 highly expressed transcripts from the milk transcriptome were identified and annotated.

Phylogenetic analysis

To investigate the evolutionary relationship between the late lactation proteins (LLP) of the devil and the various marsupial species, LLP protein sequences from the tammar wallaby (LLP-A (Genbank: AAQ15117), LLP-B (Genbank: AAL85634)), gray short-tailed opossum (Monodelphis domestica) (LLP-B1 (Genbank: XP˙007475421), LLPB-B2 (Ensemble: ENSMODG00000017471), LLP-B3 (Ensemble: ENSMODG00000017468), LLP-B4 (Ensemble: ENSMODG00000025759)), quokka (Setonix brachyurus) (Genbank: AAB33234.1) and brushtail possum (Trichosurus vulpecula) (Genbank: AAA93179.1) were obtained for phylogenetic tree construction. Tasmanian devil protein sequences for LLP-A, LLP-B, and a homolog (named LLP-like), predicted from the devil transcripts, were used in the phylogenetic tree construction for evolutionary analysis. As the transcript for LLP-like was only partial, to obtain the full sequence the missing exons were predicted from the devil genome. The region encoding LLP-like was identified using BLAST to the devil genome. This region and devil LLP-B were used as inputs in FGENESH+ (Solovyey, 2007), to identify the missing exons for the prediction. As LLP proteins belong to the lipocalin family, the opossum lipocalin-1 (Genbank: XP˙007475462.1) was used as an outgroup. MEGA 6 (Tamura et al., 2013) was used to analyse the phylogenetic evolutionary relationships between the marsupial LLP sequences. Protein sequences were aligned using MUSCLE (Edgar, 2004) using default settings (see File S1). Using the Model Selection tool in MEGA6 the JTT model (Jones, Taylor & Thornton, 1992) was identified as the best fit and a phylogenetic tree was constructed using the maximum likelihood method based on the JTT model with 1,000 bootstrap replicates. Additional alignments were produced for novel devil genes using the ClustalW algorithm (Thompson, Higgins & Gibson, 1994) in BioEdit (Hall, 1999).

Identification of immune transcripts

A list of immune transcripts in the milk transcriptome was generated by searching the milk transcriptome with proteins from the Immunome Database for Marsupials and Monotremes (IDMM) (Wong, Papenfuss & Belov, 2011) using tBLASTn. IDMM is a database of immune genes obtained from a number of marsupials including the tammar wallaby, gray-short tail opossum, brushtail possum, northern brown bandicoot (Isoodon macrourus), and the monotremes platypus (Ornithorhynchus anatinus) and echidna (Tachyglossus aculeatus). Additionally, the milk transcriptome was searched with a range of devil specific immune genes identified from the devil genome using BLAST. This included devil cytokines, immunoglobulin constant regions (Morris et al., 2015), NK cell receptors (Van Der Kraan et al., 2013), defensins (EA Jones, 2015, unpublished data) and cathelicidins (E Peel, 2015, unpublished data). The most highly expressed transcripts were selected based on FPKM expression estimates. Transcripts that had poor BLAST hits (E-value >1 × 10⁻¹⁰) to marsupial sequences in IDMM were verified using tBLASTx to the GenBank nucleotide collection or the devil genome on ENSEMBL.

Transcriptome overview

A transcriptome was constructed and annotated from the total milk cells of a Tasmanian devil milk sample obtained at the end of mid-lactation. We note that Tasmanian devils are an endangered species and access to milk samples is opportunistic. In this case, a female devil was given veterinary treatment due to an injury and milk could be collected while she was anaesthetised. The total number of transcripts, including transcript variants, that were expressed in the milk transcriptome was 233,660. Excluding transcript variants, the number of transcripts was 101,399, and of these, 17,827 sequences had BLAST hits to the SwissProt non-redundant database. This number provides an estimate of the number of protein-coding genes within the Tasmanian devil milk transcriptome. Transcripts that were not protein-coding genes included non-coding RNAs, transposons, and retroelements. The transcriptome included 845 immune genes representing 4.7% of the protein-coding genes and accounted for 6.6% of the transcript expression in the transcriptome.

Of the 17,827 sequences with BLAST hits, GO terms were assigned to 16,437 transcripts. A total of 51 level 2 GO terms were assigned (Fig. 1) (see Table S2). Within the molecular function category, binding (60.6%) and catalytic activity (34.4%) were the most common functions. Of the biological processes, genes categorised as being involved in cellular processes (75.4%) were the most common, followed by metabolic processes (53.2%), biological regulation (40.0%), and developmental processes (27.0%).

A large number of transcripts are also classified as having a role in immune system processes (8.2% or 1,341 transcripts). This number is slightly larger than the number of immune genes identified by BLAST to the IDMM database; this is likely due to a broader range of genes being classified as having a role in immune system processes within the GO annotation. Within the transcripts in the immune system process category, 14 sub-categories were represented (Fig. 2) (see Table S3). Of these transcripts, immune response genes were most common (55.5%), followed by immune system development (39.0%), leukocyte activation (37.5%) and immune effector process (24.9%).

Most highly expressed transcripts

The 200 most highly expressed transcripts in the Tasmanian devil milk transcriptome encode a range of nutritional milk proteins and immune proteins (Fig. 3) (see Table S4). A large proportion of most highly expressed transcripts were milk protein transcripts common across most mammals, including α-, β-, and κ-caseins, α-lactalbumin, β-lactoglobulin and whey acidic protein (WAP), as well as the marsupial-specific milk proteins, early lactation protein (ELP) and late lactation proteins (LLP-A and LLP-B). These together accounted for 50.21% of the total gene expression and made up the majority of the ten most highly expressed transcripts (Table 1). These proteins are mostly involved in nutrition, providing amino acids and minerals to the young, although some also have potential immune roles (Oftedal, 2012). Transcripts for proteins associated with energy metabolism, such as nicotinamide adenine dinucleotide (NADH) dehydrogenase, adenosine triphosphate synthase lipid-binding protein, and cytochrome oxidases, are present in the top 200 highly expressed transcripts.

The highest expressed transcript was ELP, accounting for 21.93% of all transcripts. LLP-A and LLP-B were also highly expressed, as the fourth and sixth most highly expressed transcripts respectively, representing 8.66% and 3.91% of the total gene expression respectively (Table 1). The transcriptome contained a third transcript with homology to LLP-A and -B, named LLP-like, which may represent a novel Tasmanian devil LLP. To investigate the relationships between LLP-like and other Tasmanian devil and marsupial LLPs, a phylogenetic tree was constructed (Fig. 4). The number of LLP genes appears to differ between species; while a single LLP sequence has been identified in quokka and brushtail possum, two have been identified in wallaby. Through BLAST searches to the opossum genome, four opossum LLP homologs were identified in the opossum. The evolutionary relationship between Tasmanian devil LLP-like with other marsupial LLP sequences could not be definitively resolved due to weak bootstrap support. However, LLP-like does group with the other marsupial LLP sequences, and shares substantial amino acid sequence identity with other marsupial LLP sequences (23.2–44.6%), thus this sequence is likely to represent an additional LLP locus in the devil.

Novel transcripts

Two novel gene transcripts were identified in the top 200, ranking 27th and 58th and accounting for 0.26% and 0.02% of the total gene expression respectively (see Table S5). Conserved motifs could not be identified in either gene. The first novel gene, which we designate here as novel gene 1, contained three exons that aligned to the Tasmanian devil genome, however there was no gene prediction made by the ENSEMBL annotation in this region. It did contain a predicted signal peptide cleavage site and a polyA tail, suggesting that it encodes a protein. A putative ortholog to novel gene 1 was identified in the tammar wallaby mammary gland transcriptome (Genbank: EX196900.1). The nucleotide and protein alignments of the devil and wallaby sequences are shown in Figs. 5 and 6. The two sequences have 64% and 68% nucleotide and amino acid identity respectively. Given that a homolog could not be identified in eutherians or non-mammals and that transcripts of this gene could only be identified in milk or mammary transcriptomes, we propose that novel gene 1 may have a marsupial-specific role in lactation.

The second novel gene did not contain a signal peptide cleavage site. There were no gene predictions in the region encoding this sequence in the devil genome, nor homologs identified through HMMER searches. Additionally, the transcript is 2,072 base pairs long but does not contain any open reading frames greater than 90 residues, thus it seems unlikely to be protein-coding. Based on the length and the lack of open reading frames we speculate that it may be a long regulatory RNA. Two possible homologous sequences were identified in the tammar testis (E value: 4 × 10¹⁶¹) and tammar uterus (E-value: 0.00) transcriptomes through BLAST. Their alignment against the Tasmanian devil nucleotide sequence is shown in Fig. 7. The devil sequence has a sequence identity of 28.9% to the tammar testis and 27.8% to the tammar uterus sequences. Although this is quite low overall, there are regions within the sequences with very high identity, for example bases 1,156 to 1,499 in the devil sequence have 86.8 and 87.4% identity to the wallaby uterus and testis sequences respectively.

Immune transcripts in the milk

There were 846 immune gene transcripts identified in the milk transcriptome, representing 6.6% of the total gene expression. The top ten are listed in Table 2, and the relative expression of the top 200 immune transcripts is shown in Fig. 8 (see Table S6). The most highly expressed immune transcripts include those encoding lysozyme C, WAP, ferritin, MHC I, S100A (calcium-binding) proteins, and CCL25.

All four isotypes of marsupial Igs were present in the devil milk transcriptome, consistent with the gene expression profile in the tammar wallaby during mid-lactation (Lefevre et al., 2007) (Table 3). The total relative expression of all Ig transcripts was 0.01%. Transcripts encoding IgM, IgG and IgE were also present in the milk transcriptome but at lower abundance, representing 9.99 × 10⁻⁴%, 4.29 × 10⁻⁴% and 2.21 × 10⁻⁵% of the total expression respectively. Two Ig light chain isotypes, Igκ and Igλ1 were identified, with Igκ being the most highly expressed (Table 3). Four Igλ light chains have been identified in the Tasmanian devil (Morris et al., 2015), however only Igλ1 was identified in the milk. Additional receptors involved in transfer and protection of Igs were also identified. β-2 microglobulin light chain (β-2), a receptor important for efficient transfer and uptake of IgG across the PY gut (Daly et al., 2007; Joss et al., 2009) was highly expressed in the devil milk transcriptome. Polymeric immunoglobulin receptor (pIgR), which protects Ig in the gut and has antimicrobial properties (Kaetzel, 2005), was also very highly expressed.

A range of immune cell receptor genes were expressed in the milk transcriptome, providing insights into the cell types that are likely to have been in the sample at the time of collection. Phagocytes, macrophages, dendritic cells, monocytes, granulocytes, T helper and cytotoxic T cells were likely to be present based on the expression of two MHC class II β transcripts and MHC II α transcripts (Ting & Trowsdale, 2002), numerous toll-like receptors (1, 4, 5, 7, 8 and 9) (Akira, Takeda & Kaisho, 2001), CD14, (Taylor et al., 2005) CD4 (Doyle & Strominger, 1987) and CD8 (Gibbings & Befus, 2009).

Ten natural killer receptor (NKR) transcripts (Table 4) and 18 chemokine gene transcripts were identified in the milk transcriptome (Table 5). CCL25, IL8, CXCL1L and CCL28 were the most highly expressed (Table 5). Antimicrobial peptides were also present (Table 6). This included four cathelicidins and three β-defensins.