Factors Influencing the Faunal Recolonization of Restored Thornscrub Forest Habitats ^†

Abstract

Tamaulipan thornscrub forests (thornforests) have high ecological and economic values, yet over 90% of these forests have been lost, and they remain threatened, making them a conservation hotspot. For decades, federal, state, NGO, and corporate entities have been acquiring land and actively or passively restoring these forests, but results have been mixed and seldom monitored. This study characterized and quantified faunal communities of restored thornforest habitats in south Texas and examined the relationships between restored faunal communities and key site characteristics and environmental factors. We surveyed and analyzed mammal, bird, Lepidoptera, and herptile communities within 12 restored sites in the Lower Rio Grande Valley of southernmost Texas, USA. The site and environmental factors that influenced animal community composition, richness, diversity, and abundance varied widely among taxa. Time since restoration began and method of restoration influenced many community metrics, whereas patch size and extent of isolation influenced few. Several aspects of restored plant communities were influential, especially ground layer diversity, and high invasive plant cover negatively impacted many animal community metrics. If actively restoring a site, efforts to control invasive plants, foster native plant diversity, and provide a nearby water source are likely the most effective ways to promote faunal recolonization.

1. Introduction

Tamaulipan thornscrub forests (thornforests) are ecologically and economically vital habitats in the Lower Rio Grande Valley (LRGV) of southernmost Texas, USA, and northern Tamaulipas, Mexico. The LRGV is a fertile deltaic plain that slopes away from the Rio Grande River, with a semi-arid and subtropical climate; annual rainfall ranges from 38 to 76 cm, mean daily highs in July range from 35.0 to 36.7 °C, and mean daily lows in January range from 8.9 to 10.6 °C [1]. Vegetation is typically drought resistant and mostly characterized by spiny shrubs, stunted trees, and dense brush (Figure 1) [1]. Thornforests provide forage and habitat for many native reptiles and mammals, including endangered ocelots (Leopardus pardalis), and support a high diversity and abundance of beetles, bees, and resident and migratory birds and butterflies, with the latter driving a hundred-million-dollar per year regional ecotourism industry [2,3].

The LRGV has a long history of human modification, especially since the arrival of Europeans in the 1700s [4]. Key modifications include overgrazing by sheep in the 18th and 19th centuries [1], hydrologic changes to the Rio Grande River in the 20th century [1], habitat clearing and expansion of row cropping in the 1950s and 1960s [4], the introduction of several invasive Old World grasses in the late 19th and 20th centuries [1,4], and most recently, land use changes associated with increased urbanization and industrialization arising from rapid regional population and economic growth. These changes have resulted in the loss of over 90% of historic thornforests and significantly fragmentated relic habitats [1]. The high biodiversity and ecological significance of thornforests, combined with their extensive historic habitat loss and ongoing threats from human development, make Tamaulipan thornforests a global conservation hotspot and have made their conservation and restoration in the LRGV a national and regional priority [5].

Since the 1980s, federal, state, NGO, and corporate entities have been working to restore Tamaulipan thornscrub habitats, especially thornforests, in south Texas by acquiring properties that had been converted to farmland or rangeland and taking measures to return their ecological structure and function to reference conditions through either passive or active restoration. Following acquisition by the United States Fish and Wildlife Service (USFWS) or the Texas Parks and Wildlife Department (TPWD), land was either simply protected, allowing for unassisted regeneration (passive restoration), or it was actively planted with native seeds or seedlings (active restoration). To date, more than 16,000 acres have been actively restored and even more land has been protected and passively restored [6]. In nearly all cases, however, a focus on maximizing the number of acres restored has meant that there have been insufficient resources for monitoring. This has made it challenging to evaluate the success of restoration efforts on Tamaulipan thornscrub habitats or to identify which methods of restoration have been the most effective. Recently, land managers have been reconsidering restoration strategies; unfortunately, many important questions remain unanswered, especially regarding differences between actively and passively restored sites and how patch size, degree of isolation, and edge effects influence restoration outcomes.

Very few studies have evaluated Tamaulipan thornscrub restoration outcomes, and they have generally focused on woody plant communities. Ewing and Best [4] evaluated and compared woody species performance in recently planted restoration sites and in a remnant community. Alexander et al. [7] evaluate how well certain restoration methods helped transplanted woody seedlings to overcome the common stressors of competition, herbivory, and drought. Alanís-Rodríguez et al. [8] compared woody plant communities in reference areas and areas with assisted and unassisted ecological succession. Perez et al. [9] measured woody species regeneration in abandoned plots with different land-use histories. Other studies simply examined Tamaulipan thornscrub plant communities after disturbance without reference comparisons, specifically after livestock were removed [10] and after fires [11].

This research has informed thornscrub planting strategies, but none of these studies investigated the faunal recolonization of areas undergoing restoration. Animals are key components of these ecosystems and core to the mission of managing agencies like the TPWD and the USFWS. To the best of our knowledge, no multi-taxa studies have been conducted to quantify and characterize the wildlife communities present in restored Tamaulipan thornscrub forests.

Surveying multiple taxa of wildlife is important. Studies show that surveying only indicator species results in inaccurate pictures of wildlife communities because different faunal groups return at different rates and select sites based on different factors [12]. Conversely, monitoring approaches that consider multiple species or taxa provide a more robust characterization of animals communities [13]. Including multiple taxa may be particularly important where habitats have been modified because of the greater potential for variation in animal responses to such changes [14]. Understanding these variable responses is particularly important in our study region. For these reasons, we surveyed the community structure of four key animal taxa, specifically mammals, birds, Lepidoptera, and herpetofauna, and quantified the richness, abundance, and diversity of these focal taxa at 12 different thornforest restoration sites within the LRGV.

To inform restoration efforts going forward, beyond assessing restoration outcomes, we must understand which factors influence these outcomes and how. Grman et al. [15] posit that four classes of drivers affect restored communities: management decisions, site characteristics, landscape context, and historical factors. For example, management actions, like planting seedlings, can be overwhelmed by site characteristics, like soil conditions [15,16]; soil conditions can be affected by historical factors [17]; and historical factors, like prior invasive plant dominance, can impact site characteristics and management decisions [18].

It is well established that site characteristics and environmental factors impact the recolonization of wildlife. One key factor is time. A meta-analysis of 83 terrestrial restoration studies found that mean biodiversity increased as restorations aged [19]. However, the relationship between time and diversity is more complex. Recolonization time varies among faunal groups; for example, one study found generalist foraging mammals recolonized much more rapidly than reptiles [12]. Similarly, a long-term study in post-mining forests found that the strength of ‘filters’ impeding wildlife recolonization and reducing population persistence varied over time and differed among taxa [20]. Also key to animal community composition and core to the ‘field of dreams’ restoration hypothesis are characteristics of a site’s plant community, such as the diversity or structure of different forest layers [21,22,23,24]. Environmental factors, like moisture and temperature, clearly impact animal communities both directly and indirectly, e.g., by influencing plant communities [25].

Other site characteristics influence restored animal communities but have less straightforward impacts. Invasive species have had many documented effects on restoration outcomes, especially animal species richness and diversity, but these effects are strongly context-dependent [26]. Alternatively, characteristics like patch size and shape may have strongly taxon-specific effects. A study in human-impacted habitats in the LRGV found lower plant and Lepidoptera richness in larger patches but higher plant and Lepidoptera richness in patches with higher edge-to-interior ratios [27]. Conversely, some butterfly species were more likely to colonize larger patches in Surrey, UK [28], and most bird species were positively but weakly associated with larger restored patches in an agriculturally fragmented landscape in Illinois, USA [29]. However, patch size had no effect on bird species richness in patches within logged areas in Victoria, Australia [30].

Regarding landscape context, many studies demonstrate how reduced connectivity, which we approximate as spatial isolation in this study, can impact the composition and abundance of wildlife. Isolation can restrict the movement and dispersal of wildlife, creating sinks around reserves [31] and increasing genetic differentiation between neighboring populations [32,33]. Also important for animals is the distance to essential resources like water, both permanent and ephemeral. Riparian zones often support disproportionately high animal species richness and abundance [34], and ephemeral wetlands can contribute to landscape connectivity [35] and enhance vertebrate activity and diversity in certain landscapes [36]. It logically follows that landscape effects on a given taxon depend on its movement or dispersal ability and, thus, would vary among taxa. However, some research found that the degree of isolation was a poor predictor of patch occupancy for most species and that the properties of the intervening matrix were more important [37].

The method of restoration can also affect faunal recolonization, though research varies widely on this topic. For example, Russell et al. [38] found that the diversity and abundance of wildlife were greater in actively restored dunes, whereas Trujillo-Miranda et al. [39] found equivalent richness and diversity in actively and passively restored cloud forests. A meta-analysis of 150 studies on the effectiveness of active versus passive restoration of forests found that simply ending human land use, regardless of prior modifications, was sufficient for most forest communities to recover and that actively planting trees in former agricultural land (as in our study) did not result in consistently faster or more complete recovery than passive restoration [40]. We know thornforest plant communities have differed significantly between actively and passively restored sites [7,8,41], but whether this translates into differences in thornforest animal communities has not been studied.

