Recovery of horse fly populations in Louisiana marshes following the Deepwater Horizon oil spill

The Deepwater Horizon oil spill in April 2010 had unprecedented impact on the Gulf of Mexico. We established the greenhead horse fly (Tabanus nigrovittatus Macquart) as a bioindicator of marsh health. This species is bound to coastal marshes, since its larvae develop as top invertebrate predators in the marsh soil. Immediately after the oil spill (2010–2011), populations of this horse fly declined in oiled areas of Louisiana marshes with significant impacts on genetic structure. In this follow-up study five years after the catastrophic event (2015–2016), we now report signs of recovery of populations in formerly oiled areas. Fly numbers increased compared to previous counts. Previously detected genetic bottlenecks in oiled populations have disappeared. Migration into oiled areas began to replenish formerly depleted horse fly populations in impacted regions with populations from non-oiled areas as an important source of migrants. Parameters of family structure that had been impacted by the oil spill (number of breeding parents, effective population size, number of family clusters) rebounded to levels similar to or exceeding those in non-oiled control areas.

. We combined samples that were not genetically differentiated and collected within the same location. This process resulted in eight genetically independent populations for 2016. For comparison, these genotyping data were analyzed together with data collected immediately after the oil spill 18 . Population genetic analysis via TESS assigned the individual tabanids from the 26 sample populations spanning the years from immediately after the oil spill (2010 and 2011) to six years after (2015 and 2016) to eight major genetic clusters (Fig. S1 in Supporting Information). The plot of the membership coefficients ( Fig. 2) that assign each individual horse fly to its predominant genetic cluster(s) visualizes the change in population genetic structure of T. nigrovittatus in oiled and non-oiled areas over a period of six years following the DWH oil spill.
The 2010 populations collected immediately after the oil spill of oiled and non-oiled areas fell into distinctly separated genetic clusters with none of the individuals assigned with membership coefficients of >80% to any region other than their own. In the first two years after the oil spill, the control populations were homogeneous suggesting high gene flow among those populations. In 2011, genetic signatures from the control populations (purple, Fig. 2) started to show up in the oiled regions indicating beginning migration from non-oiled to oiled regions, but not vice versa. Only one individual in one non-oiled population (RWR) showed affiliation to genetic signatures seen in oiled regions (gold, Fig. 2). In 2015, populations from formerly oiled regions showed still distinct genetic signatures with individuals belonging to genetic clusters not found in the control populations (cyan, Fig. 2). Similar to 2011, however, there was a considerable influx of genotypes typically found in the control areas (red, Fig. 2).
In 2016, essentially the same genetic signatures were found in individuals and populations from oiled and non-oiled regions (blue and cream, Fig. 2). Only some individuals from Elmer's Isle (EI) retained a somewhat distinct genetic signature (red and orange, Fig. 2). Populations from oiled regions in 2016 showed a greater degree of heterogeneity, possibly reflecting immigration from different genetic clusters (blue and cream, Fig. 2) from the control populations.
Genetic distances and migration rates. Genetic distances (F ST ) among non-oiled populations showed no significant difference among the years (all P ≥ 0.10, two-tailed difference of means test, 300 permutations, Table S1).  (Table S1). The genetic distances between populations from non-oiled vs. those from oiled regions, which were considerably high in the two years immediately after the oil spill, decreased after 2011 with 2016 F ST -values (0.163, SD = 0.095) being significantly lower than those from previous years (range 0.295 to 0.389, P < 0.001, Table S1).
Bayesian analyses demonstrated that recent migration rates (over the last few generations) among populations collected in 2015 and 2016 ranged from 0.1 to 25% with 68-99% being derived from the source populations in each generation (Table 2). This range is almost identical to what was recorded in 2010 and 2011 18 . However, there were changes in source and direction of migration patterns. Immediately after the spill event in 2010, only migration rates from and into non-oiled regions (see Table S4 in 18 ) exceeded the minimum value required for informative levels obtained from simulations in BAYESASS 21 . None of the populations from oiled areas showed meaningful levels of migration rates in 2010. Gene flow among geographically close oiled areas connecting Grand Isle and Grand Bayou populations with migration rates ranging from 0.026 to 0.261, was first detected in 2011, one year after the spill 18   Migration from non-oiled into oiled areas and vice versa was mostly below the detection threshold in 2010 and 2011 except for one non-oiled western population (SC) that contributed small proportions of migrants (m) to oiled areas in the east (GIP: m = 0.017 and GB3: 0.073) in 2011 18 . In 2015, however, migration from non-oiled regions into formerly oiled regions increased with CP being the foremost source population (into GIP: m = 0.082, into EI: m = 0.119). In 2016, migration was observed from most of the non-oiled regions into all formerly oiled regions with reciprocity in one instance (from oiled GI into SC/RWR: 0.065, Table 2).
The average emigration rate from non-oiled regions almost doubled in 2015 and 2016 and was significantly higher than in the two years after the spill (Table S2). Emigration from oiled areas spiked in 2015. In 2016, emigration from oiled areas was still higher than immediately after the oil spill but the difference was not significant due to large variance among populations (Table S2). Immigration into non-oiled areas also showed considerable variance among populations and, thus, did not change significantly across the years, except for an increase in 2016 compared to 2011 (P = 0.01).
The main sources of repopulation from non-oiled into oiled areas shifted across the years. In 2010, no migration from non-oiled into oiled areas was detected. In 2011, SC showed minor contributions in terms of migrant proportions to GIP (m = 0.017) and GB2 (m = 0.073) 18 . In 2015, increasing proportions of migrants from CP were detected in GIP (m = 0.082) and EI (m = 0.119). Finally, in 2016, migrant proportions ranging from 0.03 to 0.14 originating from non-oiled areas (SC/RWR, CP, and CYRC) contributed to all formerly oiled populations (EI, GI, GB) ( Table 2).
Previously confirmed genetic bottlenecks in oiled populations disappeared five years after the oil spill. In contrast to populations collected immediately after the oil spill 18 , none of the populations collected in 2015 and 2016 showed signs of recent genetic bottlenecks. No heterozygote excess resulting from a population crash induced rapid loss of allele numbers was detected under the three mutation models (IAM, TPM, SMM, Table 3), regardless whether the area was oiled in 2010 or not. On the contrary, the majority of the populations, including those from oiled areas, showed significant homozygote excess in at least one of the mutation models.

