Land Cover Trends in South Texas (1987–2050): Potential Implications for Wild Felids

: The Rio Grande Delta and surrounding rangelands in Texas has become one of the fastest urbanizing regions in the United States over the last 35 years. We assessed how land cover trends contributed to the large-scale processes that have driven land cover change since 1987. We classified LANDSAT imagery from 1987 to 2016 to quantify different rates of land cover change and used housing density scenarios to project changes in the amount and spatial distribution of woody cover until 2050 and its potential impact on wild felid habitat. Since 1987, woody cover increased from 3.9% along with patch and edge density, whereas mean patch area and Euclidean nearest neighbor decreased. Closer inspection revealed that woody encroachment of small patches (<1 ha) was the leading cause of woody cover increase by a magnitude of 4, with an observed significant skewness and kurtosis in the frequency distribution of patch size across years. By 2050, urbanization will be the dominant landscape type and at least 200 km 2 of woody cover may be lost, thereby affecting felid populations in South Texas. These results provide important information for predicting future woody cover fragmentation and its potential impact on the connectivity of wild felid populations.


Introduction
Deforestation and degradation of native vegetation for agricultural land use and urbanization have had profound global effects on wildlife [1][2]. The human footprint on the natural landscape is unprecedented with more than 40% of the land surface affected by agriculture (i.e., crops and livestock) [1,[3][4]. The global extent of urban lands increased by 58,000 km 2 from 1970 to 2000, with the largest area of change occurring within North America [4]. Today, more than three billion people live in urbanized areas worldwide, with mid-century projections expected to surpass four billion [5][6]. As habitat fragmentation and degradation rates increase, the ability of certain wildlife species to survive or disperse into new areas will be affected [7].
In North America, the rate of urbanization and habitat fragmentation is currently outpacing the rate of land preservation and human population growth rate [5]. As habitats become fragmented, generalist carnivores such as coyotes (Canis latrans), raccoons (Procyon lotor), and red foxes (Vulpes vulpes) can exploit rapidly changing areas [5,8]. However, carnivores that have specialized habitat requirements such as the ocelot (Leopardus pardalis) and Canada lynx (Lynx canadensis) [9][10][11] and large home ranges (e.g., mountain lion [Puma concolor] and jaguar [Panthera onca]; [12][13][14] are often most at risk. In southern California, increasing urban development, increased traffic volumes, and increased road construction has been linked to genetic isolation and lack of gene flow in mountain lion populations [15][16]. Further, persistence and connectivity among jaguar populations have been severely limited due to deforestation caused by agriculture, human development, and rangeland conversion [7,14,17]. In the early 1900s, five felid species were known to occupy the Rio Grande Delta and adjacent rangelands in South Texas: ocelot, bobcat (Lynx rufus), jaguar, mountain lion and jaguarundi (Puma yagouaroundi) [18]. The population statistics for these species in the early to mid-1900s are not known, but all felids were extensively hunted and trapped [18]. Today, ocelots are now endangered, and less than 80 remain in two small isolated populations in Kenedy, Willacy, and Cameron counties in South Texas [10,21]. Despite continued opportunistic trapping and hunting on private rangelands, bobcats are the most abundant felid whereas mountain lions are rare in many parts of South Texas [19][20][21].
Ocelots, bobcats, and mountain lions in Texas are dependent on large patches of woody cover [19][20]22]. Bobcats are generalists, and occur in mixed forests, thornshrub, riparian corridors and floodplains, coastal prairies, pastures, and developed areas [8,19,23]. Ocelots occur across a wide range of ecosystems [24], but they are considered a forest interior species and prefer larger patches (mean patch area = 9.2 ha) of dense woody vegetation communities with 95% vertical cover and 85% horizontal canopy cover in South Texas [19,22]. However, unlike bobcats and mountain lions, ocelots are not commonly observed near urban areas in South Texas [22], possibly due to lack of natural corridors and preferred habitat types.
The Rio Grande Delta region, including the cities of Harlingen, Brownsville, and McAllen, are among the most rapidly growing urban centers in the United States [21,[25][26][27]. Historically, increases in row-crop agriculture, livestock production, and urbanization spurred the removal of native woody vegetation, which has transformed the Rio Grande Delta into the third most productive cropland in the United States [21]. An estimated 95% of native vegetation was lost during the early to midtwentieth century in the Lower Rio Grande Valley [28]. From 1930 to 1983, Cameron County lost 91% of its native woodlands [29]. From 1983 to 2003, native woodlands in Cameron and Hidalgo counties remained unchanged despite increases in urbanization [27]. However, no region-wide assessment of woody cover has been conducted since the mid-1980s.
Conversion of natural land cover (e.g., native rangelands) to cropland and pasture has been linked to losses in genetic diversity and population decline for ocelots across their northern geographic range [9,17,21]. Rapid expansions of road networks are leading to increases in roadcaused mortality of wildlife, especially felids in South Texas [21,30]. To help address these issues, conservation organizations, academic institutions, federal and state agencies, and private landowners often work together to identify drivers of habitat loss and to preserve habitat in certain areas where it will assist populations and increase gene flow [7,11,15,21].
Assessing land cover trends over the past 30 years will contribute to our understanding of the large-scale processes that have driven land cover change in the Rio Grande Delta and surrounding rangelands since the late 1980s. Long-term availability of LANDSAT imagery [27,31] combined with housing density projections [32] enables an estimate of past, present and future land cover change and future land cover change and their potential impact on felid populations in these areas. Our results will provide biologists and decision-makers with a baseline to develop landscape-level strategies to address habitat requirements for wild felids in Texas, as well as a wide variety of other species (e.g. raccoons, coyotes, and grey fox [Urocyon cinereoargenteus]). Therefore, the specific goals of this research are to: 1) Perform a land use land cover analysis to assess the extent and change of woody cover in the Rio Grande Delta and surrounding rangelands of Texas from 1987 to 2016.
2) Quantify the spatial and temporal distribution of woody cover to determine the extent of potential fragmentation that has occurred since 1987.
3) Predict future trends in land cover change from 2020 to 2050 based on housing density projection models and discuss its potential effects on native felid populations.

