A Metagenomic Investigation into Apteryx rowi Dermatosis Identifies Multiple Novel Viruses and a Highly Abundant Nematode Open Access

Ōkārito Rowi (Apteryx rowi) are one of five species of kiwi forming the family Apterygidae (Fig. 1a). Kiwi are small ratites: flightless birds with a flat, raft-like sternum (Harshman et al. 2008; Smith et al. 2012). Typically, ratites are large birds with acute eyesight, with members including Ostrich (Struthio), Cassowaries (Casuarius), and the extinct Moa (Dinornithiformes; Orrico and Gonzalez 2022). Apteryx spp. are the smallest members of this group at 0.3–0.6 m tall and, unlike other ratites, are nocturnal and rely on their sense of smell because of poor eyesight (Kummrow 2015). Kiwi are endemic to New Zealand and live a reclusive lifestyle in rainforests with high humidity. Each species occupies a different geographical region of New Zealand, with the geographical distance between populations and the georestrictive constraints on movement preventing intraspecies interactions. All species are subject to considerable conservation efforts, with four of the five species classified as nationally threatened. The wild Rowi population is distributed around the Ōkārito region of the South Island of New Zealand (Fig. 1b).

Crèches, whether on the mainland or remote off-shore islands, are predator-free conservation sanctuaries that are used to stabilize and rescue endangered species. In 2013, half of the juvenile Rowi from a crèche island were reported to have dermatitis of unknown etiology (Gartrell et al. 2015). The hematology and microbiology findings from this outbreak included leucocytosis and the detection of several potential species of nematode. This was supported by the histological identification of cutaneous larval migrans, fecal floats, and 18S ribosomal RNA (rRNA) for the Trichostrongylus genus. French et al. (2020) classified the adult and eggs in skin biopsy samples as consistent with Capillaria sensu lato, a group of related species within the family Capillaridae, and PCR confirmed the genus as Eucoleus. With reproductively active adults and suspected larval forms found within the histology of skin biopsies, French et al. (2020) suggested cutaneous capillariasis as the probable diagnosis.

A parallel study identified a full-length circo-like virus sequence in Rowi feces (White et al. 2016). Circoviruses are abundant in environmental samples and found in a wide range of animals. Avian circoviruses have been associated with immunosuppression and feather abnormalities, including psittacine beak and feather disease that is prevalent throughout Australasia. The significance of circoviruses in Rowi feces is unknown; necropsy examinations into the deaths of juvenile Rowi in 2013 suggested that they died from unrelated factors, including aspergillosis and pulmonary complications linked to spore-contaminated tree bark (Gartrell et al. 2015).

In 2017, a follow-up investigation was conducted on four wild Rowi that developed dermatitis with features seen in 2013. The clinical signs classified as lesions located at the apteria that extended outwards from this location. Most cases examined had chronic changes to the skin, described as thickening and lichenification of the affected area, along with feather loss. Histology from skin biopsies showed mild to moderate epidermal hyperplasia and hyperkeratosis. Because Rowi are captured once annually, there is limited opportunity to capture animals in an acute stage of infection, when identification of potential causative pathogens and epidemiological parameters would be more likely.

In 2023, a typical presentation of dermatosis was detected in wild Rowi. We employed a metatranscriptomic approach to identify the viruses, bacteria, fungi, and other parasitic species present and assess their possible association with this dermatosis. To this end, we collected samples from Rowi suffering from mild, moderate, and severe disease, as well as healthy animals (Fig. 1c).

A total of 34 Rowi were tracked using radiotransmitters and captured by hand for health checks. During a full physical examination and replacement of radio transmitters, skin abnormalities were noted and photographed. A severity scale pictorial was used to score and grade the dermatitis (grade: healthy, active, or chronic; and score: mild, moderate, or severe) based on clinical signs. Both oral and cloacal swabs were collected from each Rowi during the examination and stored in RNAlater (Thermo Fisher Scientific, Waltham, Massachusetts, USA). Full-thickness skin biopsies were collected from Rowi with dermatitis, including all individuals classified as severe (n=3) and moderate (n=1) plus 8/16 of those classified as mild. All samples were sent to the University of Otago, Dunedin, New Zealand, and stored at −80 C before RNA extraction.