Comparison of family structure between populations from non-oiled and oiled areas collected across six years.
Year had no significant effect in the previous data sets from 2010/2011 (GLM) 18   Effective population size and number of parents contributing offspring to population samples were equally high in non-oiled areas across all years (Fig. 3). In the years immediately after the oil spill, the effective population size and the number of parents were significantly lower in oiled areas, consistent with the decline in adult populations (details and P-values in 18

Discussion
A holistic understanding of the progression of long-term effects of oil spills as well as the timelines for recovery of salt marshes is essential for assessing the resilience of the marshes' vital ecosystem services and fundamental ecological patterns and processes. Native infauna, i.e., the meio-and macrofauna developing in the salt marsh soil, are postulated to be the most informative indicators of ecosystem recovery after oil spills due to their important interactions as part of the food web. Here we discuss the observed recovery of the greenhead horse fly 5-6 years after the Deepwater Horizon oil spill in the framework of the multi-layered recovery of infauna with different life history, diversity and function. We postulated that species native to the marsh, like T. nigrovittatus, with carnivorous sediment dwelling developmental stages placed at the top of the invertebrate food chain, should be most sensitive 22 and, thus, most valuable bioindicators of marsh health. The DWH oil spill caused a severe drop in the numbers of horse fly larvae and adults in oiled areas, accompanied by underlying changes in population genetic and family structure 18 . This led to the question if, when, and how the process of recovery might mitigate or even restore the genetic make-up of tabanid populations colonizing or residing in oiled areas.
Six years after the oil spill, adult census numbers and the frequency of larvae being found in the soil were increased in formerly oiled areas compared to our findings immediately after the oil spill. Adult numbers as stand-alone, however, should be interpreted with caution. Census data often overestimate effective population size 23 . Trap catches vary due to seasonal fluctuations 18 , weather conditions and the occasional presence of horse fly predators such as dragonflies and bembix sand wasps, which cause horse flies to avoid traps 24 . Although increases in adult numbers in formerly oiled areas were small, the incidence rate of detecting larvae in the marsh soil clearly increased. The larva numbers (0.9 larvae per sample) in the formerly oiled Grand Bayou location reached the range of larva recovered from control areas. Since the larvae are cannibalistic, finding approximately one per sample would be expected in a productive habitat. The other formerly oiled areas (Elmer's Isle and Grand Isle) were heavily impacted by clean-up and reconstruction efforts, which may explain why the ecosystem has not yet rebounded to sustain peak larval development. Although there were no pre-event data to know the "normal" baseline number of horse flies in specific areas before the oiling, the rising census numbers in formerly oiled areas could be interpreted as the first sign of population recovery, especially since recovery subsequently was confirmed six years after the spill by genetic measures.
The observed increase in effective population size, number of breeders and number of family clusters in formerly oiled areas, which reached the level of control populations in 2016, supported the rise in census numbers and explains the mitigation of genetic bottlenecks recorded after the oil spill. The quick (≤5 years) recovery from bottlenecks of tabanid populations from oiled areas is expected to reduce negative consequences on evolutionary potential and resilience of the impacted populations 25,26 that are typically predicted and observed in slow recovering populations 27,28 . Most of the populations of 2015 and 2016 showed significant homozygote excess in at least one of the mutation models. While inbreeding is part of the life history of T. nigrovittatus 19 , the novel rise of heterozygote deficiency in populations of the formerly oiled areas likely indicates a reverse bottleneck effect, i.e., a recent population expansion. Population expansion in formerly oiled areas could be caused by 1. successful migration into oiled areas, 2. sufficient rate of reproduction in oiled areas to overcome residual larval mortality, if any, and 3. increased survival of larvae developing in formerly oiled areas due to decreased oil residue toxicity to larvae and/or increased availability of their food web.
1. In terms of migration capacity, tabanids rank as strong fliers and can disperse across considerable distances 18,29,30 . In the present study we recorded shrinking genetic distance (F ST ) over the years between populations from non-oiled and oiled areas and among populations from oiled areas, which suggests increasing migration rates compared to the low gene flow among populations after the oil spill. Therefore, we tested the hypothesis that populations from non-oiled areas were an important source for repopulation of formerly depleted oiled areas.
The lack of migration of a typically strong flyer immediately after the oil spill in 2010 could be explained by reduced dispersal capacity of tabanids from oiled areas (comparable to the reduced swimming performance in fish after sublethal oil exposure 31,32 ; and/or mortality of migrants into oiled areas, because flies were attracted by polarized reflections from oil sheens while searching for fresh water 33,34 . However, in 2011 (one year after the spill), some degree of migration was detected among geographically close oiled areas (Grand Isle and Grand Bayou) suggesting sufficient reduction of oil in the water and/or in the sediment to sustain survival of migrants. Still, larvae were almost entirely absent from formerly oiled areas in 2011 18 . The most important parameter in terms of recovery was the significant increase in migration from non-oiled areas into oiled areas over the years with 2016 showing migration from most of the non-oiled regions into all formerly oiled regions and even occasional contributions of migrants from oiled areas into non-oiled areas. This flow of migrants from populations of non-oiled regions (or closely related unsampled populations with the same genetic signature) into oiled areas was reflected in the increasing homogenization of population structure in the later years, which stands in contrast to the distinct genetic signatures of populations immediately after the oil spill in 2010.
These results highlight dispersal ability as an important driver of recovery. As a general rule, populations of organisms with rapidly dispersing life stages were the first to recover as soon as oil contamination and soil conditions were sufficiently improved to sustain survival. Similar patterns were also found in long-term colonization and succession studies of invertebrates in newly constructed salt marshes 35 . Taxa with dispersing larval stages, e.g. some annelids, copepods and nematodes, quickly colonized newly created marshes or re-colonized formerly oiled sites within three to five years to a level equivalent to healthy reference sites 16,35,36 . In contrast, slow recovery (>4 years) was shown for taxa that lack larval dispersal and those whose later stages are also poor dispersers due to their association with the marsh sediment 16,35 .
In comparison to these macro-and meiofauna invertebrates that are similarly native and tightly bound to the marshes as T. nigrovittatus, the dispersal and recovery of the latter was quicker than expected despite the fact that this tabanid has sediment-dwelling carnivorous (top of the food chain) larvae. The larvae are accomplished swimmers and can be observed hunting in the tidal zone, which likely facilitates dispersal at the larval stage to some degree (LF personal observation, https://www.youtube.com/watch?v=GAdm_JqAuoM). In addition, adults are excellent fliers, which explains the observed gene flow from unoiled into oiled areas as soon as one year after the spill.
2. While dispersal capacity was certainly a strong factor to facilitate population recovery of T. nigrovittatus, successful recolonization of formerly oiled areas also requires reproductive success. Invertebrates with slow recovery include the tube-dwelling oligochaetes, polychaetes and tanaids, whose relatively small offspring numbers develop confined within the maternal tubes 37 . In contrast, quick recovering invertebrates (annelids, nematodes, and copepods) produce and release high numbers of progeny into the water column 16 . Tabanus nigrovittatus is known to lay >200 eggs per batch and is an autogenous species, i.e., females can produce the first batch without the need for a blood meal and produce subsequent batches after blood meals have been obtained. The pursuit of a blood meal and subsequent oviposition would ensure dispersal of individuals and recovery of depleted areas. Genetic data suggest that males and females might be able to mate multiple times (polygamy) and have multiple offspring with each partner 18 . The high reproductive and dispersal capacity of tabanids should translate into rapid population growth and re-colonization of formerly oiled areas, but only if larvae can develop in the sediment. The observed rebound in effective population size of adults, number of successful breeders and families in formerly oiled areas required not only a certain degree of decontamination of formerly oiled sediments, but also sufficient secondary (heterotrophic) production and food web support for larval development.
3. The process of recovery of the food web had started bottom-up almost immediately after the oil spill setting the stage for the return of horse flies to oiled marsh regions. Microbial bioremediation of oiled marsh sites 9,10,38 paved the way for vegetation recovery. Even in heavily oiled areas with near complete plant mortality initially, some recovery of was reported six years after the spill albeit total live aboveground biomass was <50% of reference marshes. In moderately oiled areas above ground vegetation recovered within 2.5 years 3,39 . Most microalgae and 90% of meiofauna (e.g., nematodes, copepods and most annelids), which form the base of the food web in the marsh sediments, recovered within three years after the DWH oil spill following closely the re-growth of Spartina marshes except for certain polychaetes, amphipods, kinorhynchs, ostracods and gastropods, that had not reached levels of reference marshes after 4 years 15,16 . According to McCall and Pennings' study 17 , the terrestrial arthropod community, including mostly herbivorous insects, and intertidal crabs had largely recovered one year after the spill despite suppression by oil exposure a year earlier. The guts of tabanid larvae contain mostly insect DNA as revealed by a recent 18S metagenome sequencing study (unpublished data). The quick return of insects to oiled marshes was, therefore, vital for the tabanid recovery.
Tabanus nigrovittatus populations' sensitivity to oil contamination evidenced by the immediate population crash following landfall of DWH oil 18 and their recovery following in succession with regrowth of Spartina and recovery of 90% of the meiofauna 16 makes this species a valuable bioindicator to assess the process of marsh reconstitution (Table 4). To provide more detail about what is required for the greenhead horse fly to colonize newly established wetlands or recolonize marshes after a catastrophic event, we are currently describing the food web of T. nigrovittatus larvae via 18S rRNA gene metagenomic sequencing of larval gut contents.

Material and Methods
Collections. For population census data, adult female flies were collected using canopy traps baited with dry ice 40 every other week from April through October in 2016, at the same four locations with 4-5 traps per site as described before (Fig. 1 19 . The same microsatellite loci developed and used in previous population genetic studies of the greenhead horse fly 18,19 were screened as described in 19 before being employed to compare population genetic parameters in 2015 and 2016 to those measured after the oil spill (2010, 2011) 18 . Since a formerly monomorphic locus (7FA) 19 showed multiple alleles in several samples of 2016, the number of polymorphic loci was increased from 10 (2010, 2011, 2015) to 11 (2016). Locus characteristics, GenBank accession numbers, and other summary statistics have been published previously 18,19 .  43 . For exploratory analysis, the simulations were run 10 times with and without admixture for each Kmax ranging from 2-26 for 1,200 sweeps with 200 sweeps burn-in. The maximum number of genetic clusters (Kmax) was determined from changes in the deviance information criterion (DIC, Fig. S1) 44 . Since Kmax and cluster patterns were similar with both ancestry models, and population differentiation in previous studies made it less likely that the majority of individuals have recent ancestors from different populations, a full analysis was completed with 100 runs with 50,000 sweeps and 10,000 burn-in using the no admixture model. The 20 runs with the lowest DIC (20% filter) were averaged using CLUMPP 1.1.2 45 . Finally, the estimated membership coefficients of each individual multilocus genotype to each genetic cluster were plotted with STRUCTURE PLOT Version 2 46 .
Migration rates. Bayesian statistics were employed (BayesASS v. 1.3) to estimate the proportion of migrants (Table 2) and immigration and emigration rates (Table S2), because these techniques allow for deviation from migration-drift balance and Hardy-Weinberg equilibrium 21 . Preliminary runs were used to assure that delta values for allele frequency, migration rate and inbreeding stayed between 20-60% of the total chain length 21 . The final settings were 3,000,000 iterations, a sample frequency of 2000, a burn-in period of 999,999 and delta values of 0.25.
Detection of genetic bottlenecks. Immediately after a genetic bottleneck, allele numbers are reduced faster than heterozygosity. To detect signs of bottlenecks, genotypes of individuals in each population were tested for heterozygote excess across loci using a Wilcoxon sign-rank test under different mutation models (IAM, i.e., infinite allele model, TPM, i.e., two-phased model of mutation with 80% single step mutations and 20% multi-step mutations, SMM, i.e., stepwise mutation model), as implemented in BOTTLENECK v. 1.2.02 47 . In general, SMM and TPM are recommended when using microsatellite data since they are taking evolutionary relationships among alleles into account 47 . However, in our study over a limited geographical and time scale migration likely exceeded mutation as source of genetic variation and, therefore, we also presented results for the IAM, an implicit model that makes no assumptions on ancestry of alleles and origin of diversity. Mating structure and effective population sizes. Pedigree structure in each population was inferred using the full pedigree likelihood method implemented in COLONY v.2.0.3.1 48 . Program choice and parameter settings are described in detail in Husseneder et al. 18 . To infer the best configuration of relationship structure, populations were subjected to full maximum likelihood calculations 48,49 . Based on the Best Maximum Likelihood Configuration output of COLONY, we determined the number of inferred parents contributing offspring to each population (i.e., a measure of the breeding population size) and the number of family clusters in each population. Values were corrected for uneven sample sizes by calculating each for a sample size of 30 individuals genotyped per population. We also estimated the number of partners that individuals mated with and the number of offspring produced per parent. The percentages of full-sib and half-sib pairs in each population were inferred from sibship assignment plots 48 . Effective population size (N e ) was inferred with the sibship assignment method in COLONY 48 . The sibship assignment method showed superior accuracy for "real-world" data sets compared to other methods since it does not assume isolated populations or random mating and performs comparatively well even with small sample sizes and numbers of loci 50 .
General Linear Model (GLM) analyses were used to test for effects of "year" and "condition" (oiled vs non-oiled) on the variables of mating structure and population size. For each variable pairwise comparisons were performed between populations from oiled and non-oiled regions and between each year using two-tailed t-tests (SPSS, IBM Corp. Released 2012. IBM SPSS Statistics for Windows, Version 21.0. Armonk, NY: IBM Corp.). P-values ≤ 0.05 were considered significant and P-values up to 0.10 were considered marginal.