Present-Day Distribution and Potential Spread of the Invasive Green Alga Avrainvillea amadelpha Around the Main Hawaiian Islands

Algal assemblages are critical components of marine ecosystems from the intertidal to mesophotic depths; they act as primary producers, nutrient cyclers


INTRODUCTION
Disturbances can dramatically alter the community composition of coral reef ecosystems. The replacement of coral cover by macroalgae is known as a "phase shift" (Littler and Littler, 1984;MacManus and Polsenberg, 2004;Bruno et al., 2009), which can occur in response to pressures such as storm events (Rogers and Miller, 2006), coral disease outbreaks (Aronson and Precht, 2001;Porter et al., 2001), nutrient input (McCook, 1999), or removal of herbivores (Hughes et al., 2007). Often, these stressors occur simultaneously or in sequence as one stress event lowers the threshold of coral resilience, thereby magnifying total coral mortality (Brandt et al., 2013;Redding et al., 2013;Casey et al., 2014;Vega-Thurber et al., 2014). Shallow coral reef ecosystems endure common reef stresses, including overuse and physical degradation (i.e., trampling), overextraction, chemical degradation from land-based runoff, effluent, or pollutants, and invasion by harmful bacteria, viruses, or competitors (Rodgers et al., 2003;Friedlander et al., 2017;Stamoulis et al., 2017;Takesue and Storlazzi, 2017). Invasive algae are an ecological and economic concern in Hawaii ( Smith et al., 2002;Schaffelke et al., 2006;Westbrook et al., 2015). Invasive species can become opportunistic when some environmental change facilitates their spread and allows them to outcompete native species. In Hawaii, eutrophication due to runoff likely influenced the spread of the non-native and invasive red alga Hypnea musciformis (Smith et al., 2002) and blooms of the native green alga Dictyosphaeria cavernosa (Stimson et al., 1996).
One example of a particularly troublesome competitor that has garnered statewide attention in the past several decades is the invasive green macroalga Avrainvillea amadelpha. This macroalga, also known as "leather mudweed" due to its thick, paddle-like blades (Figure 1a), was first discovered in Hawaii in 1981 (Brostoff, 1989) and has been observed from the intertidal to the mesophotic zone (30-90 m) (Spalding, 2012;Pyle et al., 2016;Cox et al., 2017). Its origin and phylogeny remain uncertain, but it is not present in historical records of shallow marine habitats around the main Hawaiian Islands (MHI) (Brostoff, 1989). To date, this highly successful macroalga has spread across Oahu (Figure 2) and has been observed sporadically around Kauai (Smith et al., 2002;Wade and Sherwood, 2018). In the mesophotic zone, video records show A. amadelpha beds competing with meadows of an undescribed Udotea sp. (Spalding, 2012) (Figures 1b,c).
Avrainvillea amadelpha is of particular concern in Hawaii because of its weedy characteristics, such as benthos alteration (Martinez et al., 2009;Sansone et al., 2017), tolerance to environmental extremes (e.g., light availability, exposure, etc.), and unpalatability to native herbivores (Van Heukelem, 2016). The alga is also presumed to have successful vegetative propagation via fragmentation and holdfast siphon viability like many other siphonous green algae (Walters and Smith, 1994). This alga may also change the community diversity and structure across invaded reefs. Langston and Spalding (2017) recorded a higher abundance of fishes above open sand versus within A. amadelpha canopy; Longenecker et al. (2011) determined that removal of A. amadelpha returned invertebrate communities to compositions observed in unaffected regions. All of these traits combine to make A. amadelpha a formidable competitive presence across reefs.
The state has spent over $2 million to remove 1.32 million kilograms of A. amadelpha across a 90,000 m 2 swath of Maunalua Bay, Oahu (Kittinger et al., 2016). Much effort has been invested in observing and recording A. amadelpha invasions, as well as stymying the progress of ongoing invasions (Hawaii Department of Aquatic Resources, unpublished data). However, to date, no projections exist regarding the potential spread of this invasive macroalga across shallow reef systems. When used as a management tool, habitat suitability models may provide critical information about the current or probable future distribution of opportunistic competitors (Guisan et al., 2017). Habitat suitability models can also improve understanding of dispersal methods and provide early warning about regions that are particularly vulnerable, but not yet invaded. This work sheds light on areas of particular concern across the state and presents a transferable, quantitative framework that resource managers may utilize as one aspect of monitoring for invasive species.

Observational Data
Our observational data were sourced from video obtained during submersible and remotely operated vehicle (ROV) dives and shore-based or snorkeling surveys (Figure 2 and Table 1 Table S1). We collected 3,796 new snapshots of the benthos by pausing each video track in 30 s intervals. From these, 2,441 photo stills were excluded from our mesophotic data set due to blurriness, stationarity of the vessel, or positioning outside the depth maximum (>90 m). We collected shallow (∼5-10 m) subtidal observations of A. amadelpha during multiple snorkeling and free-diving expeditions conducted from July 2015 to November 2017. We combined these data with two existing subtidal survey data sets, which included data from 2015 to 2017 provided by the Hawaii Department of Aquatic Resources. Intertidal shorebased survey data were acquired through transect sampling at low tides (0-to −0.15 m) by the Our Project in Hawaii's Intertidal school-based monitoring program from March to June 2017 1 .
We removed 8,040 observations from the raw data set of 33,286 points due to overlap within pixels or failure to fall within 1 https://opihi.crdg.hawaii.edu our depth range (0-90 m). All shoreline survey point data were inspected for geographical precision and manually moved to the closest classified marine pixel if needed. All observational data were combined into one data set spanning depths from 0 to 90 m and covering all coastlines across Oahu. Data were resampled to fit a grid size of 250 m × 250 m and duplicate observations within the same grid cell were removed, resulting in a final data set of 276 observations of A. amadelpha presence (207) or absence (69). We note that the mesophotic video data used to build this model were gathered approximately one decade before the intertidal and shallow observations. This gap in data collection is due to the costly nature of submersible dives, and we stress that these mesophotic data comprise the most extensive, up-to-date catalog of A. amadelpha occurrence below the shallow (>30 m) zone.

Environmental Data
We considered 15 environmental predictor variables for inclusion in our model (Supplementary Table S2). We represented these predictor variables, or covariates, as 250 m raster grid layers across the study domain, which extended approximately 1-2 km from the coast and encompassed depths from 0 to 90 m. We broadly categorized predictors into three classifications: topographic, oceanographic, and anthropogenic. All covariates were selected for consideration based on their representation in the literature as potentially influential in determining the distribution of A. amadelpha. These covariates may have a direct influence on the metabolic constraints and, therefore, the occurrence of A. amadelpha. Though A. amadelpha has, to date, only been observed around Oahu and sporadically in Kauai, we applied our model to the extended domain of the MHI (Hawaii, Maui, Lanai, Molokai, Kahoolawe, Oahu, Kauai, and Niihau). We expect our model to identify regions that are most suitable for A. amadelpha habitation, which we can pinpoint as at-risk areas that managers may prioritize during monitoring efforts.
Considered topographic variables included seafloor slope, aspect (i.e., compass direction of seafloor), rugosity (i.e., roughness), and bathymetric position index, which is a metric that combines the bathymetry and seafloor slope within a 5 km neighborhood of a given point. Considered oceanographic variables included seafloor current velocity (annual mean), seafloor temperature (annual mean), surface chlorophyll a concentration (annual mean), and surface turbidity (annual mean). Lastly, considered anthropogenic variables included indices describing fishing pressure, nearshore coastal development, anomalous sea surface warming, annual estimated ship crossings, nutrient and sediment inputs, and proximity to human population (Falinski, 2016;Lecky, 2016;Wedding et al., 2018). Fishing pressure and coastal development were scaled from 0 to 1 by dividing a given value by the maximum value in the respective categories. Anomalous sea surface warming, represented as degree heating weeks (DHW), is a metric which couples the duration and intensity of abnormally warm ocean surface temperatures (Liu et al., 2013). Supplementary Table S2 includes further details about considered environmental covariates.
We first examined the spatial coverage of all covariate layers. We excluded the layers displaying proximity to human population and sediment influx from inclusion in our analyses due to their insufficient spatial overlap with mapped observational data and patchy coverage across the considered depth range surrounding the remaining MHI. We chose not to extrapolate these two layers based on our uncertainty about introducing significant error in our predictor data. The data sets detailing nutrient flux and coastal development were more spatially extensive, allowing us to extrapolate the closest offshore values for each variable and apply these values to approximately 5% of offshore observational data points outside the coverage of the layers.
Collinearity of predictor variables may preclude suitable model fitting (i.e., preventing the model from converging or contributing to high standard errors in coefficient estimates). We created a correlation scatterplot to examine the relationships between all covariates and the response variable (i.e., A. amadelpha occurrence) (Supplementary Figure S1). Per Tabachnick and Fidell (1994) and Dancey and Reidy (2004), we removed covariates that met one or more of the following conditions: lack of obvious visual correlation with the response variable, correlation of ≥70% with other predictors that were more strongly correlated with the response variable, and/or lack of a visible or quantitative correlation with the response variable. We included seven covariates in our initial model: seafloor slope, annual mean seafloor current velocity, maximum DHW, nutrient flux, shipping traffic, coastal development, and annual fishing catch.

Model Development
Boosted regression tree (BRT) models fit an ensemble of statistical models by combining the power of boosting (a technique that blends multiple models to improve predictions) and regression trees (models that identify links between response variables and covariates via recursive binary splits; Elith et al., 2008). Three key parameters determine the type and shape of the decision trees that compose the model: bag fraction, learning rate, and tree complexity. To start, we retained stochasticity in our model by setting the subsampling rate, or the bag fraction, equal to 50%. This setting indicates that when training each new tree, the model randomly selects 50% of the data from the training data set. The remaining fraction of the training data not used in the bag fraction is used during that iteration for cross-validation and evaluation of the model. Next, we specified a slow learning rate of 0.0025, which indicates the contribution or weight of each tree in the overall model. The selection of a slow learning rate minimized the contributions of each tree to the final model, minimized the number of very small, random subsets of data selected for fitting each tree, and decreased the chance of fitting abnormal trees. Lastly, we set our tree complexity parameter to 5, which permits the model to create decision trees with 5 nodes. The number of nodes corresponds to the maximum number of predictor interactions the model may identify; we methodically varied all of these parameters to find the optimal combination that would explain the most variance within the data (Elith et al., 2008).
Loss functions are used to calculate deviance, which is a goodness-of-fit statistic indicating how well the model captures variance within the data. Minimized deviance of a model indicates high explanatory, or predictive, power. To improve the performance of our model, we used binomial deviance, which is suggested in the literature due to its robustness when used with noisy or potentially misclassified data (Hastie et al., 2001;Elith et al., 2008). While we do not anticipate the misclassification of a large portion of our data, we chose to use this loss function given that our observational data sets originated from different sources, including researchers, school groups, and state workers. The loss function allowed us to calculate changes in residual deviance, or goodness-of-fit, as we systematically varied key model parameters, which guided our final selection of parameter settings .
We determined the relative importance of each predictor variable using formulae within the "gbm" library (Friedman, 2001;Friedman and Meulman, 2003). These values are determined by the number of times each variable is selected in a split and weighted by the subsequent model improvement based on that split (weights are squared to give higher value to predictors that improve the model more considerably). Predictor contributions were scaled such that the sum equated to 100, with higher values indicating greater influence of that predictor on the response variable, i.e., occurrence of A. amadelpha.
We allocated 70% of our data to develop and train our model. The remaining 30% of our data was set aside for model testing. Per Elith et al. (2008), we initially set the number of trees (nt) to nt0 = 50 and ran 10 unique BRTs. We subsequently increased the value of nt by 10 (i.e., nt1 = 60, nt2 = 70, and so on) and ran new models until we found a combination that produced higher average mean predictive performance and lower average standard errors than the previous set of models. Our lr was kept very low (lr = 0.0025) to ensure precise estimation of the starting number of trees (nt = 60).

Model Evaluation
Our initial model (193 observations, 7 predictors) fit 1,680 trees. We trained the model using k-fold cross-validation (k = 10, indicating 10 iterations during which half the training data was randomly selected to develop the model). We evaluated model predictive performance based on the cross-validated area under the curve (AUC) statistic (Hosmer and Lemeshow, 2004), percent of deviance explained by the model, and the standard error of the model. Performance metrics were recorded when we applied the model to the validation data set during training and when we applied the model to the withheld testing data set after model finalization.
We constructed a dotplot to compare known training data set observations to present-day model-predicted occurrence probabilities (Figure 3). The categorical means-indicated by red dots-differ between the two groups; for those pixels where A. amadelpha was not observed, the mean probability was 0.45, while a mean of 0.95 was calculated for pixels known to contain A. amadelpha. The standard deviation of the presence category was much smaller than that of the absence category (i.e., 0.075 and 0.4, respectively), suggesting that the model is prone to false positives. Additionally, we plotted a curve of probabilities predicted by logistic regression; this curve indicated that probabilities below 50% indicate regions fairly unsuitable for A. amadelpha colonization. Risk increases exponentially beyond this probability, with the highest risk of invasion predicted for regions with habitat suitability values >0.95 (Supplementary Figure S2).
To construct confidence intervals for a given predictor, we bootstrapped the model 100 times, which required us to resample the training data with replacement (note that this differs from the specification of bag fraction during model training, which necessitates subsampling without replacement) (Efron, 1979;Hillis and Bull, 1993). Next, we simulated dropping up to five different variables to evaluate change in predictive deviance; due to the negligible (−1%) change in deviance, we chose not to simplify the model, and retained all 7 predictors for further analyses. We applied our model to the withheld test data (83 observations) and recorded all evaluation metrics. Finally, we tested for spatial autocorrelation within the data to ensure the independence of the response variable (Moran's I coefficient = 0.088).
Degree heating weeks is a metric that measures cumulative regional heat stress, quantified as the number of weeks in a 12 weeks period during which sea surface temperatures exceed a 1 • C increase over the summer mean; a threshold of 4 DHW Frontiers in Marine Science | www.frontiersin.org has been identified as a trigger for probable coral bleaching events (Liu et al., 2008). We applied our model to forecast the change in invasion vulnerability across Hawaii's coastline given a 25% increase in annual DHW. All other predictors were held constant to isolate the effect of this change. This conservative change was based on the dire projections communicated in the fifth Intergovernmental Panel on Climate Change report, which predicts a global average sea surface temperature change of up to +2.0 • C relative to the mean temperatures observed from 1986 to 2005 (IPCC, 2014). A temperature increase of this magnitude will catastrophically destabilize marine ecosystems on a global scale, making them more vulnerable to colonization by invasive species such as A. amadelpha. Regional increases of up to 4 DHW are no longer unimaginably improbable given these grave predictions.
All code is available in the Supplementary Material (Av_code.R). Data will be archived at Dryad. We used R statistical software (R Core Team, 2017) and ESRI ArcGIS v.10.4 (ESRI, 2017) to perform all data analyses. Our BRTs were specified using the "dismo" package version 1.1-4 (Hijmans et al., 2017) and visualized using the "ggBRT" package (Jouffray et al., 2019).

Model Performance
Our initial area under the curve (AUC) statistic (0.99) indicated possible overfitting of the model to the training data. However, our cross-validated AUC (0.94) indicated high predictive accuracy during this initial training stage. Our cross-validated standard error was 2.1% and the model explained 62% of the validation data deviance. When applied to the testing data set, AUC slightly increased (0.96), with an average of 87% of presences and absences identified correctly by the model. The model explained 63% of the deviance within the testing data set, demonstrating its predictive capability when applied to novel data. Further model performance details may be found in Table 2.

Influential Environmental Covariates
Our model indicated that seven variables that exhibit some influence on the probability of A. amadelpha occurrence across Hawaii. These predictor variables include percentage of shoreline development (30%), maximum DHW (26%), fishing pressure (14%), mean current velocity at the seafloor (11%), shipping traffic (9%), seafloor slope (6%), and nutrient flux (3%). Our model scaled the relative influence of each predictor based on the total number of covariates in the final model to identify the most important predictor variables; this delineation is indicated by a dashed red line (Figure 4). The model identified shoreline development, thermal intensity and exposure (i.e., DHW), and fishing pressure as the most influential predictors of A. amadelpha occurrence.
To visualize the effect of each predictor on the probability of A. amadelpha occurrence, we bootstrapped our predictions and generated partial dependence plots with confidence intervals (Figure 5). Our model indicated the importance of coastal development in determining the probability of A. amadelpha occurrence. Our plot shows that regions of the coastline that experienced a minimal amount of recent development (5-10% from 2005 to 2010) may be more prone to colonization. Additionally, we note higher vulnerability for invasion in flat regions exposed to minimal seafloor current, at least 400 ship crossings annually, moderate fishing pressure, and at least 5 DHW. Nutrient flux was identified as the least influential covariate, and our plot shows only a mild increase in regional vulnerability when exposed to greater nutrient inputs. The confidence intervals show some model uncertainty regarding higher levels of fishing or shoreline development, but are tightest around certain spikes in the data, highlighting the important effect of certain thresholds for these predictors in relation to increased vulnerability for A. amadelpha incursion.

Areas Highlighted as Hotspots
Areas around Oahu, identified by our model as most susceptible to A. amadelpha invasion include north Oahu, eastern Oahu, and much of the Honolulu metro area (Figure 6, inset A). Patches of these coastlines have been surveyed and are already invaded by A. amadelpha (Figure 2). We extended our model to identify environmentally suitable regions at risk for invasion across the surrounding MHI. Our results indicate that managers should concentrate surveys for A. amadelpha along the southern and western coasts of Maui Nui, western and eastern Hawaii, northern and eastern Kauai, and south Niihau (Figure 6). Model predictions extend approximately 1-2 km offshore and encompass the shoreline to the mesophotic zone (0-90 m). The model indicated that approximately one quarter (423 km 2 ) of the modeled area (1,691 km 2 ) is highly suitable A. amadelpha habitat (θ ≥ 0.89, or third quantile cutoff value).

DISCUSSION
This is the first study to predictively model and map the presentday and potential distribution of invasive leather mudweed across Oahu and extrapolate results to surrounding islands in the event of potential future climate conditions. We highlight this approach as one possible management tool for prioritizing management and conservation efforts and in the face of changing anthropogenic and environmental drivers and future climate conditions.

Marine Habitat Susceptibility: Oceanographic and Topographic Factors
Shallow (0-30 m) marine regions identified as most susceptible to A. amadelpha colonization are relatively flat and experience minimal annual mean current flux. In other words: regions that meet or exceed the mean predicted present day occurrence probability of 63%, annual benthic current velocity = 0.01 m/s and mean slope = 1.4 • . Littler et al. (2004) hypothesized that A. amadelpha mounds in Belize maximize productivity when they occur in protected, calm, and shallow embayments. We expect that relatively low seafloor current flux may indicate a lower mass of sunlight-blocking suspended matter. Peyton (2009) noted that A. amadelpha hosts epiphytes in its canopy that may act as "sunscreen" during periods of high light exposure, which suggests some level of high light intolerance in this alga.
In the mesophotic zone, A. amadelpha is most likely to colonize similarly flat, calm regions (annual seafloor current velocity = 0.01 m/s, mean slope = 3.7 • for regions that meet or exceed the mean occurrence probability of 77%). Our analyses do not clearly reveal any distinct differences in the habitat preferences of A. amadelpha across depth gradients relating specifically to oceanographic drivers. If A. amadelpha continues to proliferate across the MHI, additional observations to the data set and subsequent updating of the model may tease out depthspecific differences in the settlement patterns of this invader.

Marine Habitat Susceptibility: Anthropogenic Factors
Our model suggested that anthropogenic stresses on land and in the sea may increase the possibility of A. amadelpha invasion. Regions adjacent to moderate, recent coastal development (about 10% more development from 2005 to 2010) or exposed to at least 400 ship crossings annually were identified as more susceptible to colonization by the invasive alga. Shoreline development did not appear to drive colonization patterns beyond the threshold of 10% relative construction. This metric includes sediment dispersal in the calculation of the relative value, and we theorize that sediment influx may negatively impact A. amadelpha metabolism beyond some minimal to moderate level, which may be why the trend flattens beyond the aforementioned threshold. Shipping traffic also shows some effect on habitat preference up to a certain point, but this relationship is clearly not linear, and further study may be warranted which contrasts ports known to contain A. amadelpha with those that are presently uninvaded.
The model also suggests that regions exposed to moderate fishing pressure are more prone to A. amadelpha invasion. Prior survey data (K. Peyton, Hawaii Conservation Conference; unpublished data) shows some correlation with the removal of herbivorous fish and urchins and the subsequent appearance of A. amadelpha. Though phase shifts may occur across overfished reefs (e.g., Hughes et al., 2007), researchers are uncertain about the role of grazing pressure in controlling the spread of an unpalatable invader like A. amadelpha (Van Heukelem, 2016). However, fishing catch was one of the factors that best explains the presence of A. amadelpha, which argues that the role of reduced herbivory must also be considered. Although herbivorous coral reef fish populations around Oahu are generally low (Williams et al., 2016), and the influence of reduced herbivory on the increase of A. amadelpha is unknown, the correlation between fishing pressure and spread of this alien invasive merits further study. We note that the partial dependency plot for fishing pressure (Figure 5) shows some reduction in probability of A. amadelpha occurrence for those regions exposed most intensely to fishing effort. We posit that this drop in occurrence probability as fishing stress intensifies may be akin to the same "Goldilocks"type trend we observe in relation to coastal development and A. amadelpha occurrence: a little bit of stress-not too much, not too little-is just right for opening the door to this invader.

Invasion Risk Under a Changing Climate
Our model also indicated some influence of DHW on the susceptibility of habitat to A. amadelpha colonization. Eakin et al. (2010) note that thresholds of 8 and 12 DHW indicate high and very high risks of bleaching-related coral mortality, respectively. Our model indicated that regions experiencing at least 5 DHW during 2013-2016 were more suitable for A. amadelpha invasion. Cox et al. (2017) speculate that recent warming events in 2015 caused a die-back in abundant native intertidal macroalgae, facilitating the rapid spread of leather mudweed in this intertidal environment. They observed an opportunistic increase in abundance of A. amadelpha following the reduction in brown algae (Padina sanctae-crucis and Dictyota spp.), which dominate multiple intertidal sites around Oahu (Cox et al., 2017). This invasion notably followed a period of unseasonably warm, calm seas and a series of severe bleaching events from 2014 to 2017. The impact of warmer ocean temperatures and the resulting die-back of native macroalgae may also increase the success of more thermally resistant invasive species. We hypothesize that reefs destabilized from repeated heat stress are more susceptible to incursion of invasive species like A. amadelpha. More ecophysiological research is needed on the thermal tolerance and growth response of A. amadelpha.
In addition to using our model to highlight important stressors and the most vulnerable coastal regions island-wide based on present-day climatological conditions, we used our framework to forecast how reef vulnerability might change following a 25% increase in statewide maximum DHW. With all other predictors held constant (i.e., held at present-day values), our model predicted changes in habitat suitability per 250 m pixel ranging from decreases of up to 25% to increases as high as 61%. The west coast of Oahu and much of the southern coast experienced little to no increases, or even slight decreases, in overall vulnerability to A. amadelpha invasion. Portions of the northern and eastern coasts and patches of the Honolulu metro area displayed the most dramatic increases in vulnerability around Oahu (Figure 7, inset A). Across the other MHI, we observed increased vulnerability around the southern and western shores of Kauai and Maui Nui, and around northwestern and eastern Hawaii, including Hilo Bay (Figure 7).

Implications for Management
Our analyses did not reveal any obvious differences in A. amadelpha colonization patterns along a depth gradient. We caution that the absence of a clear difference between the shallow and mesophotic populations of the alga do not guarantee the effectiveness of a uniform management approach when dealing with the removal of this invader. We note that A. amadelpha is expensive, time-consuming, and-in the case of the mesophotic populations-potentially impossible to entirely eliminate for a given region. The aim of this study is specifically to guide management efforts during ongoing monitoring by identifying the portions of the coastline around the state that are most vulnerable to invasion in the near-term based on certain environmental characteristics.
We suggest that managers used the maps we have generated here to prioritize regions of concern. These maps may serve as early warnings for regions that provide recreation for residents and tourists, serve as important sustenance fisheries, or maintain inherent value due to their cultural significance. State agencies that monitor, protect, and preserve valuable natural resources across the land and sea, particularly in the face of a changing climate, must also contend with limited manpower and financial support. We hope that these maps may simplify and expedite the decision-making process of managers regarding monitoring and testing efforts to mitigate damage due to the incursion of invasive species.

CONCLUSION
Our model pinpointed one-quarter of the coastline of the MHI as particularly susceptible to A. amadelpha invasion (Figure 5). Our model indicates that A. amadelpha is a "Goldilocks"-style opportunist that favors regions exposed to a moderate amount of stress, particularly in the form of fishing pressure, coastal development, and warming seas.

DATA AVAILABILITY
The cleaned data sheet and R code used for these analyses are included in the Supplementary Material. Requests for additional raw data files and products may be submitted to the corresponding author.

AUTHOR CONTRIBUTIONS
LV: study design. RW and HS: data sources. LV and OW: coding. LV, RT, RW, and HS: analysis. LV: wrote the first draft of the manuscript. All authors wrote subsections, contributed edits to revisions, and approved the submitted version of the manuscript.