USING MULTIVARIATE ANALYSES TO EXPLORE DISEASE PROGRESSION OF FINCH MYCOPLASMOSIS Open Access

Many diseases that threaten wildlife populations cause external pathology, providing a minimally invasive, cost-effective means by which to assess individual health in the field. Researchers have developed linear, multipoint scoring systems to quantify lesion severity for a range of wildlife diseases, including fibropapillomatosis (Work and Balazs 1999), finch mycoplasmosis (Roberts et al. 2001), chytridiomycosis (Obendorf 2005), devil facial tumor disease (Lachish et al. 2007), white-nose syndrome (Reichard and Kunz 2009), and snake fungal disease (McCoy et al. 2017). However, pathology can vary in myriad ways during a single infection, among individuals, and across populations (Adelman et al. 2013). Therefore, assessment of external lesions by severity scores alone risks overlooking important trends. We tested multivariate methods for quantifying external pathology to complement a widely used scoring system in a common wildlife disease model, finch mycoplasmosis, or House Finch eye disease.

Finch mycoplasmosis is an infectious disease that causes severe inflammation of the conjunctiva. It is caused by Mycoplasma gallisepticum (MG), a bacterial pathogen of poultry that was first isolated from songbirds in 1994 during a die-off of House Finches (Haemorhous mexicanus) on the east coast of the US (Ley et al. 1996; Hochachka and Dhondt 2000). Visual scoring systems were developed soon after this pathogen's emergence in finches and subsequently refined (Roberts et al. 2001; Kollias et al. 2004; Sydenstricker et al. 2006). Although earlier scales used more comprehensive criteria, most research groups now use a four-point scale from 0 (no pathology) to 3 (severe pathology), with half-scores permitting seven levels.

Even the most nuanced scales (Kollias et al. 2004) assign different weights to different aspects of pathology, resulting in scores that are neither consistent nor reproducible across studies. This proves especially problematic in surveys evaluating the host response across time or populations. We lose information by collapsing complex, multidimensional trait space (conjunctival pathology) into a numerical score. To capture these differences, we needed new tools that can work with multivariate data. We used this framework to ask targeted questions relating conjunctival pathology to stage of infection and more broadly, to host-pathogen evolution.

We assessed two different multivariate techniques with the potential to improve clinical assessment of wildlife diseases causing external pathology, using finch mycoplasmosis as a model. First, we explored the unique contributions that geometric morphometrics can offer the field of disease ecology. Geometric morphometrics is a recent addition to the evolutionary and organismal biologists' toolkit, combining the power of multivariate statistics with shape analysis (Adams et al. 2013). It was developed to capture geometric variation in the morphology of species or populations over evolutionary timescales. However, in host-pathogen systems, we can use the same methods to track rapid changes in host phenotype by sampling within single infections. Second, we applied methods from community ecology to analyze binary descriptor variables that capture the presence of distinct pathological phenotypes. Specifically, to better characterize the breadth of host pathologies and capture variation in the signs of inflammation, we analyzed the presence (1) or absence (0) of six variables: blepharospasm, crusting, erythema, eversion, exudate, and swelling.

We used these techniques to answer questions regarding variation in the grossly visible conjunctival pathology caused by finch mycoplasmosis. Specifically, we applied these techniques to paired samples from the same experimentally infected individuals to describe how eyes in early versus late infection differed in pathology. By collapsing these multidimensional data to two axes, we created a spatial representation of disease space that was defined by the specific data set, eye shape or community composition. We then asked whether occupancy of distinct regions of disease space correlated with stage of infection.

Experimental infection study: Wild-caught House Finches from Ames, Iowa, were recruited into an experimental challenge study between 23 June and 26 August 2018. Plumage patterns were used to target hatch-year birds, as these individuals are less likely to have yet experienced natural infection with MG. All birds were housed separately for a minimum 2-wk quarantine, during which eyes were evaluated every 3 d for pathology. On day 14, <75 µL of blood was drawn from the ulnar vein for serology using a commercially available enzyme-linked immunosorbent assay test kit (99-06729, IDEXX Laboratories, Inc., Westbrook, Maine, USA), which was validated by Grodio et al. (2009). Any birds with measurable MG immunoglobulin Y antibodies, that is, a sample:positive ratio >0.0229 on light spectroscopy (Hawley et al. 2011), were deemed ineligible and were released at the site of capture.