For thornforests in the LRGV, given the limits imposed by historical realities and data availability, we considered three of the four classes of drivers affecting restoration identified by Grman et al. [15]. Replanting and invasive grass suppression actions have varied to some degree, but the core management decision has been whether to pursue active or passive restoration [4]. Study site characteristics were directly measured or mined from existing datasets in a parallel study [41], and landscape context was evaluated using GIS analyses of public data. Historical factors, however, cannot be considered beyond our general understanding of prior human modifications in the LRGV, as described above. Records could not be found and probably do not exist regarding the specific land use histories of our 12 study sites. We can only say with confidence that all 12 sites were clear-cut at some point and either grazed (by cattle and/or sheep), plowed and farmed for row crops, or both grazed and farmed, but we do not know when or for how long, except regarding when the land was protected and restoration began.

Thus, the objectives of this study are to (1) characterize and quantify the faunal communities of four taxa (mammals, birds, Lepidoptera, and herptiles) within 12 restored thornscrub forest habitats in south Texas and (2) examine the relationships between faunal communities and key site characteristics, such as restoration method (passive versus active), habitat patch size, time since restoration began, edge-to-interior ratio, degree of isolation, characteristics of the local plant community, and abiotic factors, like soil temperature and soil moisture content. By intentionally selecting restored sites that vary in these factors, we aimed to identify important relationships and evaluate how these factors influenced the recolonization of different animals after restoration. Our study could help inform future restoration and management efforts by informing decisions about land acquisition or restoration methodologies and by providing specific recommendations for promoting faunal return to restored forest sites.

2. Materials and Methods

2.1. Site Selection

This study was performed at 12 field sites located in Cameron County, Texas (Figure 2). Sites were selected to represent variation in the time since restoration began, restoration method, patch size, and degree of isolation. Given the ubiquity of human impacts and uncertainty regarding land use history in the LRGV, as discussed above, land use history was not a selection factor, and no clearly “pristine” sites were available within the focal region that could be used as traditional reference sites [27,42]. However, our Goat Island site is notable; its land use history includes at least 77 years of protection from most human impacts—though it has been hunted regularly in that time—and it may have experienced minimal agricultural use prior to 1946. All the other study sites were previously used heavily for agriculture but were protected for restoration 15–70 years prior to being surveyed (

Table 1. Summary of focal characteristics for our 12 study sites, ordered by restoration time.

Site Name	Restoration Method	Estimated Year Restoration Began	Degree of Isolation	Patch Size (ha)
Goat Island (GI)	Passive	1945	53%	104
Longoria Unit (LU)	Passive	1955	85%	122
Tucker Unit (TU)	Passive	1955	59%	72.7
Anacua Unit (AU)	Passive	1985	81%	56.4
Arroyo Colorado Unit (AC)	Passive	1985	78%	175
Ebony Unit (EU)	Passive	1985	29%	89.0
Duck Head (DH)	Active	1990	88%	26.8
Phillips Banco (PB)	Active	1993	61%	28.2
Garza-Cavazos (GC)	Active	1994	48%	4.51
Fish Hatchery (FH)	Active	2002	92%	18.4
Villa Nueva (VN)	Active	2004	72%	18.2
Tahuacal Banco (TB)	Active	2007	49%	5.93

). The nature and extent of each site’s alteration due to land use history is unknown, which is unfortunate because historical factors are important to faunal recolonization [15,40,43]. Six sites are owned and managed by the Texas Parks and Wildlife Department (TPWD), five by the United States Fish and Wildlife Service (USFWS), and one by the University of Texas Rio Grande Valley. The exact dates when each restoration began were not always known, but sufficiently accurate estimates were provided by land managers. Most sites used exclusively either an active or a passive restoration approach; however, two had very small planted sections but otherwise regenerated naturally, so we grouped sites based on the primary restoration method (active or passive). Isolation was calculated as the percentage of land cover within 1 km from the center of the patch that was not thornscrub habitat, excluding the patch itself (see also Supplementary S2).

We explicitly selected sites from a pool of 35 suitable candidates to represent independent variation in restoration time, method, size, and isolation (see Supplementary S2). However, variation in these traits could not be fully statistically independent given the region’s history and the choices available. Specifically, the largest patches were generally acquired and protected earlier and restored passively, whereas patches acquired later tended to be smaller and actively restored. Thus, some site characteristics were unavoidably confounded because of our limited set of options, despite explicit efforts to avoid this.

2.2. Medium and Large Mammal Surveys

To survey mammals, we established three camera traps per study site (n = 36), with each located at one of the avian sampling points. Cameras were placed ca. 1 m above the ground and at least 200 m apart. If possible, cameras were positioned along or near a game trail or clearing to enhance visibility. To ensure event independence, photographs of the same species captured within 60 min by the same camera were treated as a single occurrence unless the number of individuals increased [44]. For each camera event, the date, time of day, species, and location were recorded. We checked cameras and replaced SD cards every 2–4 weeks for a period of 3–4 months (April–August). We analyzed mammal photographs using the CameraSweet set of programs [45] (Small Wild Cat Conservation Foundation, Corrales, NM, USA). This involved manually sorting photographs by mammal species and abundances and then analyzing photograph metadata to quantify the number of camera trapping days and independent captures per species for each camera. Trapping days ranged from 28 to 118 days, so we normalized mammal abundance values based on the sampling effort. See Supplementary S2 for additional details on our sampling protocols, including justifications for particular methods used for all the focal taxa.

2.3. Avian Surveys

We conducted audiovisual point-count bird surveys at the same three locations per site (n = 36) where we sampled mammals; points were ≥200 m apart to ensure independence [46]. We performed two surveys at each sampling point, one in May and one in June 2022 (n = 72). This reduced our likelihood of missing species [47] and coincided with the breeding season, meaning more birds would be vocalizing and fewer migrants would be present. To optimize detection, all surveys were conducted between sunrise and 10:30 AM (to avoid high temperatures) and only when it was not raining and the windspeed was below 30 km/h. We varied the order in which survey points were visited between surveys to minimize the effects of time or temperature. Following protocols in Ralph et al. [48], the surveys used an unlimited radius but noted detections within 50 m and lasted 10 min, with detections partitioned into intervals of 0–3, 3–5, and 5–10 min, which allows for direct comparisons with surveys using other time intervals [49].

2.4. Lepidoptera Surveys

We sampled Lepidoptera via bait traps twice at each of the same 36 sampling points, once in May and once in July 2022, but combined both samples from each point for analyses (n = 36). Traps were ca. 60 × 60 × 90 cm and made of fine netting, with a bait bowl suspended 4 cm below an opening at the bottom [50], which was filled with overripe fermenting bananas. We hung baited traps from trees 1–1.5 m above the ground in the morning and left each deployed for 24 h. Captured individuals were carefully removed, photographed, and immediately released. Traps were not deployed if rain was forecast within the subsequent 24 h.

2.5. Herpetofauna Surveys

Reptile and amphibian surveys used artificial cover objects (ACOs) and time-constrained area searches. We installed three ACOs at each of the three sampling points within each study site (n = 108). ACOs were made of 9.5 mm thick untreated and unpainted plywood; two were 60 × 60 cm, and one was 120 × 60 cm at each point. Installation involved placing ACOs horizontally under various degrees of woody canopy cover within 10 m of the sampling point to increase the likelihood of detecting more species. We then left the ACOs undisturbed for at least 3 weeks before beginning data collection because herptiles tend not to immediately utilize new cover objects [51]. From April to August 2022, each ACO was checked a total of five times. Checks occurred every 2–4 weeks. ACOs were slowly lifted, each species present underneath was photographed, and boards were slowly set back down to avoid harming any animals.

After two months of seeing herptiles in the vicinity but observing very few under the ACOs, we additionally performed two time-constrained area searches at each of the 36 established sampling points, with one in July and one in August (n = 72) [52]. For each search, two surveyors searched the area within ca. 50 m of a sampling point for 15 min; debris or rocks that could provide cover for herpetofauna were briefly removed, leaf litter was searched, and trees were scanned. Observed herptiles were identified immediately and recorded or photographed; search time was not used to capture or identify species. To minimize extraneous variation, each set of surveys was conducted within a specific timeframe and under similar weather conditions within a one-month period [53]. As for birds, we also varied the order in which we visited sites to minimize variation related to time. See Supplementary S2 for additional details on the sampling methods used for all the taxa.

2.6. Vegetation Surveys and Environmental Variables

At each of the 36 sampling points, the vegetation was systematically surveyed in two sampling areas (6 areas per site; n = 72) using different sampling approaches tailored to different forest layers that, together, quantified forest structure and plant community composition. These methods allowed us to quantify large canopy layer trees; understory trees, shrubs, and climbing species; and all ground layer vegetation, including low-statured grasses and forbs. Data on soil temperature and moisture content were also collected in situ.

These vegetation surveys were performed as part of a parallel research project, so the sampling and analytical methods utilized are not described herein but may be found in Garrett [41]. However, we do use some of the data and results from the plant community surveys in Garrett [41] in our analyses and to inform our findings.

To better explain the observed variance not attributable to our focal site characteristics of the time since restoration began, patch size, edge-to-interior ratio, and restoration method, we additionally considered various environmental and geographic metrics that we either measured in situ at our sampling points or mined from publicly available GIS data for each sampling point or site by Garrett [41].