Swabs were placed in ZR BashingBead lysis tubes (Zymo Research, Irvine, California, USA) filled with 1 mL of stabilizing reagent (DNA/RNA Shield, Zymo Research). Lysis tubes were homogenized for 5 min using a mini-beadbeater 24 disruptor (BioSpec Products, Bartlesville, Oklahoma, USA). The homogenized samples were lysed, and RNA purified using the ZymoBIOMICS MagBead RNA kit (Zymo Research). Residual guanidine contamination was removed using three short washes with molecular grade ethanol. The final concentrations of purified RNA were quantified using a NanoDrop spectrophotometer (Thermo Fisher Scientific).

We performed total RNA purification from skin biopsies using the RNeasy Plus Mini (Qiagen, Hilden, Germany) kit. Frozen skin biopsies were partially thawed and transferred into a round bottom tube filled with Buffer RLT Plus (Qiagen) that contained 1% β-mercaptoethanol and 0.5% (v/v) Reagent DX (Qiagen). Any crystals that had formed on the sample were removed before transferring. The samples were homogenized for 1 min using a TissueRuptor (Qiagen). The homogenate was then centrifuged for 5 min at 12,000 × G, ensuring that all potential tissue residue was removed. The supernatant was transferred into a gDNA eliminator spin column and purified as per the manufacturer’s protocol. An additional 500 μL ethanol wash and 1 min spin at ≥8,000 × G (≥10,000 rpm) was introduced before the final RNA elution step to reduce guanidine contamination. Purified RNA was quantified using a NanoDrop (Thermo Fisher Scientific).

Total RNA from individuals that had all three sample types (i.e., oral swab, cloacal swab, and skin biopsy) was kept separate (Fig. 2). Total RNA from healthy Rowi samples (n=14) was pooled by equal volumes into four groups: two pools from oral swabs and two pools from cloacal swabs. The RNA from mild diseased Rowi samples that did not have a corresponding skin biopsy (n=8) was pooled into four groups similarly. Library preparation for total RNA sequencing was performed using the Stranded Total RNA Prep, Ligation with Ribo-Zero Plus kit (Illumina, San Diego, California, USA) using 16 cycles of PCR. Prepared RNA libraries underwent paired-end 150-bp sequencing on a NovaSeq X (Illumina) sequencing platform.

Each paired-end library was initially screened using FastQC v0.11.9 (Andrews et al. 2010). Paired-end reads were de novo assembled using Trinity v2.14.0 (Grabherr et al. 2011). The trimmomatic flag was set to use the built-in trimmomatic plugin v0.36 (Bolger et al. 2014) with default settings for quality trimming of reads and the removal of sequencing adapters. RNA-Seq by Expectation-Maximization v1.3.3 (Li and Dewey 2011) was used to execute an alignment-based strategy to estimate transcript-level abundance.

Sequence similarity searches to annotate transcripts were performed against the NCBI nucleotide (nt) database and the nonredundant (nr) protein database with BLASTn v2.13.0 (Camacho et al. 2009) and DIAMOND (BLASTx) v2.1.6 (Buchfink et al. 2021), respectively. A maximum expected value (e-value) threshold of 1 × 10⁻¹⁰ was used as a significance cutoff. Non-viral BLAST hits with sequence similarity to host (all Apteryx spp.) or invertebrate host-associated viruses were removed from downstream analysis. Geneious Prime v2023.2.1 (Biomatters, Auckland, New Zealand) was used to identify open reading frames (ORFs). We estimated the abundance of those viral transcripts that had high confidence of being an avian virus based on their closest genetic relatives. Abundance calculations were summarized for each viral family per library and standardized by dividing by the total number of paired-end reads in a given library and multiplying by 1 million, providing reads per million (RPM). To eliminate low-abundance transcripts that might arise due to index hopping, we set a minimum cutoff of 0.01% based on a transcript’s standardized abundance in a given library compared to its highest abundance.