Birds that passed quarantine were assigned to one of three experimental treatment groups or to a control group evaluating the effect of a nonsteroidal anti-inflammatory drug on disease progression. Treated birds received 1 mg/kg oral meloxicam (Metacam, Boehringer Ingelheim Animal Health USA, Inc., Duluth, Georgia, USA) daily for different durations: 3, 7, or 28 d. After their course of nonsteroidal anti-inflammatory drug, birds received 2 mL/kg oral water daily, which was the same dose and frequency given to control birds for a 28-d duration. Birds were pretreated for 2 d prior to experimental inoculation with MG (strain VA1994, stock ID 7994-1-7P 2/12/09; D. H. Ley, College of Veterinary Medicine, North Carolina State University, Raleigh, North Carolina, USA). We evaluated pathology (11 times), pathogen load (5 times), and antibody titers (3 times) over 38 d.

Image acquisition and selection criteria: Images of House Finches from the experimental challenge study were collected using an iPhone 6S (Apple, Inc., Cupertino, California, USA). Size markers were affixed to stanchions projecting from the walls of the laboratory animal facility where the birds were housed. This allowed a single observer (R.M.R.) to hold a bird in bander's grip flush to the wall and take photographs, such that either the left or right eye was in the same plane as the size marker. Photographs were not taken at any magnification, and the distance from the phone to the bird did not need to be standardized because all images were later scaled by the size marker. Images were annotated at the time of collection with the bird's identification number according to the aluminum band on the left or right tarsometatarsus. Day postinoculation was verified using dates retrieved from image acquisition data. Photos were collected on days 0, 5, 7, 10, 14, 17, 21, 24, 28, 34, and 38.

At each time, a separate observer (J.S.A.) scored eye pathology on a four-point scale from 0 (absent) to 3 (severe; Sydenstricker et al. 2006). We used the day postinoculation at which an individual reached its maximum eye score to define peak pathology. We then evaluated each bird's eye score data to look for days pre- and postpeak when eyes received the same severity score. We extracted the corresponding images from our photo library to include as specimens for shape analysis.

Landmark digitization: We used the program tpsDIG2w32 (version 2.30, F. J. Rohlf, Stony Brook University, Stony Brook, New York, USA) to assign landmarks to each image using a thin-plate spline (.tps) overlay built in tpsUtil64 (version 1.74, F. J. Rohlf, Stony Brook University). We used 16 landmarks delineating the inner and outer eye margins to capture conjunctival pathology (Fig. 1). Landmarks were applied in the same sequence around each eye, with left eyes digitized as mirror images of the right, so all specimens could be superimposed using package geomorph (Adams and Otárola-Castillo 2013) in R, version 3.5.2 (R Development Core Team 2018). Individual specimens were scaled according to a 5-mm size marker visible in each raw image.

Alignment, processing, and shape analysis: We created a matrix of semilandmark sliders to define the positions of landmarks relative to their nearest neighbors. Specimens were then aligned by generalized Procrustes analysis (superimposition; Fig. 2). We evaluated for main effects of status (a binary value specifying prepeak or postpeak pathology), severity score (a categorical value from 0 to 3), and their interaction on Procrustes shape coordinates using analysis of variance on a randomized residual permutation procedure (package RRPP; Collyer and Adams 2018). Additionally, we ran principal components analysis (PCA) on Procrustes shape coordinates to assess patterns in conjunctival pathology. We used principal components (PC) 1 and 2 to collapse our multidimensional data into two dimensions to visualize disease space. The picknplot() function in geomorph then allowed us to visualize the distortion of eye shapes along each axis with thin-plate spline deformation plots. Finally, plotting eyes with multiple specimens along these axes allowed us to evaluate disease trajectories and explore how pathology progressed through disease space over time.

All eyes used in the mentioned shape analysis were also represented in the first data set we used to compare early- versus late-stage infection by community analysis of pathological descriptors. We also included paired data points from an experiment carried out in 2017 using the same source population of House Finches (Ames, Iowa) inoculated with a more virulent isolate of finch MG (strain NC2006, stock ID 2006.080-5(4P) 7/26/12; D. H. Ley, North Carolina State University). However, due to potential differences in the host-pathogen interaction and resultant pathology with this different isolate, we analyzed this data set separately.

In the 2017 experiment, wild-caught House Finches were recruited into an experimental challenge study between 6 July and 17 September 2017. Following an initial quarantine and housing in groups of three, 24 naïve birds were used in a transmission experiment in which they were housed in pairs. As mentioned earlier, one bird in each pair (the index bird) was treated either with meloxicam or water daily, beginning 2 d prior to inoculation and ending on day 21. Over 44 d, we evaluated pathology (15 times), pathogen load (nine times for nonindex birds and seven times for index birds, respectively), and antibody titers (three times). Data from both index and nonindex birds were included in community analysis.