2.7. Analysis

First, we normalized any response variables that differed in sampling effort based on the actual effort exerted (e.g., operational days for trail cameras). Next, to explore relationships and potential correlations among our focal site and environmental variables, we conducted a principal component analysis (PCA) in R (version 4.3) (R Foundation for Statistical Computing, Vienna, Austria) using the PCA() function from the ‘FactoMineR’ package.

To characterize the observed faunal communities and explore the relationships among species abundances, site characteristics, and environmental variables, we performed separate multivariate analyses for the mammal, bird, Lepidoptera, herpetofauna, and ensemble (all the groups combined) communities. We used the metaMDS() function in the ‘vegan’ package in R to fit nonmetric multidimensional scaling (NMS) ordinations using Bray–Curtis dissimilarity values. NMS is ideal for comparing community compositions because it can powerfully quantify and visualize differences among samples, whereas other multivariate approaches, such as redundancy analysis (RDA) or canonical correspondence analysis (CCA), are better suited for simplifying or combining multivariate datasets. Sampling points were the observational units (n = 36) for most ordinations, with abundance values from individual surveys summed for each sampling point. For birds and herptiles, we combined the three sampling points at each site and used sites as the observational unit (n = 12) to reduce the NMS fit stress to acceptable levels (<0.2). See Supplementary S2 for additional details on the data used in each ordination.

We then used the envfit() function from the ‘vegan’ package in R to fit relevant site and environmental variables to our NMS ordinations, including the method of restoration (categorical) and continuous variables for patch size, time since restoration, degree of isolation, interior-to-edge ratio, soil moisture content, soil temperature, distance to permanent water, distance to temporary water (e.g., seasonal wetlands), invasive grass cover, canopy cover, canopy tree density, understory cover, total plant richness, total plant diversity, ground layer plant richness, ground layer plant diversity, understory layer plant richness, understory layer plant diversity, canopy plant richness, canopy plant diversity, natural log of the ratio of native-to-exotic plant covers, feral hog disturbance, and combined total exotic plant cover. For a full list of values and additional details regarding these environmental factors beyond the focal site characteristics, see Garrett [41]. Hog disturbance was plotted to illustrate patterns in the observed animal communities but was not included in subsequent PerMANCOVAs (see below). We visualized our ordinations and environmental fits in R using graphing functions in the ‘ggplot2’ package. To avoid overplotting, we plotted vectors only for the most influential species and the environmental factors most strongly associated with the spread of sites. We calculated and displayed 95% confidence ellipses around the centroids (hypothetical average community composition) for the groups defined by the method of restoration (active vs. passive).

To examine the effects of our categorical and continuous environmental variables on community composition, we used the adonis2() function from the ‘vegan’ package in R to perform a permutational multiple analysis of covariance (PerMANCOVA) for each community. Given the abundance of factors considered, careful model selection was critical. Supplementary S2 provides a detailed description of our model selection process. Essentially, we first removed the weaker predictors among any correlated independent variables and then used a stepwise model selection function in R that considered Akaike information criterion (AIC) values of alternative models with human oversight.

To further examine how site characteristics related to restoration outcomes, we then performed a series of univariate analyses for the same five taxonomic groups. For each, we fit linear or permutational linear models using the lm() or lmp() functions in R and performed ANCOVAs or multiple regressions using marginal (Type III) sums of squares to examine the effects of our focal environmental variables on three key community-level response variables: richness, abundance, and diversity. Our full models included the same variables considered in our PerMANCOVAs, so we similarly purged the least explanatory of any correlated environmental variables and then used the step() function in R with forward and backward model selection based on AIC values to prune each model.

To confirm that each PerMANCOVA, multiple regression, and ANCOVA met all the relevant model assumptions, we performed Shapiro–Wilk tests of normality on model residuals and Breusch–Pagan tests of homoscedasticity for all the linear models, and we calculated the variance inflation factor of all the model terms for all the models to quantify multicollinearity using the vif() function in R. A p-value of 0.05 was used to determine significance and test model assumptions.

3. Results

3.1. Principal Component Analysis of Environmental Factors

Principal component analysis showed that most variation among the sites was related to plant diversity (PC1, which explained 35.5% of the variance), focal site traits (PC2, 14.5%), and forest canopy layer density (PC3, 10.1%). The main contributors to PC2 were patch size (20.5%), interior-to-edge ratio (19%), and time since restoration (17.5%). Figure S1 provides details on the other principal components and the variables contributing to the first three.

Figure S2a plots the contributions of all the environmental factors to the first two principal components; clustered vectors are typically correlated to some degree. For example, patch size and time since restoration are clustered and correlated (larger sites tend to be older). Figure S2b plots the positions of sites using the same axes, so proximity denotes the similarity of environmental conditions, and groups sites based on the restoration method. Figure S3 combines Figure S2a,b to more clearly illustrate associations between particular sites and different ranges of values for environmental factors.

3.2. Mammals

3.2.1. Multivariate Analysis of Mammal Communities

In total, we observed 18 mammal species (

Table 2. Total independent observations per 30 days of sampling for each photographed mammal species at the 12 study sites. Values were normalized to account for differences in sampling effort. These site-level values are the sums of the values quantified for the three sampling points within each site; supplementary Table S1 lists mammal observations by species for all 36 sampling points. The abbreviations shown are the symbols used for the species in the ordinations. Legend: * introduced; ^† invasive; AU, Anacua Unit; AC, Arroyo Colorado Unit; DH, Duck Head; EU, Ebony Unit; FH, Fish Hatchery; GC, Garza-Cavazos; GI, Goat Island; LU, Longoria Unit; PB, Phillips Banco; TB, Tahuacal Banco; TU, Tucker Unit; VN, Villa Nueva.

Scientific Name	Common Name	Abbreviation	AU	AC	DH	EU	FH	GC	GI	LU	PB	TB	TU	VN	Species Total
Bos taurus	Domestic Cattle *	Cow	0	0	0	0	0	0	0	0	0	0	0.9	0	0.9
Boselaphus tragocamelus	Nilgai *	Nilg	0	6.6	0	0	0.3	0	0	0	0	0	0	1.6	8.4
Canis latrans	Coyote	Coyo	0.6	2.9	0.8	1.7	5.1	6.5	1.1	1.5	0.3	1.3	0.5	0.3	22.4
Canis lupus familiaris	Domestic Dog *	Dog	0	0	0.4	0	0.3	0	0	1.2	0	0	0.3	0	2.2
Chiroptera spp.	Bat species	Bat	0	0	0	0	0	0	0	0	0	0.7	0	0.3	1.0
Dasypus novemcinctus	Nine-banded Armadillo	Arma	0	0	27.5	0.6	9.5	2.1	12.8	2.3	1.4	0	0	0	56.1
Dicotyles tajacu	Collared Peccary	Jave	2.6	1.5	4.5	2.8	0.3	5.1	4.1	11.6	0.6	1.9	0.6	0.8	36.5
Didelphis virginiana	Virginia Opossum	Opos	0	0	7.8	0	0	0	1.2	0.3	0	0	0	0	9.3
Felis catus	Domestic Cat *†	Cat	0	0	0	0	0.3	0	0	0	0	0	0	0	0.3
Homo sapiens	Human	Huma	0	0	1.1	0	0	1.6	0	0.9	0.5	1.3	1.8	0.3	7.5
Leopardus pardalis	Northern Ocelot	Ocel	0	0	0	0	0	0	4	0	0	0	0	0	4.0
Lynx rufus	Bobcat	Bobc	0	2.0	4.1	0.9	0	1.2	1.5	11.6	0	0	2.0	0	23.3
Mephitis mephitis	Striped Skunk	Skun	0	0	2.9	0	0.7	0	0	0	0	0	0	0	3.5
Odocoileus virginianus	White-tailed Deer	Deer	0	40.8	0	0	0	1.8	19.4	0	0	0.2	0	0	62.2
Procyon lotor	Common Raccoon	Racc	0	0.9	31.1	0.6	9.3	32.9	10.3	2.9	2.0	0	1.8	0.8	92.6
Rattus spp.	Rat species	Rat	0	0	0	0	0	0.5	0	2.3	0	2.1	0	0	4.9
Sus scrofa	Feral Hog *†	Hog	0	0	0	0	4	0	0	0	0	0	0	0	4.0
Sylvilagus floridanus	Cottontail Rabbit	Rabb	0	0	3.3	1.8	3.9	0	3.6	1.2	0.7	0	0	0	14.5
Site Total			3.2	54.6	83.7	8.4	33.7	51.6	57.9	35.8	5.5	7.6	7.7	4.1	353.7

), including 12 thornforest species and two groups identifiable only as bats or rats. Trail cameras do not reliably detect bats or rats and seldom provide adequate detail to identify either to the species level. For comparison, the species list for the nearby Laguna Atascosa National Wildlife Refuge includes 20 thornforest mammal species excluding bats and rats and 32 including them [54]. The five most common species were Procyon lotor (Common Raccoon, 92.6 observations per 30 sampling days; OPM henceforth), Odocoileus virginianus (White-tailed Deer, 62.2 OPM), Dasypus novemcinctus (Nine-banded Armadillo, 56.1 OPM), Dicotyles tajacu (Collared Peccary also known as Javelina, 36.5 OPM), and Lynx rufus (Bobcat, 23.3 OPM).