Amino acid polymerase sequences representing the known diversity of each viral family were obtained from GenBank (Supplementary Material Table S2) and aligned with translated viral transcripts found here using MAFFT with the E-INS-i algorithm and BLOSUM62 scoring matrix (Katoh and Standley 2013). Poorly aligned regions were trimmed using trimAl v1.4.rev22 (Capella-Gutiérrez et al. 2009) and sites that contained more than 80% gaps were removed. Phylogenetic trees for each viral family were estimated using IQ-Tree v1.6.12 (Nguyen et al. 2015) using the Le and Gascuel (LG) amino acid substitution model (Le and Gascuel 2008) with branch support assessed using ultrafast bootstrap (UFBoot; Hoang et al. 2018) and Shimodaira-Hasegawa approximate likelihood ratio test (SH-aLRT; Guindon et al. 2010) with 1,000 replicates. The graphical viewer FigTree v1.4.4 (Rambaut 2006) was used to visualize and edit phylogenetic trees.

A sliding-window sequence similarity plot was used to analyze per-site window sequence similarity of the Rowi hepacivirus polyprotein and its two closest related viruses, Bald eagle hepacivirus (QFR04963) and Hepacivirus Q (USF20790). The amino acid sequences were used as input with SimPlot++ v1.3 (Samson et al. 2022). Using a simple percentage model, we generated a similarity plot with a sliding window of 200 amino acids.

The read alignment tool Bowtie2 v2.5.3 (Langmead and Salzberg 2012) was used to index the A. rowi reference genome (RefSeq GCF_003343035) and remove host sequences from the nonassembled raw reads. We used the metagenomic classification pipeline CCMetagen v1.4.0 (Marcelino et al. 2020) for the metagenomic classification of nonviral sequences, and k-mer alignment (KMA) v1.4.14 (Clausen et al. 2018) for the initial read mapping against the pre-indexed CCMetagen reference database, applying a taxonomic filtering threshold at the genus level.

Beta diversity was explored using nonmetric multidimensional scaling. The abundance values for each taxonomic classification were normalized to RPM. The vegdist function in R v4.3.2 (Oksanen et al. 2022; R Core Team 2023) was used to calculate the dissimilarity matrices based on the Bray-Curtis dissimilarity metric. The metaMDS function was used to transform the dissimilarity matrix into a lower-dimensional representation and visualized with ggplot2 v3.4.4 (Wickham 2016). Dissimilarity was plotted using the groupings of sample type and health status. The adonis2 function in R v4.3.2 was applied to perform the permutational multivariate analysis of variance (PERMANOVA) test on the distance matrix to identify categorical variables that were significantly influencing the multivariate composition of the community. All taxonomic classifications were ranked based on their total abundance. The classifications ranked within the top 20 genera were selected for further analysis.

A full-length transcript identified as Eucoleus 18S rRNA and a partial cytochrome c subunit I (COX1) sequence were aligned with those from the Capillariidae obtained from GenBank (Supplementary Material Table S4). Maximum-likelihood phylogenetic trees were estimated for both genes using IQ-TREE, employing the general time-reversible model of nucleotide substitution with a gamma-distributed rate variation and an invariant site category (GTR+I+Γ). Branch support was evaluated using UFBoot (1,000 replicates) and SH-aLRT (1,000 replicates).

In total, 76 samples were successfully obtained from 34 Rowi, consisting of 32 cloacal swabs, 32 oral swabs, and 12 skin biopsies. At the time of sampling, 20 diseased Rowi were identified. All but one of the diseased Rowi were diagnosed with chronic disease; the single active infection was classified as mild. A total of 42 libraries were subject to RNA sequencing, comprising nine severe, three moderate, 26 mild, and four pooled libraries from mild and healthy control samples (Fig. 2a and Supplementary Material Table S1). The sequencing yield from these libraries ranged from 14 to 158 million reads, with a mean of 57 million reads (Fig. 2a and Supplementary Material Table S1).