Across both experiments, a single observer (J.S.A.) recorded the presence (1) or absence (0) of a community of pathological descriptors: blepharospasm, crusting, erythema, eversion, exudate, and swelling. Eyes that received the same severity score (scale: 0.5–3) pre and postpeak pathology were then included for community analysis, such that a single bird could have been represented multiple times across eyes and scores. When birds received the same eye score consecutively, we used the sample closest to peak pathology. We had rare instances in which eyes scored the maximum eye score (3) over multiple days. To accommodate these samples, we used the first and last day an eye received this severity score.

The analytical techniques that we used were borrowed from community ecology, treating descriptors as species and eyes as unique sites that supported different communities of pathological descriptors. We used R package vegan (Oksanen et al. 2019) to calculate Jaccard's dissimilarity indices, followed by multidimensional scaling based on principal coordinates analysis with R package geomorph to visualize community differences. We then used analysis of variance to evaluate main effects of eye score, status, and their interaction on the distance matrix using a randomized residual permutation procedure (Collyer and Adams 2018). Last, we used R package indicspecies (De Cáceres and Legendre 2009) to carry out multilevel pattern analysis, which allowed us to identify indicator species (descriptors) predictive of status (pre- versus postpeak pathology).

We digitized 85 paired images of 52 unique eyes from 27 birds to evaluate the effects of lesion severity and stage of infection on eye shape. This data set included the full range of observer-assigned eye scores: 46 pairs of score 0, 19 pairs of score 0.5, 13 pairs of score 1, four pairs of score 1.5, one pair of score 2, and two pairs of score 3. We found significant main effects for eye score (0–3; F_5,164=8.92, Z=6.30, P=0.001), status (pre- or postpeak; F_1,168=5.33, Z=2.84, P=0.001), and their interaction (F_5,164=3.17, Z=3.52, P=0.001) on eye shape. However, the score by status interaction should be interpreted with caution because we had small sample sizes for eye scores >1.5. Still, pairwise comparisons across all score and status combinations suggest that prepeak scores of 1.5 have a unique shape signature, as they significantly diverge from pre- and postpeak scores of 0 and 0.5 and postpeak scores of 1 (Z>2.00, P<0.035), as well as postpeak scores of 1.5 (Z=5.33, P=0.001). All other pairwise comparisons with prepeak scores of 1.5 were marginally significant (all P=0.05–0.10). Just one other contrast (pre-0:pre-1) was significant (Z=2.27, P=0.019).

Based on PCA, we found that PC 1 and PC 2 explained 63.8% of the variation in eye shape. We then plotted specimens along these axes to visualize occupied disease space (Fig. 3). We found a clear trend in eye shape such that distortion of the eye rims increased with severity score, resulting in specimens with higher severity scores being located further from the origin. The PCA plot further revealed that eyes returned to the same relative location in disease space following the resolution of pathology (black squares, Fig. 3).

To fully interpret the path through disease space, we introduced status (pre- versus postpeak pathology) to provide some measure of temporal resolution across specimens. Because birds reached different peak pathologies on different days, using status allowed us to bin specimens by score and stage of infection to calculate mean centroid locations for birds as they moved through disease space (Fig. 4a). Again, the locations for scores >1.5 should be interpreted with caution due to small sample sizes. In general, however, infected eyes appeared to move clockwise through disease space (Fig. 4a), returning to baseline in the same approximate location they started. This pattern was consistent when we plotted individual trajectories for eyes with enough paired samples to fully capture the disease cycle (n=6; Fig. 4b).

The low-virulence data set included 39 paired samples of 32 unique eyes from 19 birds. Of these, 19 pairs had observer-assigned eye scores of 0.5, 13 had scores of 1, four had scores of 1.5, one had a score of 2, and two had scores of 3. The high-virulence data set included 29 paired samples of 14 unique eyes from nine birds. Of these, six pairs had scores of 0.5, eight had scores of 1, eight had scores of 1.5, three had scores of 2, two had scores of 2.5, and two had scores of 3.

Separately for each data set, we extracted the ordinal scores for site (eye-day) from the distance matrix and plotted them along the principal coordinate axes to evaluate trends in community composition (Fig. 5). Mean ordinal location shifted down the PCoA1 axis with lesion severity, quantified by observer-assigned eye score, in accordance with richer descriptor communities. Score (F_4,73=5.21, Z=5.57, P=0.001) and the interaction between score and status (F_4,73=1.67, Z=1.81, P=0.042) were significant sources of variation in the low-virulence distance matrix, whereas only main effects of score (F_5,52=2.73, Z=4.32, P=0.001) and status (F_1,56=3.73, Z=3.10, P=0.001) were significant sources of variation in the high-virulence distance matrix.