Table 2 lists the number of observations for each mammal species at our 12 study sites; Table S1 lists observations by species for all 36 sampling points. There was considerable variation in the observations among both species and sites. Leopardus pardalis (Ocelot; a US- and state-listed endangered species) was observed but only at Goat Island.

NMS ordinations of the observed mammal communities are shown in Figure 3 and Figure S4. All our NMS ordinations illustrate community compositions and similarities among communities, which are represented as the positions and spatial proximities of points (or labels in supplemental ordinations), respectively. Distance to permanent water was positively associated with abundances of peccaries and bobcats. Invasive grass cover and combined exotic plant cover were strongly correlated (because invasive grasses were the predominant type of exotic species observed) and were associated with greater human abundance. Total plant diversity, total plant richness, and ground layer plant diversity were associated with higher abundances of rabbits, bobcats, and feral hogs. Nilgai abundance was positively associated with the degree of isolation of restored thornforest patches.

Table 3. Results of multivariate (PerMANCOVA) and univariate (multiple regression using Type III sums of squares) analyses examining the effects of site traits on mammal community composition, species richness, diversity (Shannon index, H’), and total abundance. The effect size is the predicted change in a univariate response variable given a one unit increase in the predictor variable. More complex models were considered for each response variable prior to model pruning; factors not included herein were pruned either due to insignificance or to avoid multicollinearity. Full results tables for individual models can be found in the supplemental tables shown. Legend: d.f., degrees of freedom; ***, p < 0.001; **, 0.001 ≤ p < 0.01; *, 0.01 ≤ p < 0.05; ., 0.05 ≤ p < 0.1.

Response	Predictor	Test	Effect Size	d.f.	F	p		Results Table
Composition	Patch size	PerMANCOVA	-	1	3.54	0.0001	***	S2
Composition	Method of restoration	PerMANCOVA	-	1	3.15	0.0001	***	S2
Composition	Time since restoration	PerMANCOVA	-	1	2.46	0.0060	**	S2
Composition	Distance to permanent water	PerMANCOVA	-	1	2.30	0.0090	**	S2
Composition	Ground layer plant diversity	PerMANCOVA	-	1	1.60	0.0939	.	S2
Composition	Distance to temporary water	PerMANCOVA	-	1	1.80	0.0400	*	S2
Richness	Ground layer plant richness	Multiple regression	+0.372	1	12.27	0.0016	**	S3
Richness	Soil temperature	Multiple regression	+0.442	1	7.32	0.0115	*	S3
Richness	Isolation	Multiple regression	+0.044	1	5.54	0.0258	*	S3
Richness	Understory plant richness	Multiple regression	−1.809	1	16.81	0.0003	***	S3
Richness	Understory plant diversity	Multiple regression	+4.074	1	6.60	0.0158	*	S3
Diversity	Understory total cover	Multiple regression	−0.014	1	7.18	0.0118	*	S4
Diversity	Ground layer plant richness	Multiple regression	+0.052	1	4.58	0.0406	*	S4
Diversity	Soil temperature	Multiple regression	+0.069	1	2.88	0.1003		S4
Abundance	Ground layer plant diversity	Multiple regression	+17.420	1	4.79	0.0363	*	S5
Abundance	Total plant richness	Multiple regression	−1.630	1	2.00	0.1169		S5

and Table S2 show the PerMANCOVA results examining the effects of site traits on the mammal community compositions shown in the NMS ordination. Habitat patch size, method of restoration, time since restoration began, distance to permanent water, and distance to temporary water had significant effects on mammal community composition, whereas ground layer plant diversity had a marginal effect on composition. For all our models, most pruned (removed) factors lacked a significant relationship with the response variable, but some were removed because they were correlated with other model terms; it is possible that some of the latter had a significant effect on the response variables.

3.2.2. Mammal Community Univariate Analyses

Mammal richness averaged 3.82 ± 2.49 species across the sampling points and was influenced by ground layer plant richness, soil temperature, extent of isolation, understory plant richness, and understory plant diversity (Table 3 and Table S3; Figure 4a and Figures S5–S8). Mammal diversity (Shannon–Weiner index, H’) averaged 0.91 ± 0.59 and was influenced by understory total plant cover and ground layer plant richness (Table 3 and Table S4; Figure 4b and Figure S9). Mammal abundance averaged 10.40 ± 16.92 OPM and was associated with ground layer plant diversity only (Table 3 and Table S5; Figure S10).

3.3. Birds

3.3.1. Multivariate Analyses of Bird Communities

We observed a total of 53 bird species across all the surveys and sampling points (

Table 4. Total number of surveys (out of six performed) in which each bird species was observed at the 12 study sites. These site-level values are the sums of the values quantified for the three sampling points within each site; supplementary Table S2 lists the number of surveys (out of two) in which each bird species was observed at all 36 sampling points. The abbreviations shown are the symbols used for the species in the ordinations. Legend: *, introduced; ^†, invasive; AU, Anacua Unit; AC, Arroyo Colorado Unit; DH, Duck Head; EU, Ebony Unit; FH, Fish Hatchery; GC, Garza-Cavazos; GI, Goat Island; LU, Longoria Unit; PB, Phillips Banco; TB, Tahuacal Banco; TU, Tucker Unit; VN, Villa Nueva.

Scientific Name	Common Name	Abbreviation	AU	AC	DH	EU	FH	GC	GI	LU	PB	TB	TU	VN	Species Total
Agelaius phoeniceus	Red-winged Blackbird	Agepho	4	1	0	0	0	1	1	0	0	4	2	0	13
Amazona viridigenalis	Red-crowned Parrot	Amavir	0	0	4	0	0	0	0	0	0	0	0	0	4
Amphispiza bilineata	Black-throated Sparrow	Ampbil	0	0	0	1	0	0	0	0	0	0	0	0	1
Archilochus sp.	Hummingbird sp.	Arcsp.	0	0	0	0	0	0	0	0	0	0	0	1	1
Ardea alba	Great Egret	Ardalb	0	1	0	0	0	0	0	0	0	0	0	1	2
Arremonops rufivirgatus	Olive Sparrow	Arrruf	6	5	5	5	5	3	6	5	4	6	6	3	59
Baeolophus atricristatus	Black-crested Titmouse	Beoatr	3	0	0	0	1	0	2	0	4	0	3	3	16
Buteo lineatus	Red-Shouldered Hawk	Butlin	0	1	0	0	0	0	0	0	0	1	1	0	3
Cardinalis cardinalis	Northern Cardinal	Carcar	1	5	3	1	1	3	3	4	1	5	0	2	29
Cathartes aura	Turkey Vulture	Cataur	0	1	0	0	1	0	1	0	0	0	0	0	3
Charadrius vociferus	Killdeer	Chavoc	0	0	0	1	0	0	0	0	0	0	0	0	1
Chordeiles minor	Common Nighthawk	Chomin	1	0	0	0	0	0	1	0	0	1	0	0	3
Coccyzus americanus	Yellow-billed Cuckoo	Cocame	3	4	0	5	0	3	2	6	4	1	3	4	35
Columbina passerina	Common Ground Dove	Colpas	1	2	0	3	1	1	1	0	0	5	0	0	14
Colinus virginianus	Northern Bobwhite	Colvir	1	0	0	0	0	0	0	0	0	0	0	0	1
Crotophaga sulcirostris	Groove-billed Ani	Crosul	4	0	0	1	3	1	0	0	0	1	1	2	13
Cyanocorax yncas	Green Jay	Cyaync	0	2	1	0	3	1	0	3	2	2	1	4	19
Dendrocygna autumnalis	Black-Bellied Whistling Duck	Denaut	1	0	2	1	1	0	2	2	0	0	2	1	12
Dryobates scalaris	Ladder-backed Woodpecker	Drysca	2	3	0	0	1	0	2	4	3	2	2	3	22
Egretta thula	Snowy Egret	Egrthu	0	0	0	1	0	0	0	0	0	0	0	0	1
Empidonax sp.	Empidonax species	Empsp.	1	1	0	0	0	1	0	0	1	1	0	0	5
Geothlypis trichas	Common Yellowthroat	Geotri	0	0	1	0	0	0	0	0	0	0	0	0	1
Icterus galbula	Baltimore Oriole	Ictgal	0	0	0	0	0	0	0	0	1	0	0	0	1
Icterus gularis	Altamira Oriole	Ictgul	0	1	0	0	1	2	0	0	0	0	0	0	4
Lanius ludovicianus	Loggerhead Shrike	Lanlud	1	0	0	0	0	0	0	0	0	0	0	0	1
Leptotila verreauxi	White-tipped Dove	Lepver	6	0	0	2	2	0	2	4	2	0	2	2	22
Leucophaeus atricilla	Laughing Gull	Leuatr	0	4	0	0	0	0	5	0	1	0	1	0	11
Melanerpes aurifrons	Golden-Fronted Woodpecker	Melaur	6	2	6	4	5	2	2	6	4	6	6	5	54
Mimus polyglottos	Northern Mockingbird	Mimpol	5	0	1	2	0	1	0	0	0	1	0	0	10
Molothrus aeneus	Bronzed Cowbird	Molaen	0	1	0	0	0	0	3	0	1	0	1	0	6
Molothrus ater	Brown-headed Cowbird	Molate	0	1	0	1	0	1	0	0	0	0	0	1	4
Myiarchus cinerascens	Ash-throated Flycatcher	Myicin	1	0	0	0	0	0	0	0	0	0	1	0	2
Myiarchus tyrannulus	Brown-crested Flycatcher	Myityr	4	0	2	0	2	1	1	2	3	1	4	3	23
Nyctidromus albicollis	Common Pauraque	Nycalb	0	0	0	0	0	0	0	1	0	0	0	0	1
Nycticorax nycticorax	Black-crowned Night Heron	Nycnyc	0	0	0	0	1	0	0	0	0	0	0	0	1
Nyctanassa violacea	Yellow-crowned Night Heron	Nycvio	0	0	0	0	0	0	0	0	0	1	0	0	1
Ortalis vetula	Plain Chachalaca	Ortvet	0	4	3	0	0	0	3	0	0	4	2	1	17
Petrochelidon pyrrhonota	Cliff Swallow	Petpyr	0	0	2	0	0	0	1	0	0	0	0	0	3
Phalacrocorax sp.	Cormorant species	Phasp.	0	0	0	0	0	0	0	0	0	0	0	1	1
Pitangus sulphuratus	Great Kiskadee	Pitsul	2	1	2	3	3	3	0	4	2	3	4	6	33
Psittacara holochlorus	Green Parakeet	Psihol	0	0	1	0	0	0	0	0	0	0	0	0	1
Quiscalus mexicanus	Great-tailed Grackle	Quimex	2	4	5	1	1	2	3	6	1	1	3	1	30
Setophaga pensylvanica	Chestnut-sided Warbler	Setpen	0	0	0	0	0	0	0	1	0	0	0	0	1
Spiza americana	Dickcissel	Spiame	1	0	0	0	0	0	0	0	0	0	0	0	1
Sturnus vulgaris	European Starling *†	Stuvul	0	0	0	0	0	0	0	0	0	0	0	1	1
Thryomanes bewickii	Bewick’s Wren	Thrbew	0	3	0	0	0	0	1	0	0	0	0	0	4
Thryothorus ludovicianus	Carolina Wren	Thrlud	0	2	3	1	0	6	2	0	3	0	5	1	23
Toxostoma longirostre	Long-billed Thrasher	Toxlon	0	5	4	5	5	2	5	4	4	3	0	6	43
Turdus grayi	Clay-colored Thrush	Turgra	2	0	0	0	0	0	0	0	0	0	0	2	4
Tyrannus couchii	Couch’s Kingbird	Tyrcou	6	3	5	4	3	6	3	6	4	6	3	2	51
Tyrannus melancholicus	Tropical Kingbird	Tyrmel	2	0	0	0	0	0	0	0	0	1	0	0	3
Vireo griseus	White-eyed Vireo	Virgri	0	1	0	1	0	0	2	1	2	1	0	1	9
Zenaida macroura	Mourning Dove	Zenmac	6	5	6	5	2	5	6	6	4	6	4	4	59
Site Total			72	63	56	48	42	45	60	65	51	63	57	61	683

), which is roughly only one-third of the 150 bird species documented at the Laguna Atascosa NWR that are known to use thornforest habitat in the spring or summer [54]. The five most common species were Zenaida macroura (Mourning Dove), Arremonops rufivirgatus (Olive Sparrow), Melanerpes aurifrons (Golden-fronted Woodpecker), Tyrannus couchii (Couch’s Kingbird), and Toxostoma longirostre (Long-billed Thrasher). Sixteen species (30.2%) were encountered during only one survey. Table 4 and Table S6 list the numbers of surveys that detected each bird species at each study site or sampling point, respectively.