Materials and Methods
This study was focused in the Rio Grande Delta and surrounding rangelands of South Texas (i.e., southern Kenedy, Willacy, Cameron, southeastern Brooks, and eastern Hidalgo counties) (10,065 km 2 ) ( Figure 1). The major land cover types include low to high density urban development, urban green space (e.g., maintained parks and reserves), open water, barren land, pasture, crops, woody, and emergent herbaceous wetlands, thornshrub and live oak (Quercus virginiana) forests [33]. This region is located within the Laguna Madre Barrier Islands and Coastal Marshes, Coastal Sand Plain, and the Lower Rio Grande Valley ecoregion [34], which has a semi-arid subtropical climate (10° C to 36° C) [21]. The mean annual precipitation in eastern South Texas is 837 mm, with lesser amounts in the winter than the summer and greater amounts near the coast [35]. We classified urban and non-urban cover types in the Rio Grande Delta and surrounding rangelands for 1987, 1992, 2000, 2008, and 2016. Land cover types were classified using LANDSAT 5 satellite imagery for 1987, 1992, 2000, and 2008, and LANDSAT 8 satellite imagery for 2016 [22,27]. We downloaded one LANDSAT image during March and May of each year analyzed from the US Geological Survey Global Visualization Viewer. We selected this time of the year because the number of images with no cloud cover was the highest during this time period. We obtained U.S. Census Urban Area spatial data from the Texas Natural Resources Information System (TNRIS, Austin, TX) and then manually digitized urban areas in ArcGIS 10.5.1 to match the study area to account for differences in urban extent for each focal year (i.e., 1987, 1992, 2000, 2008, and 2016).
We conducted an unsupervised image classification in ERDAS IMAGINE (Hexagon Geospatial, Norcross, GA, USA) and classified the image into six categories: urban areas (residential, commercial, industrial areas and impervious surfaces); agricultural land (citrus, cotton, and other croplands); rangeland-herbaceous (rangelands, Texas Gulf Coast Prairie, grasslands and non-woody areas along the coast and within urbanized areas); woody cover (live oak forest, thornshrub, and forested wetlands); water (e.g., inland lagunas, streams, canals, reservoirs, estuaries, bays, tidal marshes, and the Laguna Madre); and barren lands (bare soil, beaches, inland sand dunes, mixed barren lands, oil fields, and exposed rocky areas) [36]. These categories were selected as they are relevant to wild felid studies in the region [20,22,37]. Further, ocelots and bobcats have been observed using woody communities next to coastal prairie, inland sand dunes, cropland, and urbanization [19,20,22,38]. An accuracy assessment was conducted in ArcMap 10.5.1 for each year analyzed based on 200 random points, where each image was considered satisfactory once accuracy reached 85% [22,36]. We first conducted an accuracy assessment with the 2016 image and verified it against knowledge of the landscape based on concurrent fieldwork in the region, high-resolution aerial imagery, and Google Maps following similar recently developed approaches [39,40]. Once we were able to identify the different categories in 2016, we used the same approach for the previous year and we used Google Earth, National Agriculture Imagery Program historical imagery to evaluate land cover in the study area. Once we obtained an accuracy >85% we classified the remaining imagery using 2016 as the reference image in terms of land cover types to identify land cover in previous years. We then ran the 200 random points for our accuracy assessment. We merged urban and non-urban rasters for each year analyzed into one mosaic image [36].
Land cover change analyses were conducted to determine the rate of change from 1987 to 2016 for urban areas, agricultural land, woody, and rangeland-herbaceous areas. We used the combine function in ArcMap 10.5.1 to examine trends and detect changes between each period and over 29 years. A new change map was generated to display the overall rate of change from 1987 to 2016 [41]. To assess changes in landscape structure, we conducted a spatial analysis in FRAGSTATS 4.2 and examined five class metrics: percentage of landscape (PLAND; %), patch density (PD; # patches/100 ha), edge density (ED; m/ha), mean patch area (MPA; ha) and Euclidean distance to nearest neighbor (ENN; m) for areas classified as woody, urban, rangeland-herbaceous and agriculture [7,22]. We also examined two patch metrics: patch size (PS, ha) and landscape shape index (SI) for woody cover for each year analyzed [7,42].
We tested whether the patch size and landscape shape index of woody areas were different between years by using a two-sample Kolmogorov-Smirnov Z test [42]. Following this, we compared means for each year to examine the mean patch size, standard error, and skewness and kurtosis for each yearly distribution. Skewness and kurtosis increased each year, thus, we ran a patch-level analysis and examined 98% of large and small patches. Frequency distributions were classified based on small patches (0.4-20 ha) and large patches (100-21,300 ha).
To examine future trends in the land cover change in the study area, we used five housing density scenarios, including one baseline [32]. Housing density scenarios were chosen because they integrate bi-national immigration/emigration, economic, population growth and environmental factors for each scenario, all which influence land cover change differently [32]. Scenarios used for this study include rapid economic development, slow population growth and high global integration (A1), slower rate of economic growth with reduced flow of people across international border regions (A2), globally-integrated world similar to A1 but greater emphasis on environmentally sustainable economic growth (B1), regionally-orientated world of moderate population growth and local solutions to environmental and economic problems (B2), and baseline (BC) [32]. Commercial and industrial development and low-high density housing types were merged into one urban cover type (32). Housing density scenarios were developed by Bierwagan et al. [43] based on a pair of axes: economic and environmental development and globalization vs. regionalization. Each scenario used a spatially explicit regional growth model, which estimated the quantity of additional residential units needed in each county to meet the demand specified by the growth model's population projections, which were based on a ratio of residences to census block populations. This dynamic growth model accounts for changes in demographic characteristics based on different scenarios and spatially allocated housing units compared to the county-level historic growth patterns, the spatial pattern of land ownership and travel time [43]. We merged these scenarios with 2016 classified imagery to quantify changes in landscape structure (PLAND, PD, ED, MPA, and ENN) for woody, rangeland-herbaceous, agriculture and urban land cover types every 10 years from 2020 to 2050. With these combined layers, we were able to quantify the amount of land cover lost for each land cover type between 2020 and 2050 in 10-year intervals.

Results
The mean land cover composition change was 36% for every 8-year period and 52% of the region remained unchanged from 1987 to 2016 ( Figure 2). The largest impact on woody cover was rangeland-herbaceous conversion. From 1987 to 2000, urbanization had the greatest impact on rangeland-herbaceous (30.8% change) and woody cover (4.1%). From 2000 to 2008, the greatest factor for woody cover loss and fragmentation was rangeland-herbaceous conversion (42.3%) and agricultural conversion (4.9%). From 2008 to 2016, the largest change detected was the agricultural conversion of rangeland-herbaceous areas and subsequent rangeland-herbaceous conversion of woody areas. From 1987 to 2016, woody cover increased from 1187 km 2 (11.9%) to 1519 km 2 (15.1%). Rangeland-herbaceous decreased from 4831 km 2 (48%) to 3563 km 2 (35.4%), agricultural area decreased from 2415 km 2 (24.9%) to 2013 km 2 (20.0%) and urban land cover increased from 493 km 2 (4.1%) to 1851 km 2 (18.4%). For woody cover types, ED and PD increased across the study area and ENN decreased from 201.7 m (SD = 142.2 m) to 82.0 (SD = 55.2 m) ( Figure 3). Urban land cover increased by 14.3%, but woody cover mean patch area decreased by 30.4% (Figure 3). Upon closer inspection, woody cover patch size decreased from 5.18 to 2.29 ha and shape index values decreased from 1.21 to 1.13 ha (Table 1). Woody cover structure was significantly different amongst years ( Table  2). Distribution of small woody cover patches (< 1 ha) increased from 7818 to 39,506 patches from 1987 to 1992, and to 59,335 patches in 2016 (Figure 4). From 1987 to 2008, 99.9% of large patches were <100 ha and this increased to 100% by 2016 (Figure 3).   0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 Of the five economic scenarios driving housing density, scenario A1 would have the greatest impact on woody cover habitat fragmentation and thus wild felid habitat in the Rio Grande Delta and surrounding rangeland (Figures 5, 6) In this scenario, urban PLAND is projected to increase by 6.1% from 2020 to 2050, which represents a 300% increase over the baseline scenario ( Figure 6). In the A1 scenario, agriculture will produce the greatest change in PLAND, and all cover types will decline, except for urban areas. Patch and edge density of all other cover types in A1 will decrease, while MPA and ENN will increase with the greatest changes observed in urban land cover ( Figure 6). The largest remaining patches of woody cover will be found on private lands in Kenedy and Willacy counties in scenarios A1, A2, B1, and B2, with smaller patches located in state parks, national wildlife refuges, and conservation tracts (Figure 7). However, scenario B2 is the least destructive scenario for landscape fragmentation and wild felid habitat in the Rio Grande Delta and surrounding rangelands ( Figure 6).

Discussion
Recent changes that have occurred in the Rio Grande Delta and surrounding rangelands suggest that natural vegetation land cover in the US-Mexico border area may be lost for good and may limit wildlife species that depend on them for habitat in these areas. Urban development has increased from 4.1% to 18.4% between 1987 and 2016 and by 2050 is projected to increase between 33.4% and 39.2%. This sharp increase in urban cover types over the next 30 years, will undo recent gains in woody cover that have occurred in small patches. The extent and distribution of woody cover in areas along the border are projected to continue to degrade and will likely become highly fragmented, with a few small exceptions, primarily due to urbanization and the impending United States-Mexico border fence, as well as from agriculture or rangeland-herbaceous conversion. Many of the publicand privately-owned natural areas along the border will likely be enclosed by urban development over the next 30 years.
Since 1987, woody cover has increased but agriculture conversion spurred by rapid urbanization along the United States-Mexico border has led to changes in woody cover landscape structure, particularly near urban areas in the southeastern portion of the Rio Grande Delta. We observed minor changes in landscape structure in the rangelands in the northern half as compared to the southern half of the study area. Natural revegetation and habitat restoration efforts from 1987 to 2016 may have contributed to sharp increases in small patches of woody cover since 1987, however, large patches (>100 ha) of woody cover with embedded patches of dense woody cover that are preferred by ocelots are now absent, which corroborates previous findings [10,22,44].
Woody cover in the northeastern part of the study area retains the largest remaining patches of woody cover (e.g., thornshrub and live oak forests) in the region. However, the lack of patches >100 ha indicates the need to increase the spatial distribution of woody cover in the region as this area is estimated to hold 80% of the current breeding population of ocelots in the United States. Increasing landscape connectivity and the size of large woody patches is critical for felid populations [7, 15-16, 38, 44]. Loss of large patches of natural habitat and reduced connectivity among preferred habitat patches have deleterious effects on felid populations in North and South America [12][13][45][46]. Incentives for ocelot conservation and recovery should be applied to increasing large patches of woody cover that may also support bobcats, mountain lions and other species such as gray fox, and coyotes that rely on large undisturbed habitat tracts [8,15,[45][46]. Further, these results and scenarios may be used as the baseline to support the formulation of effective conservation strategies for continued woody cover restoration for ocelots and other endangered species in the region.
We hypothesize that the Rio Grande Delta near the United States-Mexico border lacks adequate woody cover with embedded patches of preferred ocelot habitat [19] to support or enhance population viability [10,27,30] over the next 30 years. Our results indicate that gains in small patches of woody cover are potentially a result of woody plant encroachment, probably due to browsing/grazing regimes, fire, climate, or habitat restoration in the region [48]. Even if current trends continue, these gains will not be enough to offset the projected losses in woody cover due to expanding urbanization.
Losses in woody cover and embedded dense woody communities due to expanding urban development observed in the southern portions of the study area may have provided conditions for a precipitous decline in ocelots, as reported in earlier studies [9,10,27,49]. Janecka et al. [9] reported that the isolation of preferred habitat and fragmentation of surrounding woody patches has led to genetic divergence among ocelots in Texas and northeastern Mexico and a statistically significant difference in the genetic structure and diversity between the two ocelot populations found in the eastern part of the study area. To date, there are not been a single successful dispersal event between the two U.S. ocelot populations and at the current pace, urban expansion around the last remnant woody patches in the southeastern half of our study area may further reduce the likelihood of successful dispersal in the future [50]. Furthermore, bobcats occurring along the Rio Grande River have already exhibited genetic divergence compared to bobcats occurring in the northern half of the study area, likely due to reduced connectivity among habitat patches. These negative effects are likely to increase, as bobcats are sensitive to roads, elevated levels of human development, and require large patches of natural cover and riparian corridors for movement and dispersal [8,23,[44][45]51]. Habitat loss and fragmentation, whether related to expanding urban areas and road networks or woody and rangeland-herbaceous conversion to agriculture, will continue to form barriers to the dispersal and gene flow of wild felids [10,17,52].
Increases in human population, denser road networks, and related traffic volumes also will contribute to road mortality for ocelots and bobcats [21,[37][38]51]. Despite reaching high densities adjacent to urban areas, bobcats require large patches of natural habitat [8,12,23,45] and will often select against areas with roads within their home range [45][46][47]. Vehicle collisions in fragmented areas are an important cause of bobcat [45,52] and ocelot [30] mortality.
Five scenarios incorporating housing density projections from 2020 to 2050 combined with 2016 land cover classification predict increased habitat loss and fragmentation for wild felids over the same period. This study predicts that urban areas will replace agriculture as the dominant land cover type by 2050. Scenarios A1, A2, B1, and B2 predict that urban coverage in the study area (e.g., cities of Brownsville, Harlingen-San Benito, McAllen, Edinburg and the proximate communities) will merge into one extensive urbanized area. If combined with adjacent cities in the Mexican state of Tamaulipas, this urban expansion would create one of the largest metropolitan areas in North America [21]. Each scenario forecasts loss and fragmentation of woody cover, particularly in the southern half of the study area but unlike the change observed since 1987, the future driver of change will primarily be urbanization. Conversely, remaining large patches of woody cover within the northeastern rangelands of the study area are projected to remain intact, thus potentially providing a suitable condition for wild felid habitat, and partly counteracting/offsetting the effects of urban development and future woody cover loss.
Protection and creation of wildlife corridors have been occurring in the study area [21], providing corridors south towards the Rio Grande River. However, these will probably not benefit ocelot populations, as woody patches near the border are highly isolated by cropland and urbanization, which has prevented dispersal events despite two ocelots being documented along the Rio Grande River during the last 35 years [10,50]. Notably, the closest known Mexican ocelot population is in the Sierra Tamaulipas National Protected Area [53], which is more than 190 km south of the United States-Mexico border. As urbanization increases in the Rio Grande Delta, these corridors will likely become intermittent islands amongst a matrix of roads and development. Some corridors that are not isolated and maintain some conservation value may benefit from habitat restoration of areas once cleared for agriculture and rangeland to help increase connectivity and woody patch size. Additionally, we hypothesize that corridors located near urban areas may increase exposure to rodenticides [54] and will potentially increase the rate of infectious diseases from domestic animals to felids [55][56]. However, without access to large source habitats on each end of the corridors, these long linear vegetation communities may also serve as ecological traps or habitat sinks for ocelots, bobcats, and other mammals, thus increasing mortality [15]. Ecological traps and habitat sinks should be identified within the interior of the Rio Grande Delta and near the border with the Mexican state of Tamaulipas, which has also experienced intensive urbanization and agricultural conversion potentially exceeding what we observed in this study in the United States [17,21,57]. Therefore, highly urbanized areas on each side of the border are not suitable for sensitive or wide-ranging predators but may be used by generalist species such as raccoons and coyotes [8,17].

Conclusions
The Rio Grande Delta in South Texas is one of the fastest-growing urban areas. Humans have altered ecosystems in the region since the late 1800s, with rangeland conversion, agriculture, and recently, urban sprawl. Our study shows that habitat for wild felids, particularly ocelots and mountain lions, is no longer available along the United States-Mexico border within our study area. While habitat may be available in the northern rangelands surrounding the Rio Grande Delta, large tracts of habitat will be required to maintain the two viable ocelot populations in the US. Our results corroborate previous research and suggest that patch-level analyses used with class-level landscape metrics can be a powerful tool in understanding the underlying processes of land cover change, particularly with the urbanizing of wildlands.