> \*Correspondence: *Nickolay I. Hristov<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Conservation, a section of the journal Frontiers in Ecology and Evolution*

> Received: *01 January 2019* Accepted: *25 June 2019* Published: *23 July 2019*

#### Citation:

*Hristov NI, Nikolaidis D, Hubel TY and Allen LC (2019) Estimating Overwintering Monarch Butterfly Populations Using Terrestrial LiDAR Scanning. Front. Ecol. Evol. 7:266. doi: 10.3389/fevo.2019.00266* (Hristov et al., 2008; Chilson et al., 2012). Human history is rich with examples of the impact and contribution of such tools as vehicles for the exploration of new frontiers. Population size is an important variable for understanding the biology, ecology, and conservation of any organism. However, estimating the size of large group-living organisms, like the monarch butterfly, under field conditions has been notoriously challenging (Hristov et al., 2010). LiDAR technology, a 3D measurement technique, and its associated computational tools offer a promising method to help tackle the challenge and assess the abundance and density measures needed to build on decades-long efforts to monitor monarch populations in Mexico and the US.

Monarch butterflies exhibit one of the most impressive displays of animal migration, spanning multiple generations and migrating over thousands of kilometers across the continent. Monarchs east of the Rocky Mountains overwinter in dense aggregations in high-elevation Oyamel fir (Abies religiosa) forests in Mexico, while western monarchs overwinter more diffusely at numerous sites in trees on the California coast (Urquhart and Urquhart, 1976; Malcolm, 2018). Regular monitoring of the overwintering population in Mexico began in 1994 by World Wildlife Fund-Mexico in collaboration with a number of other organizations (including the Comisión Nacional De Áreas Naturales Protegidas and Secretaría De Medio Ambiente Y Recursos Naturales; Thogmartin et al., 2017a). The Mexican overwintering population is mapped as total forest area occupied in hectares which is used as an index of population size (Rendón-Salinas et al., 2017). Annual western overwintering population estimates in California began in 1997 (Malcolm, 2018) and records are kept by the Xerces Society1 , but estimates, via markrecapture methods, were made as early as the 1970s (Tuskes and Brower, 1978). Currently, the California overwintering estimates are produced via direct visual counting of clusters using an "area method" in groves along the coast2 .

While the methods of estimation differ, data show a consistent decrease in overwintering numbers for both eastern and western monarch populations (Brower et al., 2012; Semmens et al., 2016; Rendón-Salinas et al., 2017; Schultz et al., 2017; Malcolm, 2018). By some measures, monarch abundance may have declined more than 80% over the past 25 years, leading to a petition in 2014 for listing under the Endangered Species Act3 . US Fish and Wildlife Service (USFWS) is leading a monarch status assessment, to be finalized by 2020, to determine whether the monarch should be listed. The accurate assessment of the size and seasonal dynamics of overwintering monarch butterflies, therefore, is increasingly urgent.

One challenge for completing an accurate status assessment is difficulty determining accurate numbers of the eastern monarch population, due to the size of the task. While monarchs settle in up to 300 sites in coastal California, estimated millions concentrate in just a few locations in Central Mexico (Urquhart and Urquhart, 1976; Malcolm, 2018). As mentioned above, occupied area estimates have been used as a proxy for monitoring the eastern monarch population (Rendón-Salinas and Tavera-Alonso, 2014; Vidal and Rendón-Salinas, 2014; Vidal et al., 2014). There is considerable variation in density estimates at these sites, from 6.9 to 60.9 million butterflies per hectare, making it difficult to obtain accurate numbers–numbers which are critical for conservation efforts (reviewed in Thogmartin et al., 2017a). Regular monitoring has produced evidence of population decline; however, with monarchs now clustered on fewer and fewer trees, more precise and accurate estimates of their population at the individual level are critical. One of the most promising methods of estimating density is the capturerecapture method replicated over multiple years, (Calvert, 2004; Thogmartin et al., 2017a), however the method is disruptive to overwintering butterflies and is currently not allowed. To overcome this limitation, we look to technology to help build upon the decades of work to more accurately and precisely estimate the density and numbers of overwintering monarchs in a non-invasive way (Hristov et al., 2008). The solution, while more critical for use in the eastern overwintering sites, due to their concentration over a smaller area, will be useful for the entire North American population of the species.

To address this urgent need, we propose a novel approach for monarch population estimation using terrestrial LiDAR scanning (TLS) and associated digital tools. TLS creates a detailed digital survey using millions of LASER-based distance measurements from a scanning device to objects in the landscape. We have created innovative workflows and analytical algorithms for obtaining and converting TLS scans and high-resolution imagery of overwintering monarch aggregations into cumulative volumes and estimates of monarch abundance. Terrestrial LiDAR scanning and other photographic techniques have recently been used by biologists to visually count individual bats (Azmy et al., 2012; Shazali et al., 2017), and butterflies (Leong et al., 2017) in clusters. The technology has also been used in other large-scale applications like estimating tree biomass (Momo Takoudjou et al., 2018). In contrast, we expand the use of TLS to generate high-resolution 3D models and use the volumetric data in the models to estimate the number of butterflies in these overwintering sites.

In this technology report we discuss the use of TLS and imaging methods and the development and application of an innovative Subtractive Volume Estimation technique for estimating overwintering monarch butterfly populations. The Subtractive Volume Estimation (SVE) method works by comparing the detailed geometry of a site with and without the animals of interest, deriving as a result an estimate of the overall bio-volume and eventually an estimate for the number of individuals and density. The SVE method has shown great promise to enable the estimation of large clusters of animals and can help increase population count precision and accuracy through technology.

#### MATERIALS AND METHODS

#### Overwintering Sites

We collected data from four monarch overwintering sites; two in California, USA (Santa Cruz Lighthouse Field State

1www.Xerces.org/monarchs/

2www.pgmuseum.org/blog/2014/10/17/how-do-you-count-monarch-butterflies 3www.biologicaldiversity.org/species/invertebrates/pdfs/Monarch\_ESA\_Petition. pdf

FIGURE 1 | The analysis of Terrestrial LiDAR Scanning (TLS) data follows a numerical model based on shape analysis. Subtractive Volume Estimation (SVE) compares a scan with butterflies, referred to as the *reference scan,* to a scan without butterflies, known as the *derivative scan* and calculates the difference between the two shapes. That difference represents the volume of butterflies present. The final estimated colony size is based on dividing that volume by a scaling factor, *n*, which represents the volume of an individual butterfly multiplied by a correction factor f that accounts for characteristics of the TLS survey.

Park and Pismo State Beach Butterfly Grove) representing the western population of the species and two sites in central Mexico (Carpinteros and San Antonio) representing the eastern population4 . The data was collected over 3 years from 2017–2019.

# Terrestrial LiDAR Scanning

To collect TLS data, we used a pair of FARO S70 long-range laser scanners (Faro Technologies, Lake Mary, Florida, USA). This model has a maximum scanning resolution of 700 megapixels and is capable of scanning objects up to 70 m away with a ranging accuracy of ± 1 mm at that distance. The S70 also has the ability to generate a 350 megapixel, composite, color image based on 64 pictures captured with its on-board camera. The composite image can then be overlaid on top of the point cloud during postprocessing of the data to give a realistic representation of the color space of the surveyed environment. Individual scan settings ranged from ½ to ¼ of the maximum resolution, based on the complexity of the surveyed environment, with lower settings used for more cluttered environments and higher-resolution scans used in more open, uncluttered environments. Scanning quality, the number of successful readings before a distance-estimate is recorded, was set at 2x; therefore, we collected approximately 175 million and 44 million distance measurements per scan at ½ and ¼, respectively. The beam diameter at exit for the S70 was 2.12 mm (beam divergence 0.3 mrad). The vertical and horizontal step size of the scanner was 0.009◦ . On average, we positioned the scanners 5.39 meters away from butterfly clusters of interest resulting in a minimum resolvable feature diameter of 1.53 and 3.07 mm at ½ and ¼ resolutions respectively; well below the dimensions of the roosting butterflies. Each scan took between 3 and 12 min so an entire survey generally took < 2 h. We used a pair of scanners to speed up the surveys and worked early in the morning, between sunrise and 10:00 local time. Any measurements of overwintering groves must be taken before ambient temperatures reach ≈16◦C, which is the flight threshold for monarchs (Masters et al., 1988). Above that temperature, the butterflies grow increasingly active and begin to leave their clusters. TLS surveys were composed of 8– 25 scan positions depending on the size and complexity of the

4www.xerces.org/wp-content/uploads/2015/10/MonarchMap-NatureServe-10. 20.png

surveyed site. Our goal was to capture most of the groves where butterflies were roosting; from 2–10 trees on approx. 0.2–0.4 hectares at California sites and 0.6–1.7 hectares with dozens of trees in central Mexico. The resulting point-cloud datasets were represented by over 352 million points in California and 4.4 billion points in Mexico. We were able to capture every portion of the habitat where butterflies were roosting at the California sites and large sections of the habitat at the sites in Mexico. From these composite datasets we selected clusters with different densities of butterflies for further analysis.

#### Subtractive Volume Estimation

The principle proposition of our Subtractive Volume Estimation (SVE) method is to compare the volumetric differences between two TLS surveys, a reference scan that records the roosting monarch butterflies in their overwintering environment and a derivative scan, that records the same site without the butterflies. The subtraction of the two geometric volumes results in an expression of the biovolume of the clustering insects. To convert from total volume of butterflies to number of individuals, we introduce a scaling factor, **n**, which represents a species-specific volumetric estimate of the space occupied by one individual. Monarch butterflies arrange themselves differently depending on the ambient temperature, characteristics of their environment and overall size and density of the cluster (Masters et al., 1988; Vidal and Rendón-Salinas, 2014). In addition, differences in the precision of the laser survey due to distance from the scanner, angle of incidence of the laser beam and surface reflectivity require an additional variable, correction factor **f**, that accounts for these characteristics of the TLS instrument and survey process. The result is an estimate for the number of butterflies in a subset cluster of the overall group or the entire colony (**Figure 1**).

#### Data Processing

The individual, raw LiDAR scans were processed and registered together into a composite point cloud using FARO Scene5 , a proprietary software provided by the scanner manufacturer. Additionally, we applied outlier, edge and stray-point filters in FARO Scene that removed any noise in the data. The point cloud data was then exported for further manipulation by other computational tools.

Few analytical tools exist for point cloud data, therefore, to make use of the richness of LiDAR data, the point clouds must first be converted to polygonal geometry—a rather computationally intensive process (**Figure 2**). Additionally, the polygonal 3D models of these clusters must be enveloped in complete surfaces to be suitable for 3D metrology. The subtraction of volume calculation can be done in either of two ways: (1) direct shape subtraction where the surveyed objects or scenes do not change between the reference and derivative scan (e.g., tree trunk; see **Figure 11A** and **Supplementary Video 2** for example) and the subtracted metrics have a spatial meaning and (2) an arithmetic subtraction of the volumes when organic structures, such as branches, change their shape or position between the reference and derivative scans due to the weight of the butterflies (**Figure 11B** and **Supplementary Video 3**).

To analyze the clustering shapes, we categorize the roosting aggregations of butterflies into three types based on features of the clusters—(A) Type 1, tree-trunk clusters, (B) Type 2, distal branch clusters, and (C) Type 3, canopy clusters (**Figure 3** and **Supplementary Video 1**). Overwintering sites along the California coast feature mostly Type 2 clusters while sites in central Mexico display all three types.

To estimate the individual number of butterflies in the total subtracted volume of a cluster, we calculated both the scaling factor, **n**, and correction factor, **f**, as described above. To derive **n** we developed a physically accurate 3D model of an individual monarch butterfly using scaled reference images of naturally dead monarch butterflies found at the survey

5www.faro.com/products/construction-bim-cim/faro-scene

sites. There are no substantial size differences between male and female butterflies and both sexes overwinter at the same sites. The model allowed us to calculate the maximum speciesspecific volume of an individual butterfly. The scaling factor, **n,** for Danaus plexippus, is 49 cm3 (**Figure 4**). Simulations described below showed that butterfly volume **n** derived from virtual pointclouds were not constant and therefore needed additional modification. We derive the correction factor **f** based on two approaches that both require starting with known numbers of butterflies: (1) The 3D model was used as an instancing reference—a 3D modeling technique that allows the efficient generation of large, repetitive geometries based on a single reference (see below; **Figure 4**), and (2) we identified instances in the LiDAR surveys where we could count individual butterflies or small clusters of them roosting on tree surfaces and extracted their volumes. While (1) initially indicated that **n** is not constant in clusters with different densities, necessitating the formulation of **f**, (2) provided a more realistic measure based on the performance of the LiDAR scanners in natural clusters.

# Testing and Modeling

Validation is an important component of computational and simulation methods as such approaches generally produce a numerical result; but, it is important to understand if such estimates are representative of biological reality. Since invasive or destructive methods of testing and validation that require the handling of butterflies were not an option for this project, to assess the accuracy of the estimation process, we developed a virtual, physically accurate, 3D model of roosting butterflies and arranged them on (1) tree-trunks and (2) distal branches. At present, because of the size, threedimensional complexity, inability of LiDAR to see through opaque surfaces and need for additional proofing data, Type 3, canopy clusters remain a challenge and are the focus of ongoing efforts.

We describe the validation process for the tree-trunk models, but a similar procedure is used for distal branch models. We created digital 3D models of bare tree trunks and distal branches of Oyamel and Eucalyptus trees (Eucalyptus globulus), based on laser scans in California and Mexico (**Figure 5**) and combined these models with the 3D model of a monarch butterfly (described above) to create virtually simulated reference and derivative point clouds. In this virtual environment we were able to "place" known numbers of butterflies in the scene to validate the performance of the SVE algorithm. We used an adapted, virtual photogrammetric method to create point-cloud representations of these virtual 3D structures. Photogrammetry, also known as structure from motion, stitches highly overlapped digital images into an accurate 3-dimensional representation of the scene (Edward et al., 2001). The photogrammetry program, Agisoft Photoscan6 , analyzes the position and rotation difference of specific patterns between different photos and derives the shape of objects on those photos, as well as the position of the camera that took the photo in relation to those objects. Through experimentation we found that 1 degree of separation between renders was adequate to accurately and effectively calculate a point cloud. We rotated the camera so that it is always facing the center of the trunk, and we took a render at each 1◦ interval. We processed the entire image-sequence in Agisoft Photoscan to calculate a collective point cloud representation of the photographed geometry. A similar process captured the tree trunk without butterflies. The two point clouds were imported into a point-cloud editing program called MeshLab7 where we used a Screened Poisson Surface Reconstruction algorithm (Kazhdan and Hoppe, 2013) to generate a polygonal surface from adjacent points in the point cloud. The result was two cylindrical 3D models representing the bare tree, and the tree populated with butterflies. Next, the 3D models were closed to make them "watertight," having no holes in the geometry, a requirement for volumetric calculations, and the resulting objects were imported into 3D Coat8 , a digital sculpting program that treats 3D objects as voxel-based volumes. In this program there are tools that allow the accurate subtraction of one volume from another, and the export of resulting volume as a new object. The subtracted volume difference between the bare trunk and trunk with butterfly models was imported once again into MeshLab, where we calculated the volume of the end shape (**Figure 6**). The process was repeated for a variety of butterfly clusters (with varying known numbers of virtual butterflies)

6www.agisoft.com

7www.meshlab.net

8 3dcoat.com/home

reference scan, including the branch with butterflies and derivative scan of the branch without butterflies. These quantitatively accurate models are used to validate the results from the TLS estimation of real-world butterfly aggregations.

to derive a direct connection between changes in volume and butterfly numbers.

The workflow sequence presented here uses generally available software tools which make them more accessible to potential users and facilitate adoption of this new methodology in practice. Other computational and metrology environments exist such as the industry-standard, Innovmetric Polyworks9 which provide integrated and more consistent and robust editing and analytical capabilities for higher-order and larger-scale computations when larger clusters or entire sites must be estimated.

#### RESULTS

In this technology report we describe the development of an analytical workflow that begins with the acquisition of threedimensional TLS data in the field and ends with an estimate of the number of butterflies in a section of the roosting habitat.

#### Derivation of Correction Factors

When comparing the simulated 3D models of clusters with an increasing number of butterflies with their surveys derived from virtual photogrammetry, we observed, as expected, that the total volume of the subtractive volume envelopes increased from 503 cm3 for 10 butterflies to 177,780 cm3 for 2,000 butterflies (**Figures 6**, **7**). However, when dividing the total subtractive

volume for each cluster by the number of butterflies, known exactly from the simulation engine, to arrive at a representative volume for an individual butterfly, the estimated volume per butterfly was not constant as might be expected (reference values from **Figure 8**). Instead, it changed as a function of the number of butterflies in the simulated clusters in a non-linear way.

the subtracted volume shapes for butterfly clusters of 100–2,000 individuals,

indicating increasing density of the resulting shapes.

If the virtual photogrammetric surveys recorded correctly the geometry of the simulated roosting butterflies, we would have had a 1:1 relationship between the number of butterflies entered in the simulation and the number of butterflies estimated through the SVE process. However, this is not the case. This inconsistency is explained by the relationship between the number of butterflies and the distance among them and how TLS and photogrammetry represent this inter-individual spacing. Noting that when butterflies are fewer (10–100), thus

9www.innovmetric.com/en

further apart, the laser and photogrammetric processes correctly detect individual butterflies as separate geometric objects and accurately represent their volumes. On the other end of the spectrum, in dense clusters of many butterflies (≥ 2,000), individuals are arranged closely together and their collective mass is captured accurately by TLS/Photogrammetry as a single envelope. However, during intermediate densities (300–1,000), the spacing among butterflies falls into a range where the virtual photogrammetric process cannot distinguish individual butterflies as separate objects and fuses them into one, filling the intermediate gaps between neighboring butterflies–a product that we call phantom volume. As a result of this additional volume, the estimate for the calculated volume of an individual butterfly (i.e., division of total envelope by correct number of butterflies) incorrectly increases (**Figures 9**, **10**) at these intermediate values (300–1,000 butterflies). Because of these observations we computed the correction factor **f** that is included as part of the scaling factor **n** used for estimating the number of individual butterflies in the total subtracted volume from simulated and natural clusters.

## Estimation of Natural Clusters

Based on these calculations and using values of 49 cm3 and 0.29 for **n** and **f** respectively, values derived from a small patch of butterflies on a twin-trunk tree at the San Antonio site, we estimate that 192 butterflies occupy a volume of 0.002657 m3 in that cluster. Thus, 2,641 butterflies occupy 0.0365 m3 of volume for every 1-meter length of trunk with Type 1 clusters exclusively. Given the approximately 12 meter section covered by butterflies

on each trunk of the tree (24 meters combined for both trunks), we estimate that 63,378 butterflies, occupying 0.8770 m3 roosted on the entire tree (**Figure 11A** and **Supplementary Video 2**). Similarly, for an SVE analysis of a tree from the Carpinteros site, consisting of Type 1 and Type 2 clusters, we estimate that the total volume of the clusters is 0.3917 cubic meters, resulting in an estimated 36,995 butterflies for the entire tree (**Figure 11B** and **Supplementary Video 3**).

spaced butterflies are not detected correctly and phantom volume is generated, necessitating the introduction of a correction factor (see text and

#### DISCUSSION

Figure 11).

The overwintering population size of monarch butterflies in Mexico, expressed as hectares of forest, has become the most commonly cited measure of yearly monarch abundance. Millions to hundreds of millions of monarch adults may inhabit these few hectares of central Mexican forest in the winter. There is a great variation in published estimates of overwintering density (Brower et al., 2004; Calvert, 2004) making extrapolation

to the entire population difficult. Thogmartin et al. (2017a) used 6 published density estimates to develop a probability distribution for monarch overwintering density in Mexico. While they were able to calculate a median density of 21.1 million per hectare, they acknowledge that considerable uncertainty in the estimated densities remain. Such uncertainty in overwintering densities (and population counts) is a problem when building conservation interventions and particularly when trying to assess the effectiveness of such efforts.

volume of a single individual and *f* represents the derived correction factor for

each monarch density.

In this technical report, we share progress and early achievements in the development of a novel workflow and analytical methods for extracting monarch abundance and density information from TLS. We confirm the suitability and exciting potential of TLS in combination with digital analytical techniques such as SVE, to provide a reliable and non-disruptive method for more accurately estimating overwintering monarch butterfly densities and population size. We demonstrate that, 3D, polygonal reconstructions based on TLS surveys have ample resolution and accuracy to enable the separation of butterflies from the background vegetation and the estimation of butterfly clusters from a few individuals in a simple arrangement to thousands in intricate clusters. Furthermore, we have developed a unique, non-destructive method for validating our estimates, via quantitative 3D modeling and virtual environments.

With these achievements, we deliver on the primary goal of this initiative–to provide an alternative approach to monarch estimation at the individual level that can expand on the decades of monitoring work in North America. Such efforts are well suited for establishing a repeatable standard for calculating overwintering density per unit areas. The work presented here can significantly reduce the uncertainty in density estimates that are still very much part of the conversation about monarch conservation (Thogmartin et al., 2017a). Therefore,

this novel technological approach could be combined with the current area estimate methodology (Rendón-Salinas et al., 2017) for practical monitoring efforts that use the two techniques hand-in-hand.

Despite these promising results, important work remains as we continue to expand the application of the approach to larger sections of forest and clusters with more complex structures. Improved speed, accuracy, portability and new developments in analytical capabilities will help the process and fine-tune the workflows. Still we recognize that, despite its promise and important capabilities, novel application and great promise, TLS/SVE have several shortcomings: At present, we have not attempted to estimate the number of butterflies in Type 3, canopy clusters as we do not know enough about their internal structure. TLS is a surface survey tool, thus dense, canopy clusters in Mexico will require further data collection and analysis to understand the intricate structure of the interior of these massive aggregations. For example, it may be equally possible that these clusters are densely packed with butterflies or largely hollow and thus any assumption, one way or the other, would significantly affect the outcome and estimates in each case (**Figure 3A**). Additional observations, scanning techniques and detailed analyses are underway to inform the application of this method with regards to these Type 3 clusters. Similarly, TLS scanners are expensive and current computational methods for analysis of the data are

FIGURE 12 | A reflective IR photograph of a cluster of monarch butterflies at Pismo Beach State Park Monarch Butterfly Grove, California, showing the absorption of monarchs in the IR part of the electromagnetic spectrum. This effect has potential to be used to identify visually clusters of monarchs in complex environments and thus aid computer vision analysis.

tedious, time-consuming and requiring trained staff to perform both the data collection and analytics.

Nevertheless, TLS combined with SVE brings important and promising capabilities to establish the foundation for future development of individual-based whole-colony assessment. With its high inherent accuracy and smooth workflow, TLS/SVE sets the groundwork for moving this agenda forward. Perhaps the greatest promise of such future developments is the use of spectral imaging in combination with photogrammetry. Just as we described in the validation work discussed here, photogrammetric (aka Structure from Motion) reconstructions from serial sequences of images can produce detailed pointcloud data that can be used for 3D models to analyze for census estimates. This capability, combined with emerging evidence that spectral imaging, particularly IR photography, provides strong separation of monarch butterflies from background vegetation a challenging computer vision problem (**Figure 12**). The two technologies combined, could offer an even more efficient and affordable technique for monarch butterfly enumeration. Further technical developments could provide important new proxies that could develop in a similar way to bat census work over the last two decades (Betke et al., 2008; Kloepper et al., 2016; Hristov pers. obs.). In that work initial efforts were intensive and expensive (multiple years and ∼\$500K) while current work is quicker and cheaper (a few evenings and \$1–2K) which offers agencies, field workers, policy managers, and others accessible, effective tools for managing and conserving resources. As this technology matures, and more user-friendly and capable processing environments for stitching reconstructions emerge, the early work described by TLS will be superseded by these more accessible technologies and tools.

In the meantime, TLS and associated analytics offer an important formalization of the methods and practices to collect, organize and archive critical temporal and spatial information on monarch butterfly distribution and population trends. TLS not only provides important survey data that can serve as the starting point for further analysis and estimation, it serves as a detailed record in space and time to archive the overwintering behavior of the species, similar to efforts around the world to protect and archive sites of historic and cultural importance by groups like CYARK10 This is particularly important in the face of unprecedented change in distribution and population numbers for the monarch butterfly.

#### CONCLUSSION

The Species Status Assessment for the monarch is currently underway by FWS with completion required by 2020. Improved abundance estimates are important for use in monarch extinction risk models that have been created for the Assessment. An ultimate goal is to vastly improve our ability to quantify the number and density of butterflies at overwintering sites, particularly in the highly localized eastern monarch sites in central Mexico. This paper describing the establishment of a robust TLS workflow and the development of SVE-based analytical algorithms improves our chances of making that goal a reality and these efforts stand to benefit the field of animal enumeration tremendously. Ultimately, this method is one of several early steps to understanding the dynamic aggregation of monarchs at a finer scale that will allow for more precise and accurate density and populations estimations.

# AUTHOR CONTRIBUTIONS

NH conceived and designed the study. LA and NH collected the data. NH, LA, DN, and TH analyzed the data and contributed to technology adaptations. NH and LA wrote the manuscript.

#### ACKNOWLEDGMENTS

We are extremely grateful to our colleagues who helped facilitate access and work at field sites, including Ronnie Glick, Danielle Bronson, with California State Parks, and Samantha Markum with USFWS and our collaborators with CONANP, WWF-Mexico, and the stewards of the Monarch Biosphere Preserve, Mexico. We are thankful to undergraduate assistants at Winston-Salem State University, Georgina Dzikunu, Zaria Jean, and Timothy Nixon, for help with TLS processing and analysis and to Raunak Kapoor and Trent Spivey for help with data collection. We also thank Ryan Lebar for insight into spectral photography. We thank Ryan Drum and Shauna Marquardt from USFWS and Ralph Grundel and Wayne Thogmartin from

10https://www.cyark.org

#### REFERENCES

Azmy, S. N., Sah, S. A. M., Shafie, N. J., Ariffin, A., Majid, Z., Ismail, M. N. A., et al. (2012). Counting in the dark: non-intrusive laser scanning for population counting and identifying roosting bats. Sci. Rep. 2:524. doi: 10.1038/srep00524

USGS for constructive conversations that helped shape this work. We wish to thank Lisa Stack, without whose logistical expertise this work would not have been accomplished and Chris White, Tommy Maddox, and Chris Morabito from Faro Technologies for technical support with LiDAR. Two reviewers (KC and JR) improved greatly the quality and clarity of the manuscript for which we are especially grateful. This work was supported by the US Department of the Interior via DOI US Fish and Wildlife Service awards (F18AC00004 and<PHONE_NUMBER>) and the Commission for Environmental Cooperation, with additional support provided by National Science Foundation DRL-1514776. Publication costs were supported by a Winston-Salem State University Faculty Professional Development Committee research grant to LA. Funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00266/full#supplementary-material

Supplementary Video 1 | A high-resolution, three-dimensional, animated model of an Oyamel Fir tree in central Mexico covered in monarch butterflies. The model is generated based on TB LiDAR which is then converted into high-resolution, polygonal mesh. The intricate detail of the model allows the analysis of elaborate clusters and roosting patterns of the butterflies. We identify three clustering patterns on these trees: *Type 1*, tree trunk clusters; *Type 2*, distal branch clusters and *Type 3*, canopy clusters (see also Figure 3).

Supplementary Video 2 | A high-resolution, three-dimensional, animated model of a pair of Oyamel Fir trees in central Mexico (San Antonio site) covered in monarch butterflies with Type 1 clusters. The model is generated based on TB LiDAR which is then converted into high-resolution, polygonal mesh. The intricate detail of the model allows the analysis of elaborate clusters and roosting patterns of the butterflies. Segment at 00:00–00:07 shows the *point cloud without butterflies*; segment at 00:07–00:13 shows the high-resolution *polygonal mesh without butterflies*; segment at 00:13–00:19 shows the *point cloud with butterflies*; segment at 00:19–00:24 shows high-resolution *polygonal mesh with butterflies*; segment at 00:24–00:30 shows the *combined*, with and without butterflies, mesh models; segment at 00:30–00:36 shows the *subtracted volume* of the butterflies used for the number of individuals estimate (see also Figure 11A).

Supplementary Video 3 | A high-resolution, three-dimensional, animated model of an Oyamel Fir tree in central Mexico (Carpinteros site) covered in monarch butterflies with Type 1 and Type 2 clusters. The model is generated based on TB LiDAR which is then converted into high-resolution, polygonal mesh. The intricate detail of the model allows the analysis of elaborate clusters and roosting patterns of the butterflies. Segment at 00:00–00:06 shows the *point cloud with butterflies*; segment at 00:06–00:12 shows the high-resolution *polygonal mesh with butterflies*; segment at 00:12–00:18 shows the *point cloud without butterflies*; segment at 00:18–00:24 shows the high-resolution *polygonal mesh without butterflies*; segment at 00:24–00:30 shows the *combined*, with and without butterflies, mesh models; segment at 00:30–00:36 shows the *subtracted volume* of the butterflies used for the number of individuals estimate (see also Figure 11B).

Betke, M., Hirsh, D., Makris, N. C., McCracken, G. F., Procopio, M., Hristov, N. I., et al. (2008). Thermal imaging reveals significantly smaller brazilian free-tailed bat colonies than previously estimated. J. Mammal. 89, 18–24. doi: 10.1644/07-MAMM-A-011.1

Available online at: http://assets.panda.org/downloads/monitoreo\_mariposa\_ monarca\_en\_mexico\_2013\_2014.pdf (accessed December, 2018).

**Conflict of Interest Statement:** DN is the owner of Animusing Productions.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Hristov, Nikolaidis, Hubel and Allen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Size of the Canadian Breeding Population of Monarch Butterflies Is Driven by Factors Acting During Spring Migration and Recolonization

Tara L. Crewe1 \* † , Greg W. Mitchell 2 and Maxim Larrivée3

1 National Data Center, Bird Studies Canada, Port Rowan, ON, Canada, 2 Wildlife Research Division, Environment and Climate Change Canada, Ottawa, ON, Canada, 3 Insectarium de Montréal, Montréal, QC, Canada

#### Edited by:

Jay E. Diffendorfer, United States Geological Survey (USGS), United States

#### Reviewed by:

Nathan Lemoine, Colorado State University, United States Chip Taylor, University of Kansas, United States

#### \*Correspondence:

Tara L. Crewe<EMAIL_ADDRESS>
#### †Present address:

Tara L. Crewe, College of Engineering, IT and Environment, Charles Darwin University, Darwin, NT, Australia

#### Specialty section:

This article was submitted to Population and Evolutionary Dynamics, a section of the journal Frontiers in Ecology and Evolution

> Received: 15 January 2019 Accepted: 31 July 2019 Published: 14 August 2019

#### Citation:

Crewe TL, Mitchell GW and Larrivée M (2019) Size of the Canadian Breeding Population of Monarch Butterflies Is Driven by Factors Acting During Spring Migration and Recolonization. Front. Ecol. Evol. 7:308. doi: 10.3389/fevo.2019.00308 The eastern North American monarch butterfly population shows a long-term population decline. While it is hypothesized that forest loss on the wintering grounds and milkweed loss throughout the breeding range are responsible for the observed decline, there is much less certainty regarding the factors driving year-to-year variation around the current population level. Using 15 years of butterfly count data, we used a community-based approach to delineate the stage of the annual cycle during which population limiting factors are most strongly acting. We compared annual fluctuations in size of the breeding population of monarch butterflies in Canada to fluctuations in 13 additional butterfly species which either migrate long distances to Canada or are resident but breed in similar habitats to the monarch. We show that the breeding population of monarchs in southern Canada shows a higher degree of synchrony with other long-distance migrants than with breeding residents, and that annual fluctuations of all migrant butterflies show a positive correlation with the number of 21◦C days during spring migration and re-colonization. Further, we found that size of the monarch breeding population shows a higher degree of synchrony with the size of the following winter population than with the size of the previous winter population. Combined, our results suggest that the monarch population in Canada is limited by factors acting during spring migration, and that weather plays an important role in the ability of the monarch to successfully re-colonize and breed in the northern portion of their summer range each year. A predicted increase in temperature in the early spring, combined with continued loss of breeding and wintering habitat, has the potential to limit the reproductive capacity of monarchs and their ability to recover from population lows.

Keywords: monarch butterfly, Danaus plexippus, population limitation, community ecology, temporal synchrony, weather

# INTRODUCTION

Monarch butterfly (Danaus plexippus) populations are declining throughout their eastern range, during all phases of the annual cycle (fall migration: Crewe and McCracken, 2015, but see Badgett and Davis, 2015; over-wintering: Thogmartin et al., 2017a; breeding: Pleasants and Oberhauser, 2013), and its migration is considered by some to be an endangered biological phenomenon (Brower et al., 2012). While it is generally hypothesized that milkweed loss throughout the monarch's breeding range and forest loss on the monarch's wintering grounds are largely responsible for the 84% population decline observed since the mid-1990s (Stenoien et al., 2016; Malcolm, 2017; Thogmartin et al., 2017b), there is much less certainty regarding the primary factors driving year-to-year variation around the overall negative population trend. For example, Thogmartin et al. (2017b) found that early warm temperatures in the northern breeding range negatively affected population size on the wintering grounds the following winter, but that warmer temperatures later in the season had a positive effect, likely due to effects on the growing conditions for milkweed. Conversely, Badgett and Davis (2015), Ries et al. (2015), and Inamine et al. (2016) suggested a lack of correlation between fall population indices in the US and population size on the wintering grounds is evidence that limiting factors are occurring during fall migration. Given these competing, yet not mutually exclusive hypotheses, there is a need to better understand when, and where limiting factors are occurring in the annual cycle of the monarch butterfly.

In Canada, the monarch butterfly is listed as "Special Concern" but is being considered for uplisting to the status of "Endangered" under the Species at Risk Act (SARA), and occurs in greatest numbers in southern and eastern Ontario, where they are at the northern limit of their breeding range (COSEWIC, 2016; Environment Climate Change Canada, 2016). The population size of monarchs in Canada may be limited by factors acting throughout the monarch's annual cycle. For example, in the north central US, summer breeding populations are best predicted by higher amounts of precipitation, and cooler temperatures in Texas during early migration/breeding, likely reflecting ideal growing conditions for milkweed in the south (Zipkin et al., 2012; Saunders et al., 2016, 2018). Indeed, correspondence between changes in egg density and milkweed abundance in the upper mid-western United States with the number of individuals migrating in fall in Canada (Crewe and McCracken, 2015) suggests that breeding habitat and/or weather conditions experienced at breeding locations south of Canada are important. Northern breeding populations could also be limited by habitat availability and/or weather conditions as well as local density dependence on the breeding grounds (e.g., Marini and Zalucki, 2017; Thogmartin et al., 2017b). Alternatively, previous winter's population size in Mexico has been found to be positively correlated with the following winter's population size (Thogmartin et al., 2017b), suggesting a positive ripple effect through successive breeding generations that may affect breeding numbers reaching the north. Last, conditions during fall migration such as nectar availability and suitable weather conditions for migration may influence migration success and the ability of monarchs to reach the wintering grounds (Badgett and Davis, 2015; Ries et al., 2015, and Inamine et al., 2016), potentially affecting Canada's overall contribution to the winter population.

Until now, research addressing monarch population limitation has generally been conducted in isolation from the broader butterfly community. However, if weather and habitat change are acting not just on monarchs, but on the entire butterfly community, then the amount of temporal synchrony among species with similar life history strategies has the potential to reveal broad-scale relationships between observed patterns of change, and the processes driving those changes (Michel et al., 2016). Monarchs are among several long-distance migratory butterfly species that breed in southern Canada; several migratory species overwinter in the southern United States, north of the monarch's Mexican over-wintering grounds. Because these additional migrants do not overwinter in the same region as monarchs, temporal synchrony among those migrants on their Canadian summer breeding grounds would lend support to the hypothesis that migrant butterflies, including monarchs, are most limited by factors influencing the success of the spring migration and recolonization period. In contrast, many species that breed in similar habitats to the monarch are non-migratory residents in southern Canada. Greater temporal synchrony between monarchs and resident species would support the hypothesis that breeding habitat quantity or quality and/or weather conditions directly affecting survival and reproduction during breeding limit southern Canada's breeding monarch population.

We use monarch overwintering density estimates (Semmens et al., 2016; Thogmartin et al., 2017a) and butterfly survey data collected in the form of checklists by the Ontario Butterfly Atlas (2003–2017; Macnaughton et al., 2017), and eButterfly (2012–2017; Prudic et al., 2017) to test four competing hypotheses on whether monarchs breeding in southern Canada are most strongly limited by factors acting during fall migration, winter, spring migration and recolonization, or summer breeding (**Table 1**). We predicted that if monarchs are limited by factors acting during fall migration ("fall limitation hypothesis"), there will be weak correlation between relative abundance on the breeding grounds and subsequent wintering density; and if monarchs are limited by factors acting on their overwintering grounds ("winter limitation hypothesis"), there should be a strong correlation between previous winter density and relative abundance on the breeding grounds during the following summer. Alternatively, we predicted that if monarchs are limited by conditions during spring migration and re-colonization ("spring limitation hypothesis"), there will be (a) weak correlation between previous winter population density and breeding abundance, and (b) stronger temporal synchrony between monarchs and other migratory butterflies, than between monarchs and resident butterfly species. If monarchs are limited by conditions on their southern Canada breeding grounds ("summer limitation hypothesis"), we predicted that there should be (a) strong correlation between summer relative abundance and monarch densities during the following winter, and (b) stronger synchrony between monarchs and resident butterfly species than between monarchs and species which undergo long-distance migrations. By incorporating community-level comparisons of a 15 year time series dataset and comparing indices of summer and winter population densities, we differentiate between these competing limitation hypotheses, and test when during the annual cycle monarchs breeding in Ontario are most likely limited.

TABLE 1 | Competing hypotheses testing when during the annual cycle monarch butterfly populations in Canada are most limited.

Predictions for correlations between wintering and breeding populations are not mutually exclusive; by also testing predictions for temporal synchrony among migrant and resident butterfly species, we can better distinguish support for or against the competing hypotheses.

# METHODS

#### Data Collection Overwintering Population Size

Estimates of the overwintering density of monarchs (2002–2014) were obtained from Thogmartin et al. (2017a), which provides the raw observed hectares collected by the World Wildlife Fund Mexico, fitted hectares as estimated by Semmens et al. (2016), and associated predicted population size in millions of individuals (see Thogmartin et al., 2017a for details).

#### Ontario Breeding Population Size

Annual indices of population size for the monarch butterfly and additional migratory and non-migratory butterfly species that share similar breeding habitats (**Table 2**) were estimated using butterfly checklist data from the Ontario Butterfly Atlas (hereafter "Atlas") for the years 2003–2017. The Atlas includes occurrence and abundance data for Ontario from several sources including museum collections, eButterfly, Butterflies and Moths of North America (BAMONA), and iNaturalist (Macnaughton et al., 2017). eButterfly is a crowd-sourced checklist-based webplatform for gathering presence/presumed absence data from across North America (Prudic et al., 2017). Regional experts verify the validity of submitted observations on both Atlas and eButterfly platforms. Records submitted to the Atlas were excluded when: (1) they had missing day, month, or year, (2) "x" was listed as the count, or (3) "road-kill" or "specimen" were listed as the record type.

Because Atlas data are submitted as independent observations for each species, a "checklist" is defined here as observations submitted by an observer or group of observers on a particular date for a particular location (latitude/longitude, rounded to the nearest 2 decimal places, or ∼ 1-km accuracy); i.e., each unique combination of date, latitude, longitude, and observer(s) is considered a checklist. Prior to calculating total number of species observed (hereafter "list length") on a checklist, we TABLE 2 | Species included in community level analyses of temporal synchrony.

Species were classified as either long-distance migrants ("Migratory") or resident; resident species were further classified as either univoltine, two-generation, or multi-generation.

aligned the butterfly taxonomy where it differed among databases represented in the Atlas (**Supplementary Table 1**). Records of unknown or very rare species were excluded from the data, and records for individuals with unknown species but known genus (e.g., Colias sp.) were included in calculations of list length only when the Genus was not already present on a list. In total, 160 species were detected by the Ontario Atlas and included in the calculation of list length (**Supplementary Table 1**).

After calculating list length, we used the open GIS software QGIS (v.2.18.16) to select lists that fell within the predicted breeding range of the monarch, based on known milkweed and monarch distributions (Larrivée, unpubl. data; **Supplementary Figure 1**). Lists with fewer than 4 species were also removed to avoid including rare species reports (Breed et al., 2013). Because list length declined with an increase in latitude, we controlled for the effect of latitude on list length by using the residuals from a linear regression of the log of list length by latitude as a surrogate for list length (**Supplementary Figure 2**). We further filtered data for each species by excluding lists collected before or after the minimum and maximum observation date for that species across years.

#### Weather Data

We used the NCEP.gather function of the RNCEP package (Kemp et al., 2012) in R (R Core Team, 2017) to extract surface-level weather variables for the lower-, mid-, and upper United States in April, May and June, respectively, from the National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) Reanalysis data set (Kalnay et al., 1996). These data have a spatial resolution of 2.5◦ × 2.5◦ and a temporal resolution of 6 h (00, 06, 12, 18 h UTC; Kemp et al., 2012). For the month of April, weather data were extracted for the lower U.S., which we defined as the area encompassed by latitudes in the 30 to 35 degree range, and within −104 to −83 degrees longitude. In May, we extracted weather data from the mid-U.S., which we defined as the area

encompassed by latitudes in the 35 to 40 degree range, and within −104 to −83 degrees longitude. In June, we extracted data for the upper-U.S., defined as the area encompassed by latitudes in the 40 to 45 degree range, and within −104 to −83 degrees longitude. For each month, we calculated the total number of ≥ 21◦C days, and for April and May, total number of freezing (≤ 0 ◦C) days. The mean number of 21◦C days and freezing days across locations in a region was used as an index of annual "21◦ days" and "0 days," respectively. We chose the total number of ≥ 21◦C days because Thogmartin et al. (2017b) found that this number for early May was negatively correlated with proceeding winter population sizes. Data for May 0◦C days were not complete in 2009, and 2009 was excluded from those analyses.

#### Breeding Population Estimates

We estimated annual indices of population size for each butterfly species (**Table 2**) by fitting a modification of the list length analysis for checklist data described in Breed et al. (2013) and Szabo et al. (2011) to the Ontario Butterfly Atlas data. Because monarchs are large and conspicuous, probability of detection nears 0.75 or greater at the end of the breeding season during both low and high abundance years (**Supplementary Figure 3**). We therefore used abundance and assumed a negative binomial distribution for overdispersed counts in place of an occurrence model with a binomial data distribution as in Breed et al. (2013) and Szabo et al. (2011). For each species we used hierarchical linear regression to model the number of individuals observed on a list as a function of fixed effects for list length (residuals of the linear regression of log list length by latitude, as described above) and year. The parameter for list length was assumed to control for all factors that influence detection, including effort, observer skill, weather, time of day, etc. (Szabo et al., 2011; Breed et al., 2013). We re-parameterized the model by removing the intercept for year; in doing so, an effect size for each year, after accounting for random variability among years, was output for the model. These year effects provided a year-specific estimate of abundance, which we used as an index of relative abundance in the correlation and synchrony analyses described below. Random first-order autoregressive random effects for year and day of year nested within year were also included to account for temporal autocorrelation of counts among years and days of the year, respectively. All models were fit in a Bayesian framework using integrated nested Laplace approximation with the INLA package (Rue et al., 2009; Martins et al., 2013) in the R statistical programming language (R version 3.4.0; R Core Team, 2017).

# Temporal Synchrony

To tease apart which of the four competing hypotheses (**Table 1**) was best supported by the data, we tested (1) the correspondence of monarch overwintering and breeding annual indices, (2) community-level correspondence of indices among migrating and resident butterfly species, and between monarchs and other migrant or resident species that breed in similar habitats; and (3) correspondence of butterfly and weather indices for the time period during the annual life cycle deemed to be driving population fluctuations based on results from (1) and (2).

#### Correspondence of Breeding Monarch Butterflies and Overwintering Densities

We tested whether size of the breeding population was correlated with the size of the previous and following overwintering population (raw observed hectares and fitted hectares, as described above) using Spearman rank correlation for both the raw and de-trended breeding and wintering annual population indices. De-trending removes any underlying linear trend from a time series by taking the residuals of a linear model fit to the indices and allows fluctuations in indices to be compared after accounting for any linear trend over time. Indices were detrended using the detrend function of the RSEIS R package (Lees and Harris, 2008). Because temporal autocorrelation of annual indices can bias cross-correlation estimates, we adjusted the p- and t- values of the r estimate by correcting for the loss of degrees of freedom due to temporal autocorrelation, following the methods described in Kirchner (2001), Michel et al. (2016). Next, we used the peaks function of the synchrony R package (Gouhier and Guichard, 2014; R Core Team, 2019) to determine the proportion of maxima and minima that corresponded between the breeding and wintering population time series, for both raw and de-trended annual indices of abundance.

#### Synchrony Among Resident and Migratory Butterfly Species

For community level correspondence of annual indices, we first tested the correlation of de-trended annual indices among species-species pairs using Spearman rank correlation, correcting for any temporal autocorrelation among annual indices (Kirchner, 2001; Michel et al., 2016). We then used a hierarchical linear model with the correlation coefficient (r) as the dependent variable, comparison type (migrant to migrant [M-M]; migrant-resident [M-R]; and resident to resident [R-R]) as an explanatory factor, and species 1 and 2 as random effects to account for repeated measures across species in the species-species comparisons. For the 5 migrant and 9 resident species (**Table 2**), this resulted in 10 M-M, 36 R-R, and 45 M-R correlations. Next, we used the community.sync function of the synchrony R package (Gouhier and Guichard, 2014) to estimate the community-level synchrony across (1) all species, (2) migrant butterfly species, and (3) resident butterfly species. If monarchs are limited by factors acting during spring migration and recolonization, migrant species should show a higher degree of synchrony amongst each other, than with southern Canada resident butterfly species. Alternatively, if monarchs are limited by factors acting during the breeding season, there should be greater synchrony between monarchs and resident butterfly species. Finally, we compared the correspondence of maxima and minima in the time series of species-species pairs using the peaks function in the synchrony R package. Using the same modeling structure described above, we also fit a hierarchical linear model with proportion of corresponding peaks as the dependent variable. Last, we repeated the same analysis with the proportion of corresponding peaks as the dependent variable, but only considering comparisons between the monarch and each other species.

#### Correspondence of Weather Variables With Breeding Population Indices

We tested the correlation of detrended weather and breeding annual indices for migrant butterflies using Spearman rank correlation, with t- and p-values adjusted for temporal autocorrelation of counts as described above and in Kirchner (2001), Michel et al. (2016). We also used the peaks function in the synchrony R package to estimate the proportion of corresponding maxima and minima between the de-trended annual indices of each migrant and each weather variable (mean number of 21◦C days and 0◦C days in April in the lower U.S., May in the mid-U.S., and June in the upper mid-west).

# RESULTS

Between 2003 and 2017, total number of lists submitted to the Ontario Butterfly Atlas increased significantly (linear regression: DF = 1, 13, F = 13.57, p = 0.003), but mean list length did not (linear regression: DF = 1, 17,352, F = 0.28, p = 0.59). Across years, mean and median list length were 8 and 7

the Ontario Butterfly Atlas (2003–2017). Indices were estimated using hierarchical linear regression that assumed a negative binomial data distribution, with count as the dependent variable, and with residuals of list length by latitude (a measure of effort), and year (factor) as explanatory variables. Temporal autocorrelation among

years and among days within years were accounted for through the specification of random year and day effects.

TABLE 3 | Spearman correlation (N.Eff, effective sample size after accounting for temporal autocorrelation of indices; r, correlation coefficient; t, t-value; p, p-value) and proportion of corresponding maxima and minima (Prop. peaks) of breeding annual indices estimated using Ontario Butterfly Atlas data with overwintering population estimates (density, fitted Hectares, and observed hectares) for the winter previous and following the breeding season, and using raw and de-trended annual indices.

Values in bold indicate correlations that were statistically significant at p ≤ 0.05.

species, respectively. Because list length declined with increasing latitude, we used residuals from the linear regression of log list length by latitude (DF = 1, 17,352, F = 145.3, p < 0.0001; **Supplementary Figure 2**) in our trend analyses, to model the effect of list length after accounting for variation due to latitude. Estimated annual indices for the 5 migrant and 9 resident butterfly species are shown in **Figure 1**.

#### Temporal Synchrony Correspondence of Wintering and Breeding Population Estimates

Annual indices for the southern Canada monarch breeding population had a significant positive correlation only with the following winter's fitted hectares and density estimates (**Table 3**), which lends support to the summer limitation hypothesis and suggests a lack of support for the fall limitation hypothesis (**Table 1**). Weak correlation with the previous winter's estimates (**Table 3**) lends support to the spring limitation hypothesis and suggests a lack of support for the winter limitation hypothesis (**Table 1**). Proportion of corresponding maxima and minima ("Prop. Peaks," **Table 3**) in the time series was also greater between the summer breeding and following winter populations, than between the breeding and previous winter populations (**Figure 2**), again supporting the summer and spring limitation hypotheses, respectively. De-trending the time series did not increase the proportion of corresponding peaks, but did result in an increase in the Spearman correlation coefficient from 0.28 to 0.66 (p = 0.02) and 0.70 (p = 0.01) for the fitted hectares and density overwintering estimates, respectively (**Table 3**; **Figure 2**), suggesting that breeding and wintering populations show a greater correspondence in the direction of annual fluctuations than among linear trends in counts over time.

#### Correspondence Between Migrant and Resident Butterfly Species

Mean correlation of annual indices for migrant-migrant species comparisons was greater than mean correlation among migrant-resident and resident-resident comparisons, and 95% confidence intervals of migrant-migrant and migrant-resident groups did not overlap. This result suggests that the annual indices of long-distance migrants are more strongly correlated with the annual indices of other migrants, than they are with annual indices of resident butterfly species (**Table 4**, **Figure 3A**; for raw correlation estimates see **Supplementary Table 2**) and suggests stronger support for the spring limitation hypothesis than for the summer limitation hypothesis. The stronger correlation between migrant-migrant pairs compared with migrant-resident pairs was maintained when residents were classified into groups according to number of generations in a year (univoltine, 2-generation, or multi-generation; **Supplementary Figure 4**). Spearman correlation coefficients were also significantly greater when monarch annual indices were compared to those of other migrants, than when compared with Ontario breeding residents (linear regression: intercept = 0.38, p < 0.001; group (M-R) = −0.27, p = 0.01; **Figure 3B**).

For the butterfly species compared, community-wide synchrony was estimated to be 0.44 (p = 0.01) across all 14 migrant and resident butterfly species, 0.48 (p = 0.01) across the 9 resident species, and 0.62 (0.04) across the 5 migrant species. When residents were broken down into groups based on number of breeding generations, community-wide synchrony increased to 0.64 (p = 0.04) among multi-generational residents, 0.89 (p = 0.01) among 2-generational residents, and 0.49 (p = 0.04) among univoltine residents.

The spring limitation hypothesis was also supported by results of proportion of corresponding peaks, which was greatest among migrant-migrant comparisons. Confidence limits (95%) of the M-M estimate did not overlap with those of the migrant-resident and resident-resident estimates, which supports that migrants are showing more similar fluctuations in annual indices as a group than when migrants are compared with resident butterfly species (**Table 4**, **Figure 4A**; for raw proportion of corresponding peaks for each species-species pair, see **Supplementary Table 3**). As with correlation estimates, the stronger relationship between migrant-migrant pairs than migrant-resident pairs was maintained when residents were further broken down into groups according to number of annual generations (**Supplementary Figure 5**). When comparisons were restricted to those that included the monarch butterfly, proportion of corresponding maxima and minima between the monarch's annual indices and indices of other migrants was significantly greater than when monarch annual indices were compared to those of resident butterfly species (linear regression: intercept = 0.68, p < 0.0001; group (M-R) = −0.42, p < 0.0001; **Figure 4B**).

#### Correspondence of Weather Variables With Breeding Population Indices

Annual indices of migrant butterfly species tended to show greater correspondence with the number of 21◦C days in April

not shown).

TABLE 4 | Estimates for mean Spearman correlation (r) and proportion of corresponding peaks for species-species comparisons between migrant species (M-M, n = 10), between migrant and resident species (M-R; n = 45), and between resident species pairs (R-R; n = 36).

Estimates were derived from a hierarchical linear model where r or proportion of corresponding peaks was the response variables, group was a predictor variable, and species 1 and species 2 were included as random effects to account for repeated measures on each species across species-species pairs.

and May than with 21◦C days in June and the number of 0◦C days in April or May (**Figure 5A**), though the only statistically significant relationships were for the correlation between annual indices of the monarch (r = 0.68, p = 0.01) and question mark butterflies (0.56, p = 0.05) and the number of 21◦C days in April in the lower U.S. (**Figure 6**; raw correlation coefficients in **Supplementary Table 4**). Correspondence of maxima and minima in the time series was greatest for all migrant species with the number of 21◦C days in May in the mid-U.S. and June in the upper mid-west (**Figures 5B**, **7**; raw proportion of peaks in **Supplementary Table 4**). Proportion of corresponding peaks was significant between annual indices and May 21◦C days for monarch (0.70, p = 0.02), painted lady (0.73, p = 0.03) and American painted lady (0.64, 0.04) butterflies, and between annual indices and June 21◦C days for monarch (0.70, p = 0.04) and American painted lady butterflies (0.73, p = 0.02).

#### DISCUSSION

Our unique butterfly community analysis allowed us to assess and disentangle the relative importance of migration and summer breeding factors as drivers of annual and interannual variability in monarch population size in southern Canada. We show that weak correspondence between the size of the overwintering population and subsequent breeding population, along with stronger temporal population synchrony among migrants than among resident butterflies, supports the hypothesis that monarchs breeding in Canada are most

FIGURE 3 | Model derived estimates (black circles) and associated standard errors of the spearman rank correlation of species-species annual population indices for models that tested for a difference among groups (M, migrant; R, resident), for (A) a model that looked at all species-species comparisons across the 14 migrant and resident species (M-M: n = 10; M-R, n = 45; R-R, n = 36), and (B) a model that included only species-species comparisons that included the monarch butterfly (M-M, n = 4; M-R, n = 7). In both cases, the raw spearman correlation coefficient for each species-species comparison is also plotted for each group (open circles).

models that tested for a difference among groups, for (A) a model that looked at all species-species comparisons across the 14 migrant (M) and resident (R) species (M-M: n = 10; M-R, n = 45; R-R, n = 36), and (B) a model that included only species-species comparisons that included the monarch butterfly (M-M, n = 4; M-R, n = 7). In both cases, the proportion of corresponding peaks for each species-species comparison is also plotted for each group (open circles).

strongly limited by factors acting during spring migration and recolonization given contemporary levels of milkweed availability in the US (see Thogmartin et al., 2017c). The lack of synchrony between annual numbers of monarchs and resident co-occurring butterfly species suggests that the size of the monarch breeding population is not driven by factors acting during the breeding season, such as direct and indirect effects of weather (e.g., effects on nectar availability or milkweed quality). Instead, interannual synchrony of all migratory butterfly species breeding in southern Canada over the 15 year time series examined in this study suggests that factors acting during spring migration and recolonization are driving observed patterns in the year-to-year variability in monarch annual population indices.

Weak correspondence between the size of the overwintering population and subsequent southern Canada breeding population did not support the hypothesis that the size of the Ontario breeding population is predominantly limited by factors acting on the overwintering grounds. This hypothesis is also not supported by our result that monarch numbers fluctuate in sync with all long distance migrant butterflies on their southern Canada breeding grounds, but not with 9 common co-occurring resident species of butterflies that breed in the same habitats as monarchs in southern Canada. Therefore, the relationship between overwintering and subsequent breeding population sizes likely breaks down as a result of variation in reproductive potential with weather conditions during spring migration and recolonization (Zipkin et al., 2012; Saunders et al., 2016, 2018).

The correspondence of monarch (and other migrant) annual indices with the mean number of 21◦ days during spring migration is in agreement with Thogmartin et al. (2017b) and suggests that given favorable breeding conditions, monarchs have the potential to recover from low winter population sizes.

Indeed, while the pattern of ups and downs were strongly correlated between the breeding and following winter population estimates, correspondence was poor when the time series were not de-trended; i.e., the strong decline in overwintering population size, observed prior to more recent increases since population lows in 2012–2014, was not reflected in the breeding population. However, Crewe and McCracken (2015) found that the magnitude of change between peaks and lows in the number of monarchs migrating out of Ontario each fall has declined over time, which might be indicative of a reduction in reproductive and re-colonization potential. If so, the monarch population may become more susceptible to further declines over time, particularly if spring weather conditions become increasingly warm with climate change (e.g., Schwartz and Reiter, 2000; Zipkin et al., 2012; Saunders et al., 2016, 2018).

Strong correspondence between breeding and following winter population sizes did not support the hypothesis that the current breeding population of monarchs in Canada is predominantly limited by mortality during fall migration (Badgett and Davis, 2015; Ries et al., 2015; Inamine et al., 2016). Also, observed shifts in the phenology of monarch butterflies, such that monarchs are staying in southern Canada later in the year (Prudic et al., 2017, e.g., Zipf et al., 2017), could have an impact on fall migration mortality, if later migration results in a higher probability of encountering extreme weather events on the breeding grounds or during migration. Research is required to determine whether a shifting phenology on the breeding grounds is resulting in additional breeding generations each year or delayed migration, and whether individuals delaying their migration, or a new late summer or early fall generation, are likely to survive a later fall migration.

Migratory monarch butterfly populations are currently extremely low and vulnerable (Semmens et al., 2016). Loss of milkweed breeding habitat as a result of glyphosate application to genetically modified crops is thought to be largely responsible for current population levels (Pleasants and Oberhauser, 2013; Flockhart et al., 2015; Stenoien et al., 2016; Pleasants, 2017). Historical forest loss on the wintering grounds in Mexico due to logging has also been hypothesized to affect population size and has been responsible for the disappearance of a wintering colony (Brower et al., 2002; Vidal et al., 2014). Further, nectar resource availability has likely declined through time (Potts et al., 2010), potentially limiting the availability of nectar resources for monarchs to build fat reserves late in the fall migration in the southern US and Mexico, which are required to sustain them throughout the winter (Brower et al., 2006). While our study suggests that fluctuations in the number of breeding monarchs counted in southern Canada since 2003 are driven by factors acting during spring migration and recolonization, our results do not imply that factors acting during alternative phases of the annual cycle, including the availability of milkweed during the summer breeding season, nectar resources during fall migration, and forest habitat during the overwintering period, have not contributed to historical rates of population decline, or are not limiting the potential of monarchs to recover to previous population levels. Further research is needed on historical and contemporary milkweed availability in Canada to better understand how land management practices affect current breeding populations of monarchs and other pollinators. Research should also address whether increasing rates and severity of extreme weather events in the fall since the turn of the 21st century, such as hurricanes and extended drought periods, are influencing monarch survival during fall migration. Finally, while overwintering forest loss has sharply declined since 2008 (Vidal et al., 2014), increasing overwintering habitat could help increase population numbers observed in Canada, particularly in years where growing conditions for milkweed in the southern US are favorable. It should also be noted that severe weather events on the wintering grounds, which were not explicitly accounted for in our analysis, have the potential to limit next summer breeding success in northeastern North-America in any given year. For example, large die-offs of monarchs on the wintering grounds owing to severe storm events (e.g., Brower et al., 2017) might affect the number of fall migrants returning to Mexico the following winter. We suggest that further research is needed on the interactions between events on the wintering grounds and the subsequent spring migration and recolonization to fully understand contemporary population limiting processes of the migratory monarch population.

In conclusion, we show that our list-length approach using data that was largely collected by citizen scientists to develop population indices for 14 butterfly species breeding in Canada, speaks to the value of community science data and butterfly

FIGURE 6 | Detrended annual indices for breeding monarch butterflies and detrended mean number of 21◦ days in April in the lower United States, by year and, inset, showing linear relationship between detrended monarch and detrended mean April 21◦ day indices. A significant Spearman correlation of 0.71 (p = 0.01) was detected.

mid-eastern United States. Significant correlations between detrended breeding annual indices and detrended mean 21◦ days were not detected, though proportion of corresponding maxima and minima ranged between 0.60 and 0.82.

checklists and suggests our approach could easily be extended to other parts of the monarch range. Our community level analysis which allowed us to compare inter-annual population indices of several co-occurring migratory and resident butterfly species in southern Canada also revealed for the first time that contemporary populations of monarchs in Canada are most strongly limited by weather events occurring in the US during spring, a result supported by multiple analyses in the US (Zipkin et al., 2012; Saunders et al., 2016, 2018). Moving forward, it will be important to understand how conservation actions taken by all three countries throughout the annual cycle of the monarch interact with stochastic weather effects, and where those interactions are strongest, so that conservation planners can better prioritize where remedial actions may be most effective.

### AUTHOR CONTRIBUTIONS

ML, GM, and TC conceived and designed the experiments. ML, TC, and GM contributed data and analysis tools. TC performed the experiments and analyzed the data. TC wrote the manuscript, with input from GM and ML.

#### REFERENCES

# FUNDING

Funding for this work was provided by Environment and Climate Change Canada.

#### ACKNOWLEDGMENTS

We are grateful to the many individuals whose butterfly observations are incorporated into the Ontario Butterfly Atlas each year.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00308/full#supplementary-material

of monarch butterflies overwintering in Central Mexico. PeerJ. 5:e3221. doi: 10.7717/peerj.3221

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Crewe, Mitchell and Larrivée. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Monarch Butterfly Conservation Through the Social Lens: Eliciting Public Preferences for Management Strategies Across Transboundary Nations

Rodrigo Solis-Sosa1 \*, Christina A. D. Semeniuk 2 , Sergio Fernandez-Lozada1 , Kornelia Dabrowska1 , Sean Cox 1 and Wolfgang Haider 1†

*1 School of Resource and Environmental Management, Simon Fraser University, Burnaby, BC, Canada, 2 Great Lakes Institute for Environmental Research, University of Windsor, Windsor, ON, Canada*

#### Edited by:

*Jay E. Diffendorfer, United States Geological Survey (USGS), United States*

#### Reviewed by:

*Paul Cross, Bangor University, United Kingdom Jose Roberto Soto, University of Arizona, United States*

> \*Correspondence: *Rodrigo Solis-Sosa<EMAIL_ADDRESS>
> > *† In memoriam*

#### Specialty section:

*This article was submitted to Conservation, a section of the journal Frontiers in Ecology and Evolution*

Received: *23 December 2018* Accepted: *06 August 2019* Published: *22 August 2019*

#### Citation:

*Solis-Sosa R, Semeniuk CAD, Fernandez-Lozada S, Dabrowska K, Cox S and Haider W (2019) Monarch Butterfly Conservation Through the Social Lens: Eliciting Public Preferences for Management Strategies Across Transboundary Nations. Front. Ecol. Evol. 7:316. doi: 10.3389/fevo.2019.00316* The monarch butterfly *(Danaus plexippus)*, an iconic species that migrates annually across North America, has steeply declined in numbers over the past decade. Across the species' range, public, private, and non-profit organizations aim to reverse the monarch decline by engaging in conservation activities such as habitat restoration, larvae monitoring, and butterfly tagging. Urban residents can actively participate in these activities, yet their contribution can also be realized as an electorate body able to influence the design of conservation programs according to their interests. Little is known, however about their preferences toward the objectives and design of international monarch conservation policies. In this paper, we investigate these preferences via a survey design using Discrete Choice Experiments (DCEs) and Latent Class Analysis (LC) of urban residents across the main eastern migratory flyway in Ontario, Canada, and the eastern United States. Attributes in the DCE included the size and trend of overwintering butterfly colonies, the type of institution leading the conservation program, international allocation of funds, and the percentage of funds dedicated to research. From the general populace, we isolated respondents already engaged in monarch conservation activities to explore how they compare. We sent a smaller set of surveys deliberately withholding the expected-success forecast of the monarch recovery program to assess the value of information for urban residents within a conservation context. The LC distinguished three groups of respondents among urban residents: (1) the main group, labeled "Eager," accounting for 72.4% of the sample, that showed a high potential for supporting conservation policies and had remarkable similarities with the monarch enthusiasts' sample; (2) a "Pro Nation" group (18.4%) marked by their increased willingness to support conservation initiatives solely focused within their country of residence; and (3) an "Opinionated" segment (9.23%), that was highly reactive to changes of the leading institution, resources allocation, and economic contribution proposed. Key findings from this research reveal that to maximize potential support amongst urban residents in the monarch's breeding range, a conservation strategy for the monarch butterfly should be led by not-for-profit organizations, should strive for transboundary cooperation, and should include the communication of anticipated ecological outcomes.

Keywords: monarch butterfly, citizen science, choice experiment, latent class, conservation, public preferences, international cooperation, transboundary conservation

#### INTRODUCTION

The design of conservation strategies for transboundary migratory species has proven to be a challenging topic for decision makers, partly due to the presence of multiple institutions, groups of interest, administrative barriers, and political and cultural differences (Grant and Quinn, 2007). The monarch butterfly (Danaus plexippus; henceforth referred to as "monarch") is a highly migratory and globally distributed butterfly species (Oberhauser et al., 2008). Its eastern North American population has the longest migration (Brower and Missrie, 1999)—up to 4,000 km—in which butterflies across the eastern states/provinces of the US and Canada establish overwintering colonies within a few specific forest patches in Mexico (Urquhart and Urquhart, 1976). The area occupied by monarchs in these overwintering sites has decreased from an average area of 5.71 ha in 1993 to an all-time-low area of 0.67 ha during the 2013–2014 season (Vidal and Rendón-Salinas, 2014). Its current estimate is at 6.05 ha (Rendón-Salinas et al., 2019).

Habitat destruction in both overwintering and breeding areas is currently the most plausible hypothesis for the population decline (Brower et al., 2012). Overwintering monarch colonies rely on the forest canopy for protection against freezing temperatures, precipitation, and wind (Anderson and Brower, 1996). During the breeding season, monarchs depend on milkweed (Asclepias spp) for larvae feeding across the breeding range from Northern Mexico to the northeastern US and eastern Canadian provinces (Zalucki et al., 2001). Here, agricultural land transformation combined with the introduction of transgenic-specific herbicides for crop management, to which only genetically modified crops can resist, have caused a general decline in milkweed abundance across the eastern states of the US over the last decade (Pleasants and Oberhauser, 2012).

In addition to those multiple stressors, the heterogeneous sociopolitical backdrop of the monarch's decline is a considerable challenge as well. Different resource-extraction activities, socioeconomic differences (Lopez-Hoffmann et al., 2009), and distinct legislative tools and processes for its protection (Waples et al., 2013) can hinder the effectiveness and coherence of joint strategies (Scott and Collins, 1997). Moreover, the limited resources available for conservation, from governments and NGOs alike, are allocated based on national priorities, which may significantly differ between countries. For example, while the monarch butterfly is a top priority for WWF-Mexico with more than 25 scientific monarch-related reports (WWF- Mexico, 2018), its Canadian office 2017 annual report has no mention of the monarch butterfly (Miller, 2017). Additionally, since political institutions tend to be responsible for internalizing environmental externalities, with their mandates focusing on local issues, externalities at an international level are frequently overlooked (Perrings and Halkos, 2012). One example of such an externality is the potential loss of revenue that Mexican communities incur from monarch-reserve tourism caused by extensive use of herbicides in the North (Esquivel-Rios et al., 2014).

Despite this intricate sociopolitical backdrop, the monarch's decline and its widespread appeal have spurred people's interest in its protection across the migratory flyway. For example, The Monarch Waystation program, an initiative seeking to stimulate the public to provide habitat for monarchs and other pollinators, is continually increasing its presence every year, with 21,946 registered waystations up to date (Lovett, 2018). Likewise, Journey North, an entry-level citizen science platform, received 1,574 reports of egg sightings and 14,381 adult sightings during fall 2017, contrasted with 193 eggs and 3,310 adults reported in 2012. Several other citizen science hubs have witnessed that same surge of interest by the general public such as the Monarch Watch Tagging Program, eButterfly, iNaturalist, and the Monarch Larvae Monitoring Program. Such participation of private residents in citizen science and ecologically-related activities provide scientists with an extraordinary capacity of having useful, cost-effective data collected and analyzed. Ries and Oberhauser (2015) estimated that 17% of 503 papers related to the monarch since 1940 have relied on citizen science data. Moreover, Lewandowski and Oberhauser (2017) found that individuals engaged in citizen science activities are more likely to provide and protect critical habitat as well.

However, the role of the general public in protecting the monarch, as well as any other imperiled species, can go beyond data gathering and habitat provision—at which farmers could be substantially more effective (Thogmartin et al., 2017). Instead, when a conservation target is embedded within a complex network of economic and cultural interests at a transboundary level as described above, the involvement of the general public is especially needed. Agnone (2007) studied how the general public's opinion and protests have impacted the passage of environmental laws in the United States between 1960 and 1998. Several national conservation policies have been successful when the public is engaged. For example, Lutrin and Settle (1975) documented the success of passing California's Coastal Zone Conservation Act due to the active engagement of the Coastal

**Abbreviations:** ASC, Alternative Specific Constant; DCE, Discrete Choice Experiment; LC, Latent Class or Latent Class Analysis; MNL, Multinomial Logit Model; NEP, New Environmental Paradigm Scale; RI, Relative Importance; WTP, Willingness to Pay; mWTP, Marginal Willingness to Pay.

Alliance with the public and contrasted it with the failure on passing the Clean Air Initiative that same year due, most likely to a lack of public engagement. More recently, Nicholls (2017) documented the crucial role the general public had for the introduction of neonicotinoid regulations in Ontario. We expect that, just as in the national context, at the transboundary level, finding the specific levers to promote the support of the general public for conservation policies could strongly influence the design, reach, and ultimately, success of conservation policies.

The present paper seeks to improve the understanding of public preferences for transboundary conservation strategies for the monarch butterfly conservation. Its main goal is to determine whether inherent heterogeneity exists in public preferences for strategic-level characteristics of a recovery-management strategy that includes institutional leadership, international cooperation, and support for citizen science and research activities. It also aims to evidence the effect that providing a projection of the conservation program's success has on the overall willingness of people to support such a program. We used Discrete Choice Experiments (DCE) with a Latent Class Analysis (LC) to achieve such objectives.

#### METHODS

#### Sampling

The sampling area included the 35 eastern-most states of the US and Canada (Ontario), representing all areas where there is more than a 50% probability that monarch populations are present (Galindo-Leal et al., unpublished). Geographically speaking, the US-Mexican Border, parallel 49, meridian 102, and the east coast constitute the southern, northern, western and eastern limits of the sampling area, respectively (**Figure 1**). Based on the study objectives, we surveyed three different respondent samples: (1) main urban residents, (2) sample of monarch enthusiasts, and (3) modified urban residents' sample with a modified version of the survey to investigate the value of knowledge.

The surveying tool was delivered through the Survey Sampling International marketing company (Teel and Manfredo, 2010), targeting urban residents1 within the sampling area. Respondents were obtained from the panellists' database of the marketing company and were contacted directly by them based on our target demographics via email with an invitation link. The survey was sent in batches of 100, and only after analyzing their geographic and gender representativeness, the next batch of surveys was sent adjusting the target demographics to obtain a representative sample.

An invitation link was also sent through the Monarch Watch DPLEX mailing list2 , which contains subscribers, mostly citizen scientists, dedicated to the conservation of the monarch. This list is maintained by Monarch Watch, a non-profit organization hosted at the University of Kansas and dedicated to the monarch butterfly conservation (Lovett, 2018). We additionally isolated responses of individuals self-reported as engaged in monarch conservation activities from the main urban resident's sample and pooled their responses with the ones from the DPLEX mailing list to obtain a monarch enthusiasts' sample.

The use of an online survey through a marketing company, instead of a mailed or in-person survey, was due to the geographical and numeric extension of the sample. Online internet surveys have many other advantages, such as reduced cost and higher design flexibility. However, they also introduce new potential sources of bias that have to be accounted when analyzing the results such as a potential increase of self-selection processes (Olsen, 2009) and the risk of introducing "professional respondents" to the sample (Dennis, 2001).

The presence of "professional respondents" is one of the main risks associated with using marketing companies for delivering an online surveying tool. Such respondents tend to click through the survey without paying proper attention and potentially adding unwanted noise to the results (Dennis, 2001). To control for this, following Malhotra (2008), we removed individuals with a time-to-completion of two standard deviations away from the mean (individuals that averaged their responses in <5 s or above 22 min per question, including the choice experiment). This range was chosen since we could not find any evidence of primacy (Belson, 1966) or recency (Kalton et al., 1978) effects within those outliers. Most of the outliers did not complete the demographics section of the survey, for the few that did answer that section, we tested their demographics and attitudinal responses against the rest of the sample and did not find any significant differences (Malhotra, 2008). Lastly, an instructional manipulation check question was embedded (Oppenheimer et al., 2009) within the survey which read: "This question is intended to filter respondents that are not reading every question thoroughly. Please select the option 'Very Little' as your answer. Another question like this one will be placed further in the survey."

We sent out the survey to 5,750 people in Canada and the US from which we received 2,557 responses with an overall completion rate3 of 40.13%. The main sample included 916 individuals from Canada and 943 from the US, from which 302 self-reported as being monarch enthusiasts. Twentynine additional surveys were obtained through the Monarch Watch mailing list. We pooled the monarch enthusiasts from those two samples into a Monarch enthusiasts' group of 331 respondents. Finally, we sent 1,104 surveys with a variation of the DCE obtaining 625 completed surveys for the modified urban resident's sample with a completion rate of 56.51%. Completion rates varied according to the source of the respondents (marketing company = 46.4%, Monarch watch = 66.26%), and their nationality (Canada = 38.49%, US = 49.4%). Although our sample was mostly similar to the demographics of the target population, there were some significant differences, e.g., respondents from the sample with no high school diploma were 30% whereas the target population was 52%. Such demographic differences may be an effect of the way we defined urban residents in comparison to how it is stated in the census data, also, being

1We defined as Urban Resident a person that does not derive their main source of income from agriculture and owns a non-rural postal code.

2The survey was sent as an open link; however, we did not observe any duplicate IP addresses in the responses.

3Here and elsewhere, completion rate is defined as the number of surveys filled out and submitted divided by the number of surveys started.

this the reason to not balance the sample with an iterative proportional fitting or other raking procedure (Kolenikov, 2014); nevertheless, broad generalities to the target population can still be inferred. Respondents' demographics from the main urban residents and monarch enthusiasts' samples are summarized in **Tables 1**, **2**, respectively.

#### Survey Overview

Choice-experiments data were collected using a web-based survey conducted during November 2016 across Canada (Ontario) and the eastern US. The survey consisted of the following sections: (1) assessment of the individual's knowledge about the monarch, (2) video introduction for the survey and essential terminology, (3) choice experiment, (4) followup questions, (5) demographics, and (6) New Environmental Paradigm Scale (NEP) Statements. The survey also included questions on the allocation of resources and level of involvement of different organizations, which were not analyzed here but will be revisited in subsequent manuscripts.

The survey design and delivery were developed following Salant and Dillman (1994) and Dillman et al. (2014) design principles. Before giving any information about the monarch, we elicited the individual's knowledge of the monarch through three Likert-scale questions: (1) awareness of the monarch's decline, (2) level of concern about the current monarch's situation, and (3) awareness of the importance of milkweed for the monarch's survival and conservation. A short introductory video (2:32 min) followed explaining the purpose of the survey, the current decline of the monarch's population, and the definition of each DCE attribute. We used a video instead of text to avoid cognitive fatigue and to ensure respondents had a better understanding of the survey elements (Mendelson et al., 2017). Although we were unable to confirm that all respondents watched the video, they were unable to skip forward through the video to continue with the survey before it ended. The full survey and the video can be found in the **Supplementary Material**.

A demographics section was included after the DCE and, finally, the respondent was presented with the NEP Scale for the assessment of their environmental attitudes (Dunlap, 2008). The NEP scale consists of 15 environmentally-related statements to which the respondent must choose their level of agreement/disagreement. The totalled result is a score between 0 and 150, where the higher the score, the more ecologically oriented the mindset of the respondent (Dunlap and Van Liere, 1978).

# Discrete Choice Experiment (DCE)

The DCE is a stated preference valuation method that forces the respondent to make trade-offs between distinct levels and attributes ideally resembling the context in which individuals make real-life decisions. The DCE assumes that respondents' decisions follow the Random Utility Model, which states that an individual will strive to maximize utility while making choices (Manski, 1977). Under this assumption, it is possible to estimate the proportion of the sample, market share, that would choose any given program configuration (Landauer et al., 2012). By including a contribution attribute, the marginal economic value

*"Obs" is the observed percentage from the sample and "Exp" is the expected percentage of people based on census data (Statistics Canada, 2016 Census of Population; https:// www.census.gov/). The definition of urban resident of the US Census and Statistics Canada differed from ours. While their definition is based on population density, presence on urban clusters or urbanized areas, land use, distance, and population thresholds, our definition is based on main source of income and zip code. Note that some items do not add up to the total sample size due to missing data from incomplete responses.*

of the attributes can be estimated as well (Kuhfeld, 2006). The ability to explore hypothetical non-existent scenarios is another advantage of this method (Vega and Alpízar, 2011).

A DCE consists of a list of key characteristics, or attributes, describing an alternative. Each of these attributes has different values, or levels, defining the configuration of that alternative. Several alternatives, 2 or 3 at a time, are presented at the same time to respondents in a choice set. Then, respondents are asked to analyze and choose their preferred one from each choice set (Louviere et al., 2000). An orthogonal experimental design ensures that each choice set is presented to respondents enough times, allowing researchers to estimate respondents' preferences for the attributes and all the levels that defined those alternatives.

The DCE estimates the utility4 , or satisfaction that respondents derive from a choice, which, in this case, is a potential management scenario. Also, the DCE allows valuing not only the resource as a whole but also the incremental worth of its components—i.e., the marginal part-worth utility5 of its attributes (Birol et al., 2006) and their Relative Importance, RI (Vermunt and Magidson, 2005). DCEs are commonly used to forecast likely changes in behavior as a reaction to changed circumstances or to the hypothetical availability of certain goods (Louviere et al., 2000). The utility estimates from the DCE represent the utility that a level or unit of an attribute provides. When the attribute is categorical, this is measured as utility relative to the mean of the other levels from the same attribute. When the attribute is numerical, the interpretation of utility is on a "per unit" basis. The RI of an attribute, also known as Relative Maximum Effect, is the proportion of the overall utility explained by a change of one unit of that attribute when numeric, or from the difference between the least and most preferred levels of that attribute when categorical (Crouch and Louviere, 2004; Casini et al., 2016). The higher the RI value of one attribute, the

4Utility is defined as the weight of outcomes in making a decision (Ariely et al., 2003). It can also be explained as the level of short-term happiness derived from a specific material or immaterial good (Kimball and Willis, 2006). DCEs quantify utility by a mean-centered dimensionless value representing the preference associated with a particular level of an attribute compared with the reference level.

5Marginal part-worth utility is a measure of welfare that the respondent derives from a one-unit increment (all else being equal) of one attribute from the choice set (Steinke and Van Etten, 2017).

#### TABLE 2 | Demographics from the monarch enthusiasts' sample (*n* = 331).

more such attribute influences the preference of the respondent (Crouch and Louviere, 2004; Casini et al., 2016).

By including a contribution attribute within the experimental design, it is also possible to estimate a marginal willingness to pay (WTP) for each attribute (Kerr and Sharp, 2009). Taking advantage of this possibility, the estimates reported within this paper are in USD value. These estimates reflect the economic value of changing any attribute by one unit while leaving the remaining attributes fixed. The WTP presented here is a marginal WTP estimate on a per-unit basis from the baseline, which is different from the total WTP provided by other methods such as Contingent Valuation (Diffendorfer et al., 2013). Instead, the marginal WTP provided here denotes the difference in the contribution that the respondent would be willing to pay from the unweighted average of all the levels, for categorical variables (Daly et al., 2016). For numerical variables, it describes the difference of the respondent's WTP to increase one unit of a particular attribute while leaving the rest of the attributes fixed (Kerr and Sharp, 2009). This manuscript explores for the first time the marginal WTP for the monarch conservation.

#### Discrete Choice Experiment Design

We constructed choice alternatives describing potential management scenarios for the conservation of the monarch using a list of attributes that described a hypothetical ecological status of the monarch, and the strategic-level characteristics of a proposed conservation initiative. These attributes were refined using input from interviews with academics with expertise in human dimensions, conservation biology, or both. The final alternatives were made up of nine6 attributes, three of them as context and the other six as program attributes (**Table 3**). Values for the levels of each attribute were selected based on feedback from academics, two focus groups, and a pilot study with 200 respondents (100 Canadians and 100 US citizens) 2 months prior to the final version release. The first focus group (n = 8) consisted of experts on this method, with the primary objective of finding technical deficiencies. The second focus group (n = 13) was composed of graduate students of the authors' universities with differing levels of familiarity with choice experiments or the monarch and sought feedback about the size and complexity of the survey. Finally, the pilot study was directed to the same demographics as the target population of the main survey and sought to detect cognitive fatigue, such as positive WTP estimates, lack of

6The context attributes that appeared in the survey were "Trend" and "Area-Trend." The attribute "Area" did not appear in the survey, but it was used to calculate the "Area-Trend" attribute (which is an interaction between "Area" and "Trend"). Also, the "Payment Vehicle" and "Leader" are part of a single attribute in the experimental design but appear separately in the survey. See **Table 3** for details.

TABLE 3 | Attributes and levels used in the choice experiment exercise.

significance of the utility estimates, or extensive skipping of optional screens.

Each choice set (**Figure 2**) consisted of an ecological context scenario with three attributes, and three options: two alternative conservation programs, and one status quo option. Context attributes established the scenario under which the respondents would be making their choice (Tversky and Simonson, 1993; Haegeli et al., 2012). Here, the context attributes set a hypothetical situation of the overwintering colonies to investigate the change of respondents' preferences with the assumption that respondent preferences were context-dependent (Mazar et al., 2014). These context attributes remained the same for all options of the choice set and only changed between choice sets.

The program attributes included international allocation of funds, probability of success of the program, institution leading the program, monetary contribution to the described program, fund-raising mode, and percentage of funds dedicated to research and citizen science activities. These attributes varied their levels independently from each alternative so that the respondent could perceive a contrast between the options. The "status quo" option as a base alternative consisted of abstaining from contributing to any program and maintaining the current trend shown in the specific scenario. Most literature agrees that a base alternative has to be included to estimate the welfare change associated with the other alternatives (Bateman et al., 2004; Train, 2009). If the respondent chose the base alternative in any of the presented choice sets, they were asked to provide a rationale for their choice.

The experimental design for the main urban residents' survey was a 46 × 8 1 orthogonal fractional factorial design with two of those factors entered as context variables. For the modified urban

resident's survey, a new design with the same characteristics was generated but with one factor removed (45 × 8 1 ). Both designs were generated with the SAS "%MktEx" Macro (Kuhfeld, 2001) and had a D-efficiency of 100% as a measure of the design's goodness (efficiency), and orthogonality (Kuhfeld et al., 1994).

#### Statistical Analyses

All the numeric levels were standardized and centered before analyzing the DCE model. The data were analyzed using conditional logit and latent class regression with Latent Gold 3.0 software (Vermunt and Magidson, 2005), obtaining Relative Importance (RI), latent class segmentation outputs, and model performance metrics.

Latent Class Analysis (LC) was used to identify and segment heterogeneity in utility estimates among urban residents. The LC assumes that the sample constitutes a finite number of groups of individuals, also known as classes, with relatively similar preferences within their group and considerably different from each other (Birol et al., 2006). Random Parameters Logit can also identify the heterogeneity of preferences within a sample (McConnell and Tseng, 1999); however, Random Parameters Logit elicits the individual differences amongst the sample rather than grouping them (as LC does). The latter scale of analysis is considered more convenient for the design of management strategies (Boxall and Adamowicz, 2002).

The non-significantly different attributes across classes in preliminary models were constrained to be the same across all classes to prioritize the delineation of classes by the most highly variable attributes (**Table 5**). That model restriction reduced the number of parameters and improved the fit of the model (Vermunt and Magidson, 2000).

Embedding a DCE within a comprehensive survey allows descriptive data, as covariates or predictors, to define individuals by linking these with their preferences. Covariates are a posteriori explanatory variables that describe class membership and can inform the policymaker about which demographic strata can be targeted with specific actions (Boxall and Adamowicz, 2002). Covariates included in the model were the pre-survey knowledge about milkweed and the monarch's status, whether the respondent was engaged in any ecological/citizen science activity, and the age group of the respondent.

Alternatively, predictors are characteristics of the choice replication or the person and have the same value across alternatives. Predictors are part of the regression model, just like attributes, and are therefore considered a priori explanatory variables (Vermunt and Magidson, 2005). As a result, covariates can predict class membership, whereas predictors contribute to its creation. Here, the model included the level of concern about the monarch's situation as a predictor.

For the three respondent samples (main urban residents, monarch enthusiasts, and the modified sample of urban residents), we also conducted a Multinomial Logit Model (MNL) analysis to obtain a one-class model for each. These types of models are suitable for observing the main trends of the sample without accounting for heterogeneity. The MNL was used to compare the three samples and qualitatively detect any differences between the general preferences of people engaged—or not—in ecological activities (urban residents vs. monarch enthusiast's sample), or between people provided with an expected probability of success—or not—of the proposed program (main urban residents' vs. modified urban residents' sample).

To control for the uneven spacing of some of the numericvariable attributes and to achieve more interpretable results, we linearized all our numeric attributes (Kohlhardt et al., 2018). All the categorical attributes were effects coded for the interpretation of their estimates (Daly et al., 2016). Numeric data were analyzed with one-way ANOVA and a post hoc Tukey's Honest Significant Tests. For categorical data, a Pearson's chi-squared test was used. All statistical treatments were done with JMP 13 (SAS Institute Inc, 2016), and R 3.51 (R Core Team, 2013) was used to plot the results.

# RESULTS

### Latent Class Analysis of Main Urban Residents

#### Description of Classes

Preliminary models with different number of classes, covariates, predictors, and constraints (**Table 4**) were defined and evaluated using Bayesian Information Criteria (Burnham and Anderson, 2004). We also built a preliminary 2-known-class model based on nationality, and no significant differences were found between the classes regarding their preferences for the attributes presented; we pooled the data as a result. The final model was a three-class model with significantly different preferences for the geographical allocation of the resources, sensitivity toward the allocation of funds across classes, and the Alternative Specific Constant (ASC), which can be described as the utility derived from selecting any choice different from the status quo without accounting for the specific levels of the rest of the attributes. Each class was labeled based on those differences as "Eager," "Pro-Nation," and "Opinionated." The final model had the "Leader," and "Area" attributes constrained between class "Eager" and "Pro-Nation," "Research" across "Eager" and "Opinionated," and "Trend" across the three classes (**Table 5**). The "Eager" group was the largest, making up 72.4% of the overall sample. The "Pro-Nation" class was second in size (18.37%) and "Opinionated" was the smallest (9.23%)7 .

Individuals from the class "Eager" showed a large estimate for the ASC, which represents a strong motivation to support conservation initiatives regardless of the configuration of the choice set (**Table 6**). In contrast, the other two classes denoted an unwillingness to participate in any management program. People from the "Pro-Nation" class strongly based their decisions on the allocation of funds across countries. When the choice task indicated that the allocation of funds would favor the respondent's country of residence, their utility markedly rose. In contrast, when funds were allocated only to Mexico or to the "other country," i.e., the US for Canadians, or Canada for US citizens, their utility considerably decreased in comparison to the other two classes. This class had a difference between the highest and lowest valued estimates 34.4% larger than that of "Eager." Finally, the third and smallest class was labeled "Opinionated" due to the large estimates associated with the leading institution, resources allocation, and especially the economic contribution. This class also had the most negative ASC, implying that they are the most reluctant to participate in any management program.

Respondents in the "Eager" group displayed the highest NEP score, indicating that these individuals possess largely pro-environmental attitudes. They tended to be younger and had a higher level of education, where 82.6% obtained at least a bachelor's degree, furthermore, 17.1% had a graduate certificate. Their income level was also higher than the other two classes, where 62.4% of the group earned at least \$50,000 per annum and also had the largest household size. The "Eager" class had the most considerable share of people contributing to ecologically oriented NGOs and actively participating in ecological conservation meetings, protests, and lectures. However, 58.5% of the people participating in those activities did not contribute economically to any ecologically oriented NGO (**Tables 7**, **8**).

The "Pro-Nation" and "Opinionated" classes were similar in attitudinal preferences and demographics, except in the percentage of individuals contributing to environmentally related activities and in age. Also, a higher proportion of the "Pro-Nation" class contributed to ecologically oriented organizations in comparison with people from the "Opinionated" class.

Only the level of concern about the monarch situation was included as a predictor of choice in the definition of the model as it significantly improved model fit. The overall utility estimates for "Eager" and "Pro-Nation," which add up to 91% of the overall sample, were positively affected when respondents had a higher level of concern about the monarch's situation. The reaction of "Opinionated" was counterintuitive, where its overall utility was negatively affected by an increase in their level of concern.

#### Context Attributes' Estimates

Further interpretation of the classes can be made by considering the attributes themselves and their levels (for a full list of estimates refer to **Table 6** and **Figure 3**). Respondents reacted to the percentage change of the overwintering monarch colonies' size over the last 5 years, in relation to the current area, similarly negative across the three classes, and all respondents' interest in supporting management programs decreased when the monarch population trend increased.

For the current area of the overwintering colonies, the "Eager" and "Pro-Nation" classes reacted similarly. They both were significantly affected negatively by the increase of the area of the overwintering colonies, i.e., their interest in supporting management programs decreased when the current colony population was higher. For the "Opinionated" class, we found the opposite effect. All the "Area" estimates were significant only at the 10% level.

7A LC provides the posterior probability that an individual belongs to a certain class (McCutcheon, 1987). We assumed that the class membership of a respondent was dictated by the class that gave them the highest posterior probability (Pacifico and Yoo, 2013).

TABLE 4 | Model selection for the main urban resident sample (*n* = 1,859).

*The base models have no restrictions, whereas subsequent models (with the same number of classes) are variations of that first model with different combinations of constraints, covariates, and predictors. Model selection was based on the best (lowest) BIC and smaller classification error (Class. Err).* \*\*\**1% significance level with two-tailed tests.*

TABLE 5 | Definition of constraints for the 3-latent class model of the main urban resident's sample (*n* = 1,859).

*Classes with similar preferences on preliminary models for a particular attribute were assumed to be the same in the final model, so other attributes with higher variance could drive the splitting of classes. Classes with the same letter denote that they have the same estimate for that specific attribute.*

As described in **Table 3**, the "Change" attribute was an interaction attribute between the overwintering colonies' Trend and Area. Respondents from the "Eager" class derived a positive utility from this attribute, i.e., the more substantial the increase, the higher the interest in supporting management programs. "Pro-Nation" respondents derived a negative utility, and "Opinionated" respondents were not significantly affected by this attribute.

#### Program Attributes' Estimates

The estimates for the institution leading the program were equal across "Eager" and "Pro-Nation." For these two classes, International NGOs and Educational institutions were significantly positive. Alternatively, "Opinionated" respondents showed a preference for local NGOs as leaders of the program. In all cases, the least preferred leading institution was the federal government.

When the allocation of resources was distributed to the respondent's own country, the utility estimates were the highest for the "Pro-Nation" class. The utility of the "Pro-Nation" and "Eager" classes became negative when either Mexico or the counterpart North American country were the receivers of those resources. Respondents from the "Opinionated" class were only significantly negatively affected when the counterpart country was the beneficiary of the resources. When the resources were distributed equitably across the three countries, the attribute's estimates were the highest for the "Eager" and "Opinionated" classes.

Regarding the percentage of funds dedicated to research and citizen science activities, the utility was similarly negative across the "Eager" and "Opinionated" classes and not significant for "Pro-Nation." For the probability of reaching the conservation goal of a minimum size of 6 ha for the overwintering colonies in 10 years, the utility estimates for "Eager" and "Pro-Nation" were significant and positive but being the first double than the latter; "Opinionated" had no significant preferences.

Finally, the attribute asking for the amount of money that respondents would be willing to donate for supporting the selected management strategy was negative and highly significant for all three classes. However, the "Opinionated" class estimate was almost double than that of "Pro-Nation" and almost 10-fold than that of "Eager" respondents.

#### Monarch Enthusiast's Estimates

The monarch enthusiasts sample (n = 331) consisted of individuals from the main urban residents' sample that selfreported as being monarch enthusiasts, and people from the DPLEX Monarch Watch mailing list. The primary objective of this sample was to identify differences between this group TABLE 6 | Latent class (3 classes) estimates and Marginal Willingness to Pay (mWTP) for the main urban residents' sample.

Monarch Conservation Through the Social Lens

*See text for the definition of Relative Importance (RI).* \*\*\**1% significance level,* \*\**5% significance level,* \**10% significance level with two-tailed tests. The Attribute "Area-Trend" is an interaction attribute between "Area" and "Trend".*

*We tested the differences between classes with a one-way ANOVA test and post hoc Tukey's test and the class "Eager" was significantly different (P* < *0.001) to the other two classes, which were no significantly different between each other.*

\*\*\**1% significance level,* \*\**5% significance level.*

and the main urban residents' sample. For this sample, the estimates obtained from the MNL (**Table 9**) closely resembled the estimates from the "Eager" class of the main urban residents' sample with the following exceptions: this sample showed a positive utility for the type of institution leading the program only when it was an educational institution. The remaining levels did not significantly affect the monarch enthusiasts' choice, unlike "Eager" respondents that had significantly positive estimates for both international NGOs and educational institutions. Also, while in the main urban residents' sample each of the classes had significant estimates for at least one of the context attributes ("Area," "Trend," or "Area-Trend"), the monarch enthusiasts did not exhibit significant preferences for any of them. Lastly, the estimate for the monetary contribution to support the program was negative (just as with the main urban residents' sample), but the value of the attribute was noticeably smaller in magnitude.

The ANOVA test shows that the demographics of this sample were significantly different from the main sample and each one of the three classes. A more substantial proportion of monarch enthusiasts were engaged in ecologically-related activities (p < 0.001) as well as the percentage of them who contributed to ecologically-oriented NGOs (p < 0.001). The percentage of enthusiasts that were Canadian was significantly lower than the share of Canadians from the urban residents' sample (p < 0.001). Respondents from the monarch enthusiasts' sample also had a higher level of education (p < 0.001), although the income level was not significantly different. Unlike the main sample that had more females than males, the citizen scientists' sample had a significantly higher proportion of males (p = 0.007). Finally, the average age of the enthusiasts' sample averaged significantly (p < 0.001) lower than the urban resident's sample (**Table 2**).

#### Modified Urban Resident's Estimates (Success Omitted)

The attribute most influenced by the inclusion/exclusion of a success probability was the percentage of resources dedicated to research. When included, the utility estimate of contributing funds to research was negative, i.e., respondents from the main urban sample were less willing to provide funds toward research when the program specified an expected success. Conversely, with the removal of this attribute, the estimate for research became positive; i.e., contribution-support increased in the absence of knowing success. However, amongst the respondents from the modified sample, the ASC value was negative, denoting a decrease of willingness to support conservation measures overall (**Table 9** and **Figure 4**).

#### Willingness to Pay

The marginal willingness to pay (mWTP) for each of the attributes was calculated and is shown in **Table 6**. The mWTP is defined as the difference in the contribution that the respondent would be willing to pay from the mean of all the levels, for categorical variables (Daly et al., 2016) and the difference of the respondent's WTP to increase one unit of a particular attribute while leaving the rest of the attributes fixed (Kerr and Sharp, 2009). Finally, the total WTP to support a conservation program for the monarch was also estimated. The WTP was contingent on the configuration of the program8 and followed the utility estimates described in previous sections and, based on the

8The configuration of the "Best program" was defined as a program with the levels that obtained the higher utility estimate for each of the categorical attributes, with 90% success, and 20% of funds dedicated to research. Conversely the "Worst program" used the levels with lower utility, had 70% success, and also dedicated 20% of funds to research.

current area and trend of the overwintering colonies (Rendón-Salinas et al., 2019), it ranged between \$100.41 and \$141.01 for the worst and best program configurations, respectively. When analyzing each of the classes, the average WTP was \$161.76, \$76.85, and \$-5.04 for the classes "Eager," "Pro-Nation," and "Opinionated," respectively.

# DISCUSSION

The monarch butterfly is an iconic species for people from the US, Canada, and Mexico alike (Guiney and Oberhauser, 2008). As such, its conservation provides an excellent opportunity to find common points of interest and strengthen, or create, institutions of tri-national cooperation for the recovery of the monarch and other transboundary migratory species as well (Lopez-Hoffmann et al., 2009). Moreover, the monarch's plight has mobilized a considerable number of urban residents across the three countries to participate in habitat restoration and citizen science efforts to protect it (Ries and Oberhauser, 2015). The role of small habitat providers and citizen scientists that urban residents play in this context also extend to conservation-policy support. Conservation practitioners should strive to find the most effective ways to funnel this potential capacity, with that objective, and this study aimed to determine urban-resident preferences toward strategiclevel characteristics of a management strategy for monarch conservation that would generate the highest amount of support from urban residents.