Across 13 libraries, we identified RNA transcripts from three families of RNA viruses— the Coronaviridae, Picornaviridae, and Flaviviridae—that shared sequence similarity to viruses that infected avian hosts (Fig. 3 and Supplementary Material Table S1). All but one of the identified viral families were associated with mild cases of dermatitis, with only a single library from pooled healthy cloacal swabs having a very low abundance of transcripts belonging to the Coronaviridae. None of the libraries were below the cutoff threshold set for index hopping. Viral transcripts were identified across all sample types, with oral swabs representing the largest portion (54% of oral samples).

A total of 33 viral transcripts with sequence similarity to the Flaviviridae were identified. The closest genetic relatives included Bald eagle hepacivirus (QFR04963) with 42% amino acid identity, Hepacivirus Q (USF20790) with 41% amino acid identity, and Jogalong virus (QHD25538) with 38% amino acid identity across the polyprotein (Fig. 4a). These transcripts, including several full-length and partial polyproteins, were the most prevalent viral transcripts by both the number of libraries (24%) and the range of standardized abundance (0.04–14.5 RPM). The virus genome, provisionally named Rowi hepacivirus, was only detected in mild cases of dermatitis across all sample types (cloacal swab, oral swab, and skin biopsies, and in both chronic and active infections). This virus was absent from clocal and oral swabs from healthy individuals (skin samples were not obtained from healthy individuals). Phylogenetic analysis placed Rowi hepacivirus within the genus Hepacivirus in a clade containing avian hepaciviruses, including those that infect Zebra finch (Taeniopygia spp.) and ducks (Fig. 4a).

A complete genome for Rowi hepacivirus was recovered that exhibited 39.4% amino acid sequence similarity to Bald eagle hepacivirus. The 9,289 nucleotide sequence contained a single ORF that translated to a 3,064 amino acid polyprotein. A sliding-window sequence similarity analysis comparing Rowi hepacivirus and reference sequences from Bald eagle hepacivirus and Hepacivirus Q revealed considerable divergence in the novel hepacivirus (Fig. 4d). Specifically, compared to Bald eagle hepacivirus, the window of similarity ranged from 32% to 62%, being most closely related in the RNA-dependent RNA polymerase (RdRp). Similarly, Hepacivirus Q displayed a similar range of 33–59%.

Two transcripts shared sequence similarity with the ORF 1ab of White-eye coronavirus HKU16 (AFD29200) and Gentoo penguin virus (QIS87929) from the genus Deltacoronavirus within the Coronaviridae family, with 81% and 89% amino acid identity, respectively. One transcript was from an oral swab collected from a mild case (standardized viral abundance 0.1 RPM) and the other from a pool of healthy Rowi cloacal samples (standardized viral abundance 0.03 RPM). These transcripts, termed Rowi coronavirus 1 and 2, were related to other avian deltacoronaviruses (Fig. 4b). Deltacoronaviruses are one of two main avian-associated genera of coronaviruses (Woo et al. 2012; Lau et al. 2018; Miłek and Blicharz-Domańska 2018). These viral transcripts did not overlap in the alignment such that they may represent the same viral species across multiple libraries.

Sequence similarity analysis of a single transcript showed 43% amino acid sequence similarity to the polyprotein of Pigeon picornavirus B (AGJ03164), a member of the genus Sapelovirus in the family Picornaviridae. This transcript was identified in an oral swab from a mild case, with low abundance of 0.14 RPM. Tentatively termed Rowi anativirus, phylogenetic analysis placed this virus within the Anativirus genus (Fig. 4c)—this genus, containing avian-associated viruses, was recently split from the porcine-associated viruses in the Sapelovirus genus (Sunaga et al. 2019). Rowi anativirus clustered with sapeloviruses (anativiruses) associated with mute swans (Cygnus olor) and an anativirus-like virus from ducks (Anas superciliosa).

Composition of the nonviral microbiome revealed significant clustering based on bird health status and sample type (Fig. 5a). PERMANOVA analysis supported the identification of clusters associated with sample type (P=0.001) and health status (P=0.006). Oral swabs from mild cases formed a tight cluster, indicating high similarity in both presence and abundance of nonviral microbial organisms. Healthy and moderate cases from the same sample type aligned around or within clusters and maintained proximity to the tight clustering from mild cases. Overall, these results suggested that both sample type and health status were most strongly associated with the observed dissimilarity among samples.

Compositional analysis was conducted on the 20 most abundant nonviral microbial taxa (Fig. 5b). The analysis revealed a diverse range of organisms, including bacteria, fungi, archaea, protozoa, bdelloidae, and nematoda. Most identified taxa were those commonly found in environmental or plant-associated samples. Notably, several bacteria, such as those within the Bacteroidota phylum (formerly known as Bacteroidetes), are typical constituents of the gut microbiome. Malassezia, a known skin-associated fungus in avian hosts, was also prevalent across all sample types. Although its role in pathogenesis remains uncertain, it is believed to interact with the avian immune system (Preziosi et al. 2006; Eidi et al. 2016; Benmaarouf et al. 2023). Taxa relevant to pathogenesis in avian hosts, such as Mycoplasmopsis and Eucoleus spp., were of interest. Mycoplasmopsis synoviae has been associated with subclinical upper respiratory infections in chickens (synovitis); Eucoleus spp. nematodes are known to infect birds (Stenkat et al. 2013; French et al. 2020; Rehman et al. 2023). It is worth noting that although we also detected transcripts from Plasmodium spp., these were in low abundance and across all health statuses.

Comparison of taxonomic compositions across libraries based on sample type and health status revealed a shift in the abundance of Eucoleus transcripts, particularly in severe dermatitis cases. Eucoleus was detected in eight libraries, predominantly from severe dermatitis cases across all sample types. In these cases, Eucoleus constituted a significant proportion of the total abundance, ranging from 82% to 99.8% of the non-Rowi RNA (Supplementary Material Table S3). Only two libraries from severe cases did not contain Eucoleus transcripts.

Severe cases had the greatest combined standardized abundance (2,536 RPM) compared to other health statuses (healthy: 497 RPM, mild: 102 RPM, moderate: 0 RPM). The combined abundance of Eucoleus in severe cases was significantly larger than in mild (P=0.038, by Welch’s t-test) and moderate (P=0.036, by Welch’s t-test) cases (Fig. 6a).

Transcripts assembled from 10 libraries had sequence similarity to Eucoleus aerophilus and Eucoleus contortus. An assembled Eucoleus 18S rRNA consensus sequence clustered within the Eucoleus genus, sharing a common ancestor with Eucoleus perforans (Fig. 6b). The reference sequences used in this study primarily consisted of partial 18S rRNA sequences, including those from a species of Capillaria previously reported to have infected Rowi (French et al. 2020). The alignment indicated significant divergence between the Eucoleus species identified in our study and the short (270–279 nucleotides) 18S rRNA sequences obtained from prior studies. Our phylogenetic analysis supported the placement of the newly identified 18S rRNA sequence within the genus Eucoleus (Fig. 6b). The COX1 consensus sequence was phylogenetically placed with short sequences of Eucoleus previously reported to have infected Rowi (French et al. 2020; Fig. 6c).

Raw reads are available through the Aotearoa Genomics Data Repository [PENDING]. Viral transcripts are available on GenBank under accessions PP946257–PP946260. Nematode transcripts are available on GenBank under accession numbers PP951251 and PP951580–PP951585. For further analysis, the alignment files and code used for statistical analysis are on our GitHub repository (https://github.com/Ditean/rowi-dermatosis).

We employed total RNA sequencing to characterize the infectome of wild Rowi suffering from dermatosis, an approach that has played a crucial role in the a priori analysis of the transcriptome for pathogen discovery (García et al. 2020; Li et al. 2020; Rajagopala et al. 2021; Mick et al. 2023). One significant advantage of targeting RNA is its ability to facilitate the sequencing and characterization of RNA viruses, which might be missed using DNA-based approaches (Arroyo Mühr et al. 2021; Du et al. 2022).

Of the three virus families, Coronaviridae, Picornaviridae, and Flaviviridae, from which we detected viruses with sequence similarity to known avian viruses, Rowi anativirus and Rowi coronavirus 1 and 2 were detected in low relative abundance and were closely related to known avian viruses. Members of the genus Sapelovirus (Picornaviridae) have a broad vertebrate host range and could be considered ageneralists, with viruses including simian, avian (anativirus), and porcine enterolike (sapelovirus) viruses (Zell et al. 2017). Within the Coronaviridae, Deltacoronaviruses are largely associated with avian hosts, with pigs being the only nonavian host identified to date (Ye et al. 2020; Woo et al. 2023). The low viral abundance of Rowi anativirus and Rowi coronavirus 1 and 2, coupled with their infrequent detection across samples, suggested that these viruses are unlikely to play a role regarding the described dermatosis in wild Rowi.

Hepaciviruses (Flaviviridae) are nonsegmented positive-sense RNA viruses usually associated with liver tropism (Simmonds et al. 2017). These viruses have a broad host range including birds (Goldberg et al. 2019; Zhang et al. 2022), mammals (Tabata et al. 2020), reptiles, and amphibians (Harding et al. 2022). Waterfowl and Bald Eagles (Haliaeetus leucocephalus) are recent expansions of the known host range of hepaciviruses (Chu et al. 2019; Goldberg et al. 2019). The Rowi hepacivirus that we detected was most closely related to Bald eagle hepacivirus and Hepacivirus Q. Birds infected with Bald eagle hepacivirus showed acute clinical signs, with birds presenting with severe neurological deficit and liver disease; there was active viral replication in liver (Goldberg et al. 2019). Hepacivirus Q has been detected in ducks (Anatidae) across China, and is speculated to be involved with a decline in egg production (Zhang et al. 2022). The Rowi hepacivirus was both highly abundant and prevalent across many sample types—both oral and cloacal swabs as well as in skin biopsies, and in both acute and chronic infections. Nevertheless, although Rowi hepacivirus was absent from healthy individuals, it was only found in about 50% of mild cases of disease, was absent from more severe dermatitis, and hepaciviruses have not previously been associated with dermatitis. The relationship between Rowi hepacivirus and dermatitis is therefore uncertain, although it merits additional investigation.

Eucoleus transcripts were dominant in severe cases of chronic dermatitis. Although Eucoleus was absent from most samples, in samples in which it was present, Eucoleus exhibited remarkably high abundance, >90% of the total relative abundance of the nonviral microbiome. French et al. (2020) speculated that Eucoleus was probably not the primary cause of the dermatitis, rather suggesting its involvement in perpetuating the disease. Whether Eucoleus represents an opportunistic parasite during severe cases of dermatitis is yet to be determined. Based on the phylogenetic placement of the COX1 consensus, the Eucoleus species that we found is probably the same as those previously identified by French et al. (2020).

Disease investigations of endangered wildlife species are highly challenging. Because of limited, once-yearly sampling, we know little about the epidemiology of dermatosis in Rowi. The findings in our study are preliminary, but provide a basis for further investigation into the potential causes of dermatosis in Rowi, with both the Rowi hepacivirus and the probable Eucoleus sp. nematode requiring follow-up investigations.

We are grateful to Te Rūnaka o Makaawhio as kaitiaki (guardians) of Rowi and for their support of this research. The authors acknowledge the use of New Zealand eScience Infrastructure (NeSI) high-performance computing facilities, consulting support, and/or training services as part of this research. New Zealand's national facilities are provided by NeSI and funded jointly by NeSI's collaborator institutions and through the Ministry of Business, Innovation & Employment's Research Infrastructure program (https://www.nesi.org.nz). J.T.T. was funded by a University of Otago Pacific Island Master’s Research Scholarship. J.L.G. is funded by a New Zealand Royal Society Rutherford Discovery Fellowship (RDF-20-UOO-007), a Marsden Fund Fast Start (20-UOO-105) and the Webster Family Chair. J.W. is funded through a Morris Animal Foundation Fellowship. E.C.H. is funded by a National Health and Medical Research Council (Australia) Investigator Grant (GNT2017197).

Supplementary material for this article is online at http://dx.doi.org/10.7589/JWD-D-24-00115.