Multilevel pattern analysis on the low-virulence data set identified exudate as a marginally significant indicator of eyes in prepeak pathology (P=0.07). Exudate was present in 36% (14/39) of prepeak eyes compared to 15% of postpeak eyes (6/39), meaning that when exudate was present, eyes were two times more likely to be early in the disease course. Of the 20 specimens with exudate, 70% (14/20) were in prepeak pathology. In contrast, crusting was a significant indicator of eyes in prepeak pathology in the high-virulence data set (P=0.029). Crusting was present in 38% (11/29) of prepeak eyes compared to 10% of postpeak eyes (3/29), meaning that when crusting was present, eyes were three times more likely to be early in the disease course. Of the 14 specimens with crusting, 79% (11/14) were in prepeak pathology. No other descriptors separated out by status.

We introduced two novel techniques for capturing patterns in multivariate trait space using one of the best-characterized wildlife disease systems with external pathology, finch mycoplasmosis. Borrowing methods from geometric morphometrics and community ecology, we demonstrated how multivariate analyses can improve our understanding of pathology in this system.

Based on shape analysis, we found that finch eyes occupied unique areas of multidimensional disease space depending on the severity and stage of infection. Specifically, eyes with mild disease (for example, bird eye 480.L, peak eye score=1.5/3; Fig. 4b) tracked much closer to the origin than did eyes with severe disease (for example, bird eye 318.R, peak eye score=3/3; Fig. 4b). All eyes also moved through disease space in a clockwise progression. Given this consistency, it is plausible that we could characterize the stage of infections using single snapshot evaluations of pathology in wild individuals. In particular, distinguishing between recently infected animals and those in recovery could improve our temporal resolution when tracking new disease introductions or seasonal epizootics (Pepin et al. 2017).

Based on community analysis, we found that the composition of pathological descriptors changed with severity score. Specifically, as observer-assigned eye scores increased, descriptor communities occupied different areas of multidimensional disease space defined by unique suites of pathological descriptors (Fig. 5). In addition, we found that that certain phenotypes were linked to eyes in early-stage disease. Specifically, exudate was a weak indicator during infection with the low-virulence isolate, while crusting was a strong indicator during infection with the high-virulence isolate. If we consider exudate and crusting along a continuum, it follows that pathogenicity can intensify this phenotype when abundant exudate leads to crusting. Thus, evaluation for the presence of certain pathologies, as with tracking changes in eye shape, can improve the temporal resolution of disease investigations. In particular, this could lend itself to direct field applications, such as deciding when to collect swabs for pathogen isolation, because samples from recent infections tend to grow better in culture (Raka 2012). Therefore, evaluating for the presence of exudate or crusting could make for an easy screening procedure for birds with pathology due to wild-type infection. It could also identify the individuals that could benefit most from treatment or rehabilitation because they are still early in the disease course. This may be especially useful when managing infectious diseases of threatened or endangered species or other high-value populations with limited resources.

Our results build upon recent efforts to infer stage of infection and predict disease outcomes based on multivariate data collected at discrete time points. Torres et al. (2016) found that mice infected with the malarial parasite Plasmodium chabaudi progressed through disease space in fairly predictable loops. Although their analyses reflected only two univariate metrics at once, the authors found that early infection performance could predict infection outcomes—when individuals entered danger zones in disease space from which they would not recover. We saw a hint of this in bird eyes 75.R and 318.R, which originated further to the right along PC 1 and went on to experience worse pathology than their cohort (Fig. 4b). Perhaps baseline eye shape explains some of the individual variation we see in disease progression, predisposing birds to certain types of conjunctival pathology. Evaluating such differences in multivariate disease trajectories may reveal similar danger zones to those seen in malarial infections and, from a clinical perspective, improve early prognostic indicators for individuals requiring therapeutic intervention. Moreover, by leveraging baseline eye shapes to predict certain disease phenotypes, we may also be able to identify individuals that pose a higher health risk to their population. For example, eversion was a frequent wild-type pathology encountered during field capture and could enhance pathogen shedding if everted tissues are more likely to come in contact with fomites. Examining these linkages across the multivariate data sets described here may further prove illuminating.

Our results suggest that measurement of external pathology with multivariate techniques captures information that is distinct from observer-assigned severity scores, which can improve our understanding of disease outcomes and host-pathogen interactions more broadly. Such methods could open up a new branch of inquiry in disease ecology and offer unique insights into myriad diseases affecting free-ranging wildlife.

We thank Amali Stephens and Sarah Tosh for their assistance with inoculations for the low-virulence study that generated experimental infection data analyzed here. The National Science Foundation supported this work (IOS-1950307 to J.S.A.).