An NMS ordination of the observed bird communities is shown in Figure 5. Myiarchus cinerascens (Ash-throated Flycatcher), Spiza americana (Dickcissel), Tyrannus melancholicus (Tropical Kingbird), and Lanius ludovicianus (Loggerhead Shrike) were associated with two of the older, larger, passively restored sites, whereas Cathartes aura (Turkey Vulture) and Toxostoma longirostre (Long-billed Thrasher) were associated with three other older, larger sites. Higher abundances of Zenaida macroura (Mourning Dove) and Cardinalis cardinalis (Northern Cardinal) were associated with half of the actively restored sites and with greater distances to temporary water, soil moisture content, and hog disturbance. The other actively restored sites were associated with greater abundances of Leptotila verreauxi (White-tipped Dove) and higher canopy tree densities. However, the PerMANCOVA results showed that only distance to temporary water and canopy tree density were significantly associated with bird community composition (

Table 5. Results of multivariate (PerMANCOVA) and univariate (ANCOVA using Type III sums of squares) analyses examining the effects of site traits on bird community composition and species richness. The effect size is the predicted change in a univariate response variable given a unit increase in the predictor variable. More complex models were considered for each response prior to model pruning; factors not included herein were pruned due to insignificance or to avoid multicollinearity. Full results tables for individual models can be found in the supplemental tables shown. Legend: d.f., degrees of freedom; **, 0.001 ≤ p < 0.01; *, 0.01 ≤ p < 0.05; ., 0.05 ≤ p < 0.1.

Response	Predictor	Test	Effect Size	d.f.	F	p		Results Table
Composition	Distance to temporary water	PerMANCOVA	-	1	2.05	0.0010	**	S7
Composition	Canopy tree density	PerMANCOVA	-	1	1.65	0.0450	*	S7
Richness	Canopy plant diversity	ANCOVA	+2.465	1	7.15	0.0119	*	S8
Richness	Method of restoration	ANCOVA	−1.136	1	9.42	0.0044	**	S8
Richness	Invasive grass cover	ANCOVA	+0.041	1	6.16	0.0187	*	S8
Richness	Understory plant richness	ANCOVA	+0.399	1	3.64	0.0657	.	S8

and Table S7).

3.3.2. Bird Community Univariate Analyses

Bird richness averaged 13.56 ± 2.26 species and was influenced by canopy plant diversity, method of restoration, and invasive grass cover, whereas understory plant richness had a marginal effect (Table 5 and Table S8; Figures S11–S13).

We could not analyze bird diversity or abundance because we could not count or estimate the abundance for all species in all the surveys. Relationships between bird richness and canopy plant diversity, method of restoration, and invasive grass cover are shown in Figure 6, Figures S10 and S11, respectively. Canopy plant diversity and invasive grass cover had positive linear relationships with bird richness. The effect size for the restoration method suggests actively restored sites had 1.14 more species than passively restored sites, which is slightly greater than the difference between the observed mean values (passive had 13.06 ± 2.56 species versus active with 14.06 ± 1.86 species).

3.4. Lepidoptera

3.4.1. Multivariate Analyses of Lepidoptera Communities

We sampled 1641 Lepidoptera in total and observed 77 different species (eight butterflies and 69 moths), which are over half (59%) of the 130 Lepidoptera species documented at the Laguna Atascosa NWR [54]. Supplementary Tables S3 and S4 list the abundances of all the Lepidoptera species observed at each of our 12 study sites and 36 sampling points, respectively. None of the species observed were identified as invasive.

NMS ordinations of the observed Lepidoptera communities are shown in Figure 6 and Figure S14. Helia agna and Focillidia texana (Southern focillidia moth) were mainly associated with older, larger, passively restored sites, whereas a Renia species, Plusiodonta compressipalpis (moonseed moth), Bleptina caradrinalis (bent-winged owlet), and Elaphria chalcedonia (chalcedony midget) were strongly associated with younger, smaller, actively restored sites, as well as higher cover by exotic plants, especially invasive grasses. Moths in the subfamily Sterrhinae, Elaphria chalcedonia, and an Agonopterix species were associated with higher isolation and canopy cover. Metria bilineata, Libytheana carinenta (American snout), Melipotis indomita (indomitable graphic), Herminiinae morphospecies #7, Numia bicoloraria (bicolored Chloraspilates), Erebidae morphospecies #3, and Prochoerodes lineola (large maple spanworm) were associated with greater distances to permanent water and higher ground layer plant richness, as well as both methods of restoration.

Table 6. Results of multivariate (PerMANCOVA) and univariate (multiple regression and ANCOVA using Type III sums of squares) analyses examining the effects of site traits on Lepidoptera community composition, species richness, Shannon diversity (H’), and total abundance. The effect size is the predicted change in a univariate response variable given a one unit increase in the predictor variable. More complex models were considered for each response prior to model pruning; factors not included herein were pruned due to insignificance or to avoid multicollinearity. Full results tables for individual models can be found in the supplemental tables shown. Legend: d.f., degrees of freedom; ***, p < 0.001; **, 0.001 ≤ p < 0.01; *, 0.01 ≤ p < 0.05; ., 0.05 ≤ p < 0.1.

Response	Predictor	Test	Effect Size	d.f.	F	p		Results Table
Composition	Canopy cover	PerMANCOVA	-	1	3.01	0.0020	**	S11
Composition	Soil moisture content	PerMANCOVA	-	1	2.25	0.0050	**	S11
Composition	Method of restoration	PerMANCOVA	-	1	2.47	0.0060	**	S11
Composition	Invasive grass cover	PerMANCOVA	-	1	2.04	0.0220	*	S11
Composition	ln(native:exotic cover)	PerMANCOVA	-	1	2.07	0.0100	**	S11
Composition	Isolation	PerMANCOVA	-	1	2.64	0.0100	***	S11
Composition	Time since restoration	PerMANCOVA	-	1	2.47	0.0050	**	S11
Composition	Distance to temporary water	PerMANCOVA	-	1	1.56	0.0920	.	S11
Composition	Ground layer plant diversity	PerMANCOVA	-	1	2.20	0.0079	**	S11
Richness	ln(native:exotic cover)	Multiple regression	+1.320	1	6.54	0.0165	*	S12
Richness	Distance to permanent water	Multiple regression	−0.004	1	3.18	0.0860	.	S12
Richness	Canopy plant richness	Multiple regression	−1.311	1	5.33	0.0288	*	S12
Richness	Understory total cover	Multiple regression	+0.054	1	1.86	0.1838		S12
Richness	Patch size	Multiple regression	−0.023	1	2.77	0.1074		S12
Diversity	Soil moisture content	Multiple regression	+0.019	1	15.19	0.0005	***	S13
Diversity	Distance to permanent water	Multiple regression	−0.0004	1	3.29	0.0800	.	S13
Diversity	Invasive grass cover	Multiple regression	−0.004	1	2.46	0.1276		S13
Abundance	Soil moisture content	ANCOVA	−1.088	1	10.66	0.0030	**	S14
Abundance	Time since restoration	ANCOVA	+1.086	1	9.38	0.0049	**	S14
Abundance	Method of restoration	ANCOVA	+19.139	1	7.91	0.0091	**	S14
Abundance	Ground layer plant diversity	ANCOVA	−44.818	1	17.82	0.0002	***	S14
Abundance	ln(native:exotic cover)	ANCOVA	+15.230	1	14.07	0.0008	***	S14

and Table S11 show the PerMANCOVA results examining the effects of site traits on Lepidoptera community composition, as depicted in the ordinations in Figure 6 and Figure S14. The degree of isolation, canopy plant cover, soil moisture content, restoration method, natural log of the ratio of native-to-exotic plant species covers, time since restoration began, ground layer plant diversity, and invasive grass cover all had significant effects on Lepidoptera community composition; the distance to temporary water had a marginal association.

3.4.2. Lepidoptera Community Univariate Analyses

Lepidoptera richness averaged 13.97 ± 4.44 species across the sampling points and was influenced by the natural log of the ratio of native-to-exotic plant covers and canopy layer plant richness (Table 6 and Table S12; Figures S15a and S16), whereas the distance to permanent water had a marginal effect. Lepidoptera Shannon diversity (H’) averaged 1.95 ± 0.42 and was affected only by soil moisture content (Table 6 and Table S13; Figure S15b), although distance to permanent water had a marginal effect. Lepidoptera abundance averaged 50.64 ± 29.16 individuals and was influenced by soil moisture content, time since restoration began, restoration method, ground layer plant diversity, and the natural log of the ratio of native-to-exotic plant covers (Table 6 and Table S14; Figures S17–S21). Actively restored sites had significantly higher Lepidoptera abundances (53.28 ± 30.76 individuals) than passively restored sites (47.47 ± 27.83 individuals); however, the estimated effect size is considerably larger, with active sites predicted to have 19.14 more individuals (Table 6 and Table S14; Figure S19).

3.5. Herpetofauna

3.5.1. Multivariate Analyses of Herpetofauna Communities

Overall, we observed 83 herptiles representing 10 species: 4 lizard, 3 frog, 2 skink, and 1 snake species. For comparison, 33 herptiles have been documented at the Laguna Atascosa NWR, excluding marine and aquatic species [54].

Table 7. Total herptile abundance by species at each of the 12 study sites. These site-level values are the sums of the values quantified for the three sampling points within each site; supplementary Table S5 lists the abundance of each herptile species observed at all 36 sampling points. The abbreviations shown are the symbols used for the species in the ordinations. Legend: AU, Anacua Unit; AC, Arroyo Colorado Unit; DH, Duck Head; EU, Ebony Unit; FH, Fish Hatchery; GC, Garza-Cavazos; GI, Goat Island; LU, Longoria Unit; PB, Phillips Banco; TB, Tahuacal Banco; TU, Tucker Unit; VN, Villa Nueva.

Scientific Name	Common Name	Abbreviation	AU	AC	DH	EU	FH	GC	GI	LU	PB	TB	TU	VN	Species Total
Anolis carolinensis	Green Anole	Anocar	0	0	1	0	0	0	0	0	0	0	0	0	1
Anolis sagrei	Brown Anole	Anosag	1	0	15	0	18	0	0	0	1	1	0	1	37
Aspidoscelis gularis	Common Spotted Whiptail	Aspgul	0	0	0	4	1	0	0	0	0	1	1	2	9
Aspidoscelis laredoensis	Laredo Striped Whiptail	Asplar	0	1	0	9	0	0	0	1	0	0	1	0	12
Eleutherodactylus cystignathoides	Rio Grande Chirping Frog	Elecys	0	0	2	0	0	0	0	0	0	0	0	0	2
Hypopachus variolosus	Sheep Frog	Hypvar	0	0	0	0	0	1	0	0	0	0	0	0	1
Incilius nebulifer	Gulf Coast Toad	Incneb	0	0	0	0	0	2	0	0	0	0	0	0	2
Plestiodon tetragrammus	Four-lined Skink	Pletet	0	0	1	0	2	12	0	0	1	0	1	0	17
Salvadora grahamiae	Texas Patch-nosed Snake	Salgra	0	0	0	0	0	0	1	0	0	0	0	0	1
Scincella lateralis	Little Brown Skink	Scilat	0	0	0	0	0	0	0	0	0	0	1	0	1
Site Totals			1	1	19	13	21	15	1	1	2	2	4	3	83

and Table S15 list herptile abundances by species at the 12 study sites and 36 sampling points, respectively. Anolis sagrei (Brown Anole, an invasive species) was the most abundant by far, with 37 observations (44%), followed by Plestiodon tetragrammus (Four-lined Skink), with 17 observations (20%). We observed zero herptiles at 13 of 36 sampling points.

Figure 7 shows an NMS ordination of the observed herptile communities. Anolis sagrei, Incilius nebulifer (Gulf Coast Toad), Hypopachus variolosus (Northern Sheep Frog), and Plestiodon tetragrammus were associated with higher distances to temporary water sources and canopy plant richness. Aspidoscelis laredoensis (Laredo Striped Whiptail) was closely associated with higher ground layer plant diversity and the log of the ratios of native-to-exotic plant covers. Scincella lateralis (Little Brown Skink) and Aspidoscelis gularis (Common Spotted Whiptail) were mostly associated with older, larger, passively restored sites.

Table 8. Results of multivariate (PerMANCOVA) and univariate (multiple regression and ANCOVA using Type III sums of squares) analyses examining the effects of site traits on herptile community composition, species richness, and total abundance. The effect size is the predicted change in a univariate response variable given a one unit increase in the predictor variable. More complex models were considered for each response prior to model pruning; factors not included herein were pruned due to insignificance or to avoid multicollinearity. Full result tables for individual models can be found in the supplemental tables shown. Legend: d.f., degrees of freedom; ***, p < 0.001; **, 0.001 ≤ p < 0.01; *, 0.01 ≤ p < 0.05; ., 0.05 ≤ p < 0.1.

Response	Predictor	Test	Effect Size	d.f.	F	p		Results Table
Composition	Isolation	PerMANCOVA	-	1	3.13	0.0100	**	S16
Composition	Distance to temporary water	PerMANCOVA	-	1	3.63	0.0070	**	S16
Composition	Method of restoration	PerMANCOVA	-	1	3.45	0.0030	**	S16
Richness	Time since restoration	Multiple regression	−0.025	1	16.42	0.0007	***	S17
Richness	Ground layer plant diversity	Multiple regression	+0.899	1	16.85	0.0007	***	S17
Richness	Distance to permanent water	Multiple regression	−0.001	1	8.45	0.0094	**	S17
Richness	Understory plant richness	Multiple regression	−0.170	1	7.35	0.0142	*	S17
Abundance	Method of restoration	ANCOVA	+1.518	1	7.26	0.0148	*	S18
Abundance	Distance to permanent water	ANCOVA	−0.006	1	7.26	0.0245	*	S18
Abundance	Soil moisture content	ANCOVA	+0.067	1	1.21	0.2858		S18
Abundance	Invasive grass cover	ANCOVA	−0.055	1	3.34	0.0843	.	S18

and Table S16 show the PerMANCOVA results examining the effects of site traits on the observed herptile communities, as depicted in the NMS ordination shown in Figure 7. The degree of isolation, distance to temporary water, and method of restoration were significantly associated with herptile community composition.

3.5.2. Herpetofauna Community Univariate Analyses

Herptile richness averaged 0.39 ± 0.44 species and was influenced by the time since restoration, ground layer plant diversity, distance to permanent water, and understory plant richness (Table 8 and Table S17; Figures S23–S26). Mean herptile diversity (H’) was 0.39 ± 0.44 and was independent of all the site characteristics or environmental factors tested. Herptile abundance averaged (0.30 ± 0.67 individuals) and was influenced by the method of restoration and distance to permanent water, while invasive grass cover had a marginal effect (Table 8 and Table S18). Actively restored sites had higher herptile abundances (4.77 ± 3.22 individuals) than passive sites (2.10 ± 2.28 individuals).

3.6. Ensemble (Combined) Animal Communities

3.6.1. Multivariate Analyses of Ensemble Animal Communities

In total, we observed 158 animal species (77 Lepidoptera, 53 birds, 18 mammals, and 10 herptiles). Each species was given equal weight, so Lepidoptera (49% of the species) held relatively more influence than the other taxa. NMDS ordinations of the combined animal communities are shown in Figure 8 and Figure S27. Mephitis mephitis (Striped Skunk), Anolis sagrei (Brown Anole), Lesmone detrahens (detracted owlet), Marimatha nigrofimbria (black-bordered lemon moth), and Bendisodes aeolia most strongly drove separation in the direction associated with younger, smaller, actively restored sites and with higher levels of understory cover, hog disturbance, overall exotic plant cover, understory plant richness, canopy cover, and total plant diversity. Asterocampa clyton (tawny emperor), Dryobates scalaris (Ladder-backed Woodpecker), and Dicotyles tajacu (Collared Peccary/Javelina) drove separation in the direction associated with older, larger, passively restored sites and with higher values for the interior-to-edge ratio and distance to permanent water. Melanerpes aurifrons (Golden-fronted Woodpecker), Myiarchus tyrannulus (Brown-crested Flycatcher), and Leptotila verreauxi (White-tipped Dove) most strongly drove separation in the direction associated with greater habitat isolation and canopy density. Odocoileus virginianus (White-tailed Deer), Cardinalis cardinalis (Northern Cardinal), an unidentified Zale moth species, and Erebidae morphospecies 1 were associated strongly with lower habitat isolation and canopy density and weakly with higher soil moisture content.

Table 9. Results of multivariate (PerMANCOVA) and univariate (ANCOVA using Type III sums of squares) analyses examining the effects of site traits on ensemble community composition and species richness. The effect size is the predicted change in the response variable given a one unit increase in the predictor variable. More complex models were considered for each response prior to model pruning; factors not included herein were pruned due to insignificance or to avoid multicollinearity. Full results tables for individual models can be found in the supplemental tables shown. Legend: d.f., degrees of freedom; ***, p < 0.001; **, 0.001 ≤ p < 0.01; *, 0.01 ≤ p < 0.05; ., 0.05 ≤ p < 0.1.

Response	Predictor	Test	Effect Size	d.f.	F	p		Results Table
Composition	Time since restoration	PerMANCOVA	-	1	1.95	0.0010	**	S19
Composition	Soil temperature	PerMANCOVA	-	1	2.13	<0.0001	***	S19
Composition	Total plant richness	PerMANCOVA	-	1	1.55	0.0380	*	S19
Composition	Distance to permanent water	PerMANCOVA	-	1	1.99	0.0011	**	S19
Composition	Total exotic plant cover	PerMANCOVA	-	1	1.53	0.0338	*	S19
Composition	Understory plant cover	PerMANCOVA	-	1	1.40	0.0807	.	S19
Richness	Invasive grass cover	ANCOVA	−0.142	1	14.94	0.0005	***	S20
Richness	Method of restoration	ANCOVA	+2.689	1	8.22	0.0073	**	S20
Richness	Canopy plant richness	ANCOVA	−2.045	1	7.48	0.0101	*	S20

and Table S19 show the PerMANCOVA results examining the associations between site traits and ensemble animal communities. The time since restoration began, soil temperature, total plant richness, distance to permanent water, and total exotic plant cover significantly influenced ensemble community composition, whereas understory plant cover had a marginal effect.

3.6.2. Ensemble Communities Univariate Analyses

We had to simplify all the taxa into presence–absence format because our avian data were binary; thus, we could not analyze ensemble diversity or abundance. The ensemble richness averaged (31.44 ± 6.01 species) and was influenced by invasive grass cover, method of restoration, and canopy plant richness (Table 9 and Table S20; Figures S28 and S29). Actively restored sites had significantly higher richness (32.72 ± 6.17 species) than passively restored sites (30.17 ± 5.73 species) (Figure S28b).

4. Discussion

In this study, we surveyed the animal communities of four taxa (mammals, birds, Lepidoptera, and herptiles) within 12 restored Tamaulipan thornscrub forest habitats in the Lower Rio Grande Valley (LRGV) of southernmost Texas. We compared our observed richness values to a reference list of species observed in the Laguna Atascosa NWR to obtain a sense of the potential species richness in the region. We did not expect to observe all the listed species, given (a) the size of the NWR (44,515 ha) compared to those of our study sites (5–175 ha), (b) the land use histories of our study sites compared to that of the NWR (relatively undisturbed), and (c) our sampling protocols, which were limited in scope and sometimes biased toward a subset of species. Nevertheless, even recognizing that missing species are not necessarily locally extinct, there is value in assessing these shortfalls and identifying potential candidate species for assisted recolonization.

We hypothesized that site characteristics influenced the composition, richness, diversity, and abundance of animal communities and that the nature and strength of these drivers would vary based on the taxon and ecological context. We found that each focal characteristic used for site selection (patch size, restoration time, isolation, and restoration method) significantly influenced at least one response metric, but most were less influential than expected. Other site characteristics, such as ground layer plant diversity, distance to permanent or temporary water, and soil temperature, influenced many community metrics and sometimes explained more variance than our focal site traits. Of our focal site traits, restoration time and method were most often significant, followed by isolation and then patch size. However, as discussed above and below, some focal factors were unavoidably confounded, so we cannot rule out the importance of those site traits, like patch size.

We consider our results in the context of several broader ecological concepts core to ecological restoration, namely, island biogeography, community assembly and succession, resource availability, and abiotic filters. In doing so, we aim to better understand the mechanisms underlying the observed patterns and link our results to the growing body of ecological research that continues to refine our understanding of these broader concepts.

4.1. Island Biogeography Concepts

Our study sites can be considered as thornforest islands surrounded by a matrix of non-forest land cover types that are separated from other thornforest islands to varying degrees. Habitat islands are more permeable than true islands, so predictions of island biogeography, like a positive species–area relationship [55], do not always directly translate [56]. We found no significant species–area relationship for any taxa, although patch size did influence mammal community composition (

Table 10. Table summarizing relationships between site traits and the community composition (C), species richness (R), diversity (D), and total abundance (A) of the four faunal communities surveyed plus the combined ensemble community (“all” below). The factors shown include all the site characteristics and environmental variables used in our analyses. Diversity and abundance could not be evaluated for birds or the ensemble community. Legend: ***, p < 0.001; **, 0.001 ≤ p < 0.01; *, 0.01 ≤ p < 0.05; ., 0.05 ≤ p < 0.1.

Factor	Mammals				Birds		Lepidoptera				Herptiles				All
	C	R	D	A	C	R	C	R	D	A	C	R	D	A	C	R
Time since restoration	.						**			**		***			**
Patch size	***
Extent of isolation		*					***				**
Method of restoration	***					**	**			**	**			*		**
Interior-to-edge ratio
Distance to permanent water	*							.	.			**		*	**
Distance to temporary water	*				**		.				**
Soil moisture content							**		***	**
Soil temperature		*													***
Canopy plant cover							**
Canopy stem density					*
Canopy plant richness								*								*
Canopy plant diversity						*
Understory plant cover			*												.
Understory plant richness		***				.						*
Understory plant diversity		*
Ground layer plant richness		**	*
Ground layer plant diversity	.			*			**			***		***
Total plant richness															*
Total plant diversity
Invasive grass cover						*	*							.		***
ln (native cover/exotic cover)							**	*		***
Total exotic plant cover															*

). It is possible that patch size did influence richness but was overshadowed by restoration time or method, which correlated with patch size. Alternatively, patch size may not be important for birds and Lepidoptera because flight allows them to utilize multiple patches more easily or not for herptiles because most have smaller home ranges than medium and large mammals.

The negative effects of habitat fragmentation and reduced patch size on mammal richness and diversity are well documented [57,58] but contrast our findings, and their effects on birds are mixed. Some studies found that patch size affected birds’ functional diversity [59] or abundance [60], but others found that factors like tree presence superseded patch size [61], and patterns vary among major climate zones. Our results aligned with those of Mellink et al. [61], who found that canopy tree density influenced bird richness but that patch size had no effect. These relationships likely depend on the identity of the species present and their individual requirements [60,61]. A future study could use our data to analyze the impacts of site traits on individual bird species to help resolve this uncertainty. For Lepidoptera, another south Texas study found that patch size did not influence richness or community composition but that site characteristics, like plant diversity, did [27]. We also found that several other site traits influenced Lepidoptera richness and composition more strongly than patch size (Table 10). Our finding that patch size did not affect any herptile community metric disagrees with much of the herptile literature, which often found that larger patches contained higher herptile richness, diversity, or abundance or otherwise differed in community composition [62,63].

If the effects of patch size operate more at the species level, they should depend on species-specific factors, like ecological needs, body size, mobility, and home range. Consideration of these factors for individual species is merited but beyond the scope of this study. The scale of the patch sizes considered (5–175 ha) is below the range many studies consider as large and may limit our capacity to detect the effects of patch size. Lastly, the influence of patch size may be masked by correlated factors that explain more variance, namely, the time since restoration and restoration method. Further insight may be gained by separately analyzing larger, older, passive sites and smaller, younger, active sites.

Island biogeography theory also posits a negative species–distance relationship; richness is lower on islands farther from a colonization source [55]. We considered the degree of isolation, rather than distance, because so few thornforests remain, and their capacities to function as a species source are unclear. We found that isolation influenced mammal richness and Lepidoptera and herptile community compositions (Table 10). For every 1% increase in non-thornforest land cover within 1 km of a site center, mammal richness increased by 0.04 species; i.e., 1 more species per 25% isolation (Table 3). This positive relationship conflicts with predictions of island biogeography, and the underlying mechanism is unknown. Isolated patches might have fewer predators, allowing for a greater variety of prey species to proliferate. Alternatively, isolation could have a concentrating effect, with thornforest-dependent species in an area having fewer options for where to forage or den.

Similarly unclear are the mechanisms driving the relationship between isolation and community composition for Lepidoptera and herptiles but not for mammals or birds. This pattern is likely related to dispersal ability; herptiles are generally less able to move between patches than birds or larger mammals. Landscape distribution also matters; common reptile species present in more areas (more sources) were impacted less by isolation [64]. However, Lepidoptera, being capable of flight, are less subject to dispersal challenges related to isolation [65] and may respond more to habitat quality than patch size or isolation [66]. We found evidence of the latter; isolation affected only composition, whereas several factors related to habitat quality (canopy cover, invasive grass cover, ground layer plant diversity, canopy plant richness, and log of the ratio of native-to-exotic plant covers) influenced composition, richness, and abundance (Table 10). Yet, we cannot rule out dispersal. Larger Lepidoptera generally have greater dispersal ability [67], and many Lepidoptera surveyed were very small, except that this pattern is often reversed for frugivorous butterflies [68], which our bait traps attracted. When community metrics have negative relationships with isolation, animal communities tend to be more vulnerable to habitat loss. We did not see any such relationship for our focal taxa, but that could be because sensitive species were already lost.

4.2. Assembly and Succession Patterns

Many studies on animal community assembly and succession in forests have found that, after disturbance, richness and diversity increase over time as dominance shifts from habitat generalists to forest specialists [69]. A meta-analysis found that vertebrate richness and diversity increased at later successional stages, but recovery was influenced by site traits and varied by taxa [70]. We expected diversity to generally increase with higher restoration time, but this was not the case. However, restoration time was associated with differences in the compositions of the mammal (marginally), Lepidoptera, and ensemble communities (Table 10). Restoration time also influenced Lepidoptera abundance and herptile richness. The mechanisms underlying these relationships are relatively clear, but why bird and herptile compositions lack such relationships is unclear.

Flight may give birds dispersal capabilities sufficient to recolonize more rapidly than other taxa, but the independence of herptile composition from restoration time is harder to explain. Our surveys may have been less likely to detect specialist species present later in succession, which would mean we mainly observed generalist or cosmopolitan species present in all successional stages. Alternatively, plant characteristics (especially those of ground layer and understory vegetation) may influence herptile composition more than time (Table 10). Yet, restoration time did influence herptile richness. It is possible that late-succession specialist herptile species are more likely to be cryptic, or at least more difficult to detect. Given our low observation rates, if even one species was missed at older restoration sites, it could explain the weakly negative relationship between herptile richness and restoration time.

Restoration methods have often impacted animal community assembly or succession [15,19,40]. In this study, the restoration method influenced the mammal, Lepidoptera, and herptile community compositions, and active restoration increased bird and ensemble richnesses and Lepidoptera and herptile abundances (Table 10). Ensemble richness was higher in actively restored sites, but bird richness was higher in passively restored sites. Passive restoration has been effective in supporting bird diversity in some cases [71], but other studies found that active restoration was more effective than passive in restoring bird communities [72], while still others saw little difference between active and passive methods [73]. At least one community metric for each focal taxon was influenced by restoration method. Our results reflect a frequent finding in the literature: The method of restoration is important, but its impacts vary by taxa. Our findings increase our understanding of the effects different restoration approaches have on animal communities in the LRGV of south Texas, but important questions remain, especially regarding mixed results observed among active restoration sites.

4.3. Resource Availability and Abiotic Filters

Finally, resource availability is fundamental to assembly and succession, and it has often been shown to impact faunal recolonization after disturbance. Factors like food and water availability can permit or preclude species establishment or persistence [74]. This is one reason why plant diversity promotes animal diversity. Ground layer plant diversity or richness influenced the mammal and Lepidoptera community compositions (Figure 3 and Figure 6) and promoted mammal and Lepidoptera abundances, mammal and herptile richnesses, and mammal diversity (Table 10). Invasive species are relevant here because they can displace native species that provide food and other resources for other species. We found that site metrics related to invasive species prevalence, which were overwhelmingly driven by introduced forage grasses, often had significant—but not always negative—relationships with animal community metrics and negative relationships with plant diversity metrics (Figures S2a and S3; Table 10). Many studies document the negative impacts of invasive plants on faunal communities [75,76,77], and this study provides more examples. Invasive animals were present but appear to be unimportant; relatively few species were documented, and they were always limited in abundance and distribution, except for Anolis sagrei in the herptile community (see species tables), and limited in their probable ecological effect.

Water availability was very important to animal communities. Distances to permanent and temporary (seasonal) water sources were frequently influential and impacted the community compositions of all the focal taxa (Table 10). For every kilometer of distance to permanent water, our models predicted herptile richness would decrease by 1 species and abundance would decrease by 6 individuals (Table 8). Distance to temporary water—which mainly represented proximity to critical seasonally flooded resacas (oxbow lakes) common in the region—influenced the community compositions of mammals, birds, and herptiles significantly and Lepidoptera marginally. The responses of birds to temporary, but not permanent, water could be because ephemeral shallow wetlands (e.g., seasonally flooded resacas) attract different species than deeper, more permanent water sources and/or because they can travel to more distant permanent water sources relatively easily if temporary sources are not available.

Soil moisture content was particularly important for Lepidoptera, significantly influencing their composition, diversity, and abundance (Table 6 and Table 10). However, it is unknown if this was due to short-term effects of the weather, which could drive higher bloom abundances, or long-term differences in moisture regimes. Studies show that soil moisture content can be very important to Lepidoptera, for example, by impacting larval pupation [78] or larval escape [79], and because chrysalis and caterpillar stages are often directly impacted by litter and soil conditions (about 25% of Lepidoptera interact directly with soil) [80].

5. Conclusions

We surveyed restored thornforests that varied in their times since restoration, restoration method, patch size, isolation, and various landscape and environmental characteristics. One or more of these factors significantly impacted every animal community metric quantified across four focal taxa, except for herptile diversity. The restoration method and time were the most influential focal site traits, affecting seven and five of the 16 community metrics analyzed, respectively, compared to patch size (one metric) and isolation (three metrics). However, patch size was confounded with both restoration time and method. Key landscape and environmental factors were distances to permanent and temporary water (six and four metrics impacted, respectively), ground layer plant diversity (five metrics), and invasive grass cover (four metrics). Soil moisture content influenced three metrics, as did the understory plant richness and the log of the ratio of native-to-exotic plant covers. Together, all 11 variables related to plant community structure influenced 19 animal metrics, and our three variables related to invasive species prevalence impacted eight metrics.

This study reinforces established notions that (a) restoration is not instant; (b) animal communities change over time in concert with plant communities as both progress through succession; and (c) invasive plants impact animal community composition, diversity, and abundance, though not always negatively for all taxa. The method of restoration was important, impacting animal composition, richness, and abundance, but taxa responded differently, so no one approach to Tamaulipan thornforest restoration is necessarily superior for all taxa. Core guiding principles, nevertheless, emerge: When restoring thornforests, invasive plant control, fostering native plant diversity, and ensuring there is a nearby water source are practical steps that can be taken to encourage faunal recolonization.

Future Directions

Many beneficial analyses could be performed using this dataset. Analyses of how individual species were affected by site variables, particularly conservation foci, like bobcats and ocelots, are worthwhile. A deeper understanding of the factors driving restoration outcomes might be gleaned by separately analyzing the actively and passively restored sites in this study. Doing so would show how the same site traits impacted restoration outcomes in the different groups and, crucially, would more fully disentangle the relative impacts of patch size, restoration time, and restoration method, which were unavoidably confounded in the current study. Explicit considerations of exotic versus native species and their patterns and prevalences across sites are merited. Analyses that consider animal functional groups, guilds, dietary niches, or other life strategies when examining the effects of site traits on restoration outcomes could improve our understanding of how individual species respond to restoration. Finally, there is value in investigating the effects of proximity to the United States–Mexico border and its barrier system (border wall), which impacts animal movement and, thus, may have important effects on restoration outcomes in the LRGV of southernmost Texas.