Plant community predicts the distribution and occurrence of thick-billed grasswren subspecies (Amytornis modestus) in a region of parapatry

Additional keywords: ecology, geographical range, hybridisation, population distribution, population ecology.

Divergent evolution in ecologically heterogeneous environments is an important component of speciation (Schluter and Grant 1984; Schluter 2001). Ecologically heterogeneous landscapes are often linked to population divergence because they provide a mechanism for inhibiting gene flow between populations. For example, fragmented habitats can act as barriers to gene flow by preventing dispersal between populations (Dudaniec et al. 2011; Amos et al. 2012), and patchy habitats may promote selection for different traits that lower the probability of interbreeding (Ravinet et al. 2016). When diverged populations come into contact, gene flow may be enhanced through hybridisation, or reduced when reproductive barriers are reinforced. In either case, parapatric populations may develop novel evolutionary trajectories given (1) the formation of novel phenotypes through hybridisation (e.g. Kleindorfer et al. 2014), (2) genetic introgression via hybridisation (e.g. Borge et al. 2005; Bierne et al. 2013), (3) divergence as a result of selection (e.g. Chapman et al. 2016) and/or (4) the expression of reproductive barriers that finalise the speciation process (e.g. Beysard and Heckel 2014). Despite widespread evidence for the idea that different ecological environments select for different phenotypes, it is not well understood how ecological factors contribute to the occurrence of diverged populations within a region of parapatry. Contemporary ecological landscapes could hinder or favour particular patterns of gene flow affecting processes that lead to speciation.

Plant species with heterogeneous distributions across the landscape can reduce population connectivity among specialist organisms as resources are patchily distributed. There is growing evidence that fragmented habitats can disrupt fine-scale population processes such as mate choice or create isolated populations with limited gene flow at larger spatial scales (Athrey et al. 2012; Harrisson et al. 2013). Populations that occur in different habitats may have low gene flow even when the different ecosystems are in close proximity due to local adaptation (Cicero 2004). Some ecosystems are particularly prone to develop heterogeneous ecological landscapes. For example, the Australian arid zone has irregular soil types and limited and sporadic rainfall that creates patchy distributions of plant communities (Tongway and Ludwig 1990). Animals that live in the arid zone may be greatly affected by these habitat characteristics, particularly in areas where the distribution of different plant communities may be affected by changing climate (Ford 1987; Martin 2006). Patchy plant distributions are likely to affect gene flow between populations through high resource patchiness, even when the populations are in contact. Understanding how habitat changes across the landscape and whether regions of parapatry are associated with distinct habitat types will help explain how the ecological landscape can affect gene flow.

The thick-billed grasswren (Amytornis modestus) (TBGW), is a cursorial songbird of the family Maluridae (Rowley and Russell 1997). This species is endemic to Australia and is well known for its restricted geographic distribution (short-range endemism) typical of poorly dispersing organisms (Austin et al. 2013). The distribution of the TBGW has been reduced, probably through habitat disturbance (Garnett et al. 2011). Today, the TBGW is found in the arid regions of South Australia and New South Wales (Black et al. 2011). The species is geographically isolated from its closest taxonomic relative, the western grasswren (Amytornis textilis) (Christidis et al. 2010), which currently occurs in Shark Bay, Western Australia (A. t. textilis), and east of the Gawler Ranges and the adjacent plains of north-eastern Eyre Peninsula (A. t. myall) (Black et al. 2010). There are seven subspecies of TBGW, which are morphologically distinct in tail length and plumage colour (Black 2011, 2016). These subspecies cluster into two mtDNA ND2 haplogroups, Eastern and Western, that diverged ~0.36 million years ago (Austin et al. 2013). Two of the seven subspecies are extinct, one from the Eastern haplogroup, A. m. inexpectatus, and one from the Western haplogroup, A. m. modestus (Garnett et al. 2011). Consequently, there is one extant subspecies, A. m. indulkanna, in the Western haplogroup and four extant subspecies, A. m. raglessi (Vulnerable, IUCN), A. m. curnamona (Near Threatened, IUCN), A. m. obscurior (Critically Endangered, IUCN), and A. m. cowarie (probably Vulnerable), in the Eastern haplogroup (Garnett et al. 2011; Austin et al. 2013; Black 2016).

Allopatric divergence is thought to be a key process in the formation of TBGW subspecies (Austin et al. 2013). In a previous study, we found that two subspecies of the TBGW (A. m. indulkanna and A. m. raglessi) are currently parapatric between the South Australian salt lakes, Lake Eyre and Lake Torrens (Slender et al. 2017). In the west, A. m. indulkanna had a larger body and longer and narrower bill than A. m. raglessi in the east, while birds within the region of parapatry had either an intermediate morphotype (males) or morphology similar to that of A. m. indulkanna (females). These subspecies have diverged at the mtDNA ND2 gene by 1.7% (Austin et al. 2013), yet closer examination of haplotype distributions revealed that there were birds from A. m. raglessi that had a Western haplotype east of Lake Eyre and Lake Torrens (Slender et al. 2017). Similarly, there was mitochondrial and phenotypic discordance in the region of parapatry where females had both morphology of A. m. raglessi and eastern haplotypes and males had intermediate morphology and both eastern and western haplotypes. Gene flow between the subspecies appears to be asymmetric from A. m. indulkanna to A. m. raglessi, which may affect contemporary processes of divergence between these TBGW subspecies.

Ecologically heterogeneous landscapes may be an important factor affecting divergence and gene flow between grasswren populations because they are habitat specialists. Most grasswren (Amytornis) species and subspecies occur in distinct plant communities that provide a specific vegetation structure (Rowley and Russell 1997). There is a phylogenetic component to the occurrence of different grasswren species in habitat with different plant species (Christidis et al. 2010). The predominant habitat of the western grasswren subspecies, A. t. textilis, is Acacia spp. or Triodia spp. (spinifex) dominated shrublands (Brooker 2000), but for A. t. myall the predominant habitat type is blackbush or Australian boxthorn (Lycium australe) low shrubland (Black et al. 2009). The habitat type of the TBGW is low shrublands with variation in dominant chenopod shrub species (Black et al. 2011). A. m. indulkanna, found west of Lake Eyre and Lake Torrens, occurs typically in habitat dominated by Oodnadatta saltbush (Atriplex nummularia omissa) and cottonbush (Maireana aphylla), while A. m. raglessi, found east of Lake Eyre and Lake Torrens to the fringes of the Flinders Ranges, occurs in habitat dominated by blackbush (M. pyramidata) and low bluebush (M. astrotricha) (Black et al. 2011). Within the region of parapatry, a dune field stretches between Lake Eyre and Lake Torrens that may be a remnant of the Eyrean Barrier, a hypothetical biogeographic barrier between divergent populations (Serventy 1972; Schodde 1982). The predominant vegetation type found on the Lake Eyre–Lake Torrens sand dunes is sandhill canegrass (Zygochloa paradoxa). The extinct TBGW subspecies A. m. modestus is the only TBGW subspecies that has previously been recorded in this vegetation type (Black 2012). Therefore, these sand dunes and their corresponding vegetation type are likely to play a role in directing patterns of gene flow between A. m. indulkanna and A. m. raglessi.

This study aimed to determine whether habitat heterogeneity within a region of parapatry could affect gene flow between A. m. indulkanna and A. m. raglessi. We hypothesised that there is an association between habitat type and subspecies distribution that could affect the occurrence and direction of gene flow. We compared (1) plant community between the ranges of the subspecies outside the region of parapatry, (2) plant community within and outside the region of parapatry, and (3) plant community in the region of parapatry at sites with and without grasswrens. We expected that habitat in the region of parapatry would be suitable for both subspecies, which may have facilitated secondary contact and that TBGWs will not be found in habitat with sandhill canegrass, which may affect patterns of gene flow.

The two TBGW subspecies that are the focus of this paper occur only in South Australia: A. m. indulkanna is found to the west of Lake Eyre and Lake Torrens and A. m. raglessi is found to the east (Fig. 1). Because the distribution of each subspecies is so large, data collection was focussed within nine ~2500-km² areas that cover the subspecies’ range. Data were collected at a total of 104 sites within these areas (Table 1). The areas are grouped into three zones that delimit the overlapping and non-overlapping distributions of the subspecies. The non-overlapping distributions are previously described by Black et al. (2011) and contain birds with distinct morphology (Slender et al. 2017): Zone A (A. m. indulkanna) and Zone B (A. m. raglessi). The overlapping distribution is based on the occurrence of birds with intermediate morphology in the region of parapatry (Zone AB) (Slender et al. 2017). Zone A includes five areas (2500 km²) associated with A. m. indulkanna (n = 32 sites), Zone B includes two areas (2500 km²) associated with A. m. raglessi (n = 43 sites), and Zone AB includes two areas (2500 km²) in the region of parapatry (n = 29 sites) (Fig. 1). In Zone AB, data were collected from sites where TBGWs were present (n = 14) or absent (n = 15).

Fig. 1.  Map showing the location of sites where vegetation was sampled for A. m. indulkanna (Zone A), A. m. raglessi (Zone B) and their region of parapatry (Zone AB) in South Australia. The distributions of other thick-billed grasswren (TBGW) subspecies are excluded. Symbols indicate sampling sites where TBGWs were present (filled shapes) in Zone A, Zone B and Zone AB or absent (hollow circles) in Zone AB. The location of dune fields (grey shading) is representative of the possible distribution of sandhill canegrass (Zygochloa paradoxa). Zones (different symbols) comprise of 2500-km² sampling areas (numbers) that contain multiple sampling sites. Square: Zone A (1, Coward Springs Railway Siding; 2, William Creek; 3, Peculiar Knob; 4, Coober Pedy; 5, Oodnadatta), Circle: Zone AB (8, Stuart Creek Station; 9, Mulgaria Station separated by a line) and Triangle: Zone B (6, Witchelina Nature Reserve; 7, Mount Lyndhurst Station). The morphology of TBGWs in Zone A was distinct from the morphology of TBGWs in Zone B whereas in Zone AB, morphology of TBGWs was either that of A. m. indulkanna (females) or was intermediate (males) between the two subspecies. The unbroken line indicates the currently known distribution of the Western haplogroup while the broken line indicates the currently known distribution of the Eastern haplogroup. The inset of Zone AB shows sites where TBGWs were absent (open circle), sites where TBGWs were present but unsampled (star), sites where TBGWs were present and had an Eastern haplotype (filled black circle), and sites where TBGWs were present and had a Western haplotype (filled grey circle). Some absence sites are hidden behind presence sites and the circle size reflects the number of sampling sites within an area. Where there are multiple sampling sites within an area, the proportion of shading of the circle reflects the frequency of the two mtDNA haplogroups.

Table 1.  Coordinates of areas and sample size of vegetation transects within each location in Zone A, Zone B and Zone AB, where TBGWs were present or absent
ND, not done

Collecting data associated with TBGWs is particularly demanding because TBGWs occur in remote areas with limited infrastructure, and are shy and difficult to observe (Rowley and Russell 1997; Black et al. 2010). Surveys to record TBGW occurrence were conducted from 2012 to 2014, within the months from July to October. A peak in TBGW breeding activity has been recorded during these months in previous years (Black et al. 2011). Observations of TBGW absence were recorded only in Zone AB because we were interested in potential habitat fragmentation only within the region of parapatry. Our survey effort was targeted at sites (1) where there had been a previous sighting (Zone A and Zone B), (2) where the habitat was suitable (Zone A, Zone B and Zone AB), or (3) that were spaced at regular intervals (Zone AB). Sites with suitable habitat contained >10 plants within 1 ha that were a chenopod shrub species and that were larger than 0.5 m in either height or diameter. This description is based on habitat features previously described for TBGWs (Black et al. 2011). In large areas of suitable habitat and in Zone AB, sites were sampled by stopping roughly every 10 km along preselected roads or opportunistically if grasswrens were observed from the road.

To identify TBGW territories, we initially surveyed sites by visually scanning the ground and shrubs within 200 m with binoculars and recording whether grasswrens were seen or heard. We did this while systematically approaching large shrubs in order to flush concealed TBGWs. These surveys were performed across 1 ha and lasted 60 min. Following TBGW surveys, we additionally used playback of TBGW song recordings at sites where grasswrens were not observed. We used mist-nets to trap birds once TBGWs were identified for analysis of morphology and genetics for an earlier study (Slender et al. 2017).

In order to assess the detectability of TBGWs we subsequently revisited sites where TBGWs were initially detected in Zone A and Zone B, as well as at sites where TBGWs were either present or absent in Zone AB. Sites were revisited up to six times on different days in the same year as initially surveyed. We estimated detection probabilities and occupancy of TBGWs in the revisited sites using the program PRESENCE (MacKenzie et al. 2002). We used a single-season model and controlled for differences in detectability due to survey-specific conditions such as differences in wind strength, time of day, air temperature, observer (authors ALS and ML) and survey length. Using a subset of data from 2014 where TBGWs were present and that excluded opportunistic sightings, grasswren presence was confirmed by visual or aural detection on average (±s.d.) 27 ± 24 min (n = 33) from commencement of the survey.

At sites where TBGWs were present, we defined the core territory as the location where (1) there was a nest, (2) there was more than one independent grasswren sighting, or (3) a grasswren was first seen. Because most sites were visited multiple times and at least one visit was performed just after sunrise, this proved to be a good estimate of core territory as the same core territories were identified during a radio-telemetric study of TBGW home range at Witchelina Nature Reserve (Louter 2016).

We measured vegetation parameters from 2013 to 2015. Vegetation parameters were not likely to change significantly over this period, as the target plant species (long-lived perennial chenopod shrubs) are slow growing, even in the absence of grazing (Osborne et al. 1935; Crisp and Lange 1976; Crisp 1978). We used Jessup transect methodology along a 4 × 100 m transect line that covered a total area of 400 m² (Jessup 1951). Plants were measured across a 400-m² area to account for the spatial heterogeneity of arid-zone vegetation that occurs across patchy soil types and under conditions of variable rainfall (Williams 1982; Morton et al. 2011). For ease of measurement, the transect line was segmented into forty 5 × 2 m sample quadrants. We targeted the core of a TBGW territory, which is estimated to be 2 ha (Louter 2016), by placing the centre of the transect within the territory core, as defined above. Two transect arms proceeded 50 m north and 50 m south from the territory centre. The cover of all plants (not defined by species) as well as the cover of adult long-lived perennial shrubs (defined by species) was recorded. Adult chenopod shrubs were regarded as plants taller or wider than 10 cm. The variables recorded were plant species (for long-lived perennial shrub species), percentage cover of each target plant species and percentage cover of all vegetation.

The dominant plant species were identified separately for each independent variable (zone or TBGW presence). Dominant plant species within the presence variable (present/absent) were analysed only within Zone AB and dominant plant species within the zone variable (Zone A, Zone B and Zone AB) were analysed only across transects where TBGWs were present. Within each site (400 m²), the percentage cover of a plant species was calculated as a proportion of the total plant cover. The values for percentage plant cover were very small; the percentage plant cover as a proportion of the total plant cover was therefore more useful for statistical analysis and allowed us to control for possible differences in plant size between years and sites due to differences in rainfall. The average percentage cover of each plant species within an area (2500 km²) was then calculated. Only plant species with an average cover greater than 1% within an area grouped by the independent variable (zone or presence) were retained for further analysis. Dominant plant species were defined as plant species that were present in more than half the areas within an independent variable group (zone or presence).

To compare the dominant plant species between (1) zones and (2) presence/absence categories, we performed a one-way MANOVA with a post hoc Tukey test and calculated 95% confidence intervals with 1000 bootstraps. Many samples had zero percentage cover of a particular plant species, so the data were inversely transformed before the analysis. The transformed variables still contained some statistical outliers, so we also performed a non-parametric Mann–Whitney U test to determine the sensitivity of the results. The MANOVA (and confirmatory Mann–Whitney U tests) were used to compare the percentage cover of dominant plant species between (1) Zone A and Zone B (presence only), (2) Zone A and Zone AB (presence only), (3) Zone B and Zone AB (presence only) and (4) presence versus absence (Zone AB only). To reduce over-parameterisation of the data, we also calculated derived factor scores with a principle component analysis (PCA) using the percentage cover of the dominant plant species with varimax rotation. We performed an independent t-test to compare factor scores across the same four tests as above. All statistics were performed using SPSS Statistics 22.0 (SPSS Inc., Chicago, USA).

Occupancy analysis was performed with 81 revisited sites where TBGWs were previously recorded as present (Zone B = 37, Zone A = 29, and Zone AB = 15). Naïve occupancy for these sites was 0.53. When the probability of occupancy was constant and the probability of detection was adjusted according to the site-specific covariates, the average likelihood of detection/site survey was 0.46. The corrected probability of occupancy (±s.e.), taking detectability into account, was 0.84 ± 0.10.

Three components were extracted from a PCA of eight plant species that were identified as dominant in any of the independent variables (Table 2). These components all had an eigenvalue greater than one that cumulatively explained 63.7% of the variation within the dataset. The plant species that loaded heavily for Component 1 (PC1, eigenvalue 2.27) were Atriplex vesicaria (bladder saltbush) and Zygochloa paradoxa (sandhill canegrass). Lower PC1 scores indicated more percentage cover from A. vesicaria and less percentage cover from Z. paradoxa. The plant species that loaded heavily for Component 2 (PC2, eigenvalue 1.54) were Maireana aphylla (cottonbush), M. astrotricha (low bluebush) and M. pyramidata (blackbush). Lower PC2 scores indicated more percentage cover from M. astrotricha and M. pyramidata and less percentage cover from M. aphylla. The plant species that loaded heavily for Component 3 (PC3, eigenvalue 1.28) were A. nummularia ssp. omissa (Oodnadatta saltbush), Acacia spp. (usually A. tetragonophylla or A. victoriae) and Rhagodia spinescens (spiny saltbush). A lower PC3 score indicated more percentage cover from Acacia spp. and R. spinescens and less percentage cover from A. n. omissa.

Table 2.  Variable reduction for dominant plant species
The loading and model contributions of the percentage plant cover of eight dominant shrub species across three zones where TBGWs were present or absent are shown. Factor scores were calculated using inverse transformed variables with varimax rotation. Shrub species with the highest factor loadings are shown in bold

In Zone A and Zone B (TBGWs present), there were five dominant plant species. The plant species A. vesicaria, M. aphylla, M. astrotricha and R. spinescens were dominant in both Zone A and Zone B, while M. pyramidata was dominant only in Zone B. Post hoc comparisons of the MANOVA showed that percentage cover in three of five dominant plant species was significantly different between Zone A and Zone B (Table 3). Zone A had more cover from M. aphylla, but less cover from M. astrotricha and M. pyramidata compared with Zone B. The plant species A. vesicaria and R. spinescens had comparable cover in Zones A and B. We found similar statistical outcomes using the non-parametric Mann–Whitney U test, but with an additional significant difference showing more cover from R. spinescens in Zone A than Zone B (three of five comparisons P < 0.005, R. spinescens P = 0.017 and A. vesicaria, P = 0.100). To control for the effects of possible covarying plant species, we explored associations between zone and plant species using derived factor scores. We found the same pattern whereby the cover of dominant plant species in Zone A and Zone B was different. All of the three derived factor scores for cover of dominant plant species were statistically different between Zone A and Zone B, PC1 (t₇₃ = 2.168, P = 0.033), PC2 (t₇₃ = 5.207, P < 0.001) and PC3 (t₇₃ = 4.929, P < 0.001) (Fig. 2).

Table 3.  Vegetation cover (mean ± s.e. [95% confidence interval]) for six dominant shrub species in Zone A (n = 32), Zone AB (n = 14) and Zone B (n = 43), where TBGWs were present in South Australia
Vegetation cover was calculated using survey transects. The data are shown as raw variables with statistical analysis of the transformed variables using MANOVA. P-values <0.05 are shown in bold

Fig. 2.  Factor scores (mean ± s.e.) representing different plant communities across zones. The derived factor scores (see Table 3) are for Atriplex vesicaria and Zygochloa paradoxa (PC1), Maireana spp. (PC2) and Atriplex nummularia omissa, Rhagodia spinescens and Acacia spp. (PC3). All plant communities (PC1, PC2, PC3) differed significantly between Zone A and Zone B. Two plant communities (PC1, PC3) differed between Zone B and Zone AB, while Zone A differed from Zone AB in only one plant community (PC2).

In Zone AB (TBGW present) there were four dominant plant species. These were A. vesicaria, M. astrotricha and M. pyramidata (also dominant in Zone A and/or Zone B) and A. n. omissa, which was dominant only within the region of parapatry (Zone AB). Two plant species that were dominant only in Zone A and Zone B, M. aphylla and R. spinescens, were also included in the analysis. The MANOVA test showed that there was variation in cover between all five dominant plant species except A. vesicaria (Table 3). Post hoc comparisons showed there was more A. n. omissa in Zone AB than in either Zone A or Zone B. Zone A had more cover from M. aphylla than Zone AB while Zone B had more cover from M. pyramidata. Of the dominant species, four of six (67%) had similar cover between Zone A and Zone AB, and four of six (67%) had similar cover in Zone B and Zone AB. Apart from one plant species, the non-parametric tests showed a similar result to the MANOVA for comparisons between Zone A and Zone AB. There was less cover from M. astrotricha in Zone A than Zone AB (all P > 0.10 except A. n. omissa, M. aphylla and M. astrotricha, for which P < 0.05). Apart from one plant species, the non-parametric test also showed similar results to the MANOVA for comparisons between Zone B and Zone AB. There was more cover from R. spinescens in Zone B than Zone AB (all P > 0.10 except A. n. omissa, R. spinescens and M. pyramidata, for which P < 0.01). The derived factor scores were significantly different between Zone A and Zone AB for PC2 (t₄₄ = 3.458, P = 0.001), and between Zone B and Zone AB for PC1 (t₅₅ = –3.031, P = 0.004) and PC3 (t₅₅ = –4.966, P < 0.001) (Fig. 2).

There were six dominant plant species in Zone AB where TBGWs were either present or absent. The plant species M. astrotricha was dominant in both present and absent groups. The plant species Acacia spp. and Z. paradoxa (sandhill canegrass) were dominant only in the absent group. The plant species A. n. omissa, M. pyramidata and A. vesicaria were dominant only in the present group. The MANOVA test showed that there was more cover from Z. paradoxa and less cover from A. vesicaria, M. pyramidata, M. astrotricha and A. n. omissa where TBGWs were absent compared with where they were present (Table 4). The non-parametric test showed similar results except there was no difference in cover of M. astrotricha and A. n. omissa between presence and absence (three of six comparisons P > 0.05 and A. vesicaria, M. pyramidata and Z. paradoxa P < 0.05). Only the derived factor score, PC1, was significantly different between presence and absence (t_23.860 = –3.886, P < 0.001) (Fig. 3).

Table 4.  Vegetation cover (mean ± s.e. [95% confidence interval]) for six dominant shrub species in Zone AB (intervening zone) where TBGWs were present (n = 14) or absent (n = 15)
Vegetation cover was calculated using survey transects. The data are shown as raw variables with statistical analysis of the transformed variables using MANOVA. P-values <0.05 are shown in bold

Fig. 3.  Factor scores (mean ± s.e.) representing different plant communities where thick-billed grasswrens were present or absent in Zone AB. The derived factor scores for the plant communities with Atriplex vesicaria and Zygochloa paradoxa (PC1) differed significantly in Zone AB where thick-billed grasswrens were present and absent while the derived factor scores for the plant communities with Maireana spp. (PC2) and Atriplex nummularia omissa, Rhagodia spinescens and Acacia spp. (PC3) were similar.

The dominant plant species in the range of each grasswren subspecies (A. m. indulkanna and A. m. raglessi) and their region of parapatry differed. Where TBGWs were present in each core area there were different dominant plant species and the region of parapatry contained a mix of dominant plant species that were also present in either core area. The geographic distribution of the subspecies on either side of Lake Eyre and Lake Torrens suggests that the Eyrean Barrier may have isolated these subspecies in the past (Slender et al. 2017). The findings from this study suggest that both A. m. indulkanna and A. m. raglessi may be locally adapted to different plant communities that also may have been isolated in the past across the Eyrean Barrier. We suggest that secondary contact within the region of parapatry may have occurred because contemporary habitat within this region is suitable for both subspecies. The distribution of TBGWs in the region of parapatry was patchy, and we did not detect grasswrens at 52% of surveyed sites. The sites without grasswrens had more cover from Z. paradoxa (sandhill canegrass), a plant species that occurs mostly along sand dunes. Therefore, the distribution of sand dunes and their associated vegetation within the region of parapatry may be affecting contemporary patterns of gene flow.

Subspecies may be associated with different plant species because of heterogeneous plant distributions across the Eyrean Barrier. The TBGW subspecies A. m. indulkanna was found in habitat dominated by M. aphylla (cottonbush) whereas A. m. raglessi was found in habitat dominated by M. pyramidata (blackbush) and M. astrotricha (low bluebush) – a finding that mirrors previous studies (Black et al. 2011). Vegetation surveys show that M. aphylla and M. pyramidata occur in different regions of the Stony Deserts as they belong to different floristic groups. The associated floristic groups were M. aphylla/Eragrostis setifolia/Astrebla pectinata/Atriplex vesicaria (Low Very Open Shrubland) and M. aphylla/Atriplex nummularia nummularia (Low Open Shrubland) or M. astrotricha/Atriplex vesicaria/M. pyramidata (Low Open Shrubland) (Brandle 1998). If the predominant chenopod species on either side of the Eyrean Barrier were different, then local adaptation to different plant species may have facilitated divergence between historically isolated TBGW subspecies. The historical rise in abundance of chenopod species (saltbush and bluebush) and their isolation from their ancestral habitat type (Acacia and Eucalyptus) has been implicated in the divergence between the TBGW and western grasswren (Norman and Christidis 2016). It is possible that contemporary divergence of TBGW subspecies, as evidenced by different morphology (Slender et al. 2017), may be affected by the availability of different chenopod shrubs through a process involving ecological opportunity (e.g. Wellborn and Langerhans 2015).

Habitat clines create opportunity for hybridisation between parapatric populations. In this study, both A. n. omissa (Oodnadatta saltbush) and M. astrotrica (low bluebush) were found within the region of parapatry (Zone AB) at sites occupied by TBGWs. These plant species were also found in the range of A. m. indulkanna (A. n. omissa) and A. m. raglessi (M. astrotrica). We previously suggested that these two TBGW subspecies have made secondary contact and are potentially interbreeding in the region of parapatry (Slender et al. 2017). The correlation between the location of the region of parapatry and the presence of the vegetation types associated with both subspecies suggests that a mixed ecotype within the region of parapatry may have facilitated secondary contact and subsequent gene flow. Selection linked with particular plant species may be preventing morphological introgression outside the region of parapatry (e.g. Endler 1977). Correlations between environmental clines and the presence of hybrid zones are common in other parapatric species and are generally facilitated by low dispersal (Hollander et al. 2015; Patel et al. 2015) and/or low hybrid fitness (Alexandrino et al. 2005; Beysard and Heckel 2014). The distribution and availability of particular plant species may be an important factor mediating TBGW subspecies distribution and gene flow across the region of parapatry.

The fragmentation of particular plant communities within the region of parapatry may also affect gene flow by restricting dispersal of TBGW subspecies. In this study, the plant community in the region of parapatry at sites without grasswrens had significantly more Z. paradoxa than sites with grasswrens. Z. paradoxa is a specialised plant that is adapted to resource-poor and highly variable (stressful) environments such as sand dunes (Moseby et al. 1999; Roda et al. 2013). We previously showed that grasswrens with A. m. indulkanna morphology and western mitochondrial haplotypes were found east of the dune fields near the region of parapatry, but A. m. raglessi morphotypes and eastern mitochondrial haplotypes were missing on the western side of the dune fields (Slender et al. 2017). These findings suggest that Z. paradoxa may prevent dispersal of A. m. raglessi to the west but not A. m. indulkanna to the east, leading to greater genetic introgression into A. m. raglessi. We suggest that this is because A. m. indulkanna may be more tolerant of habitat with Z. paradoxa. The extinct sister lineage to A. m. indulkanna, A. m. modestus, was distributed in the MacDonnell Ranges of the Northern Territory and was known to occur in habitat with Z. paradoxa (Black 2012, 2016). There have also been opportunistic observations of A. m. indulkanna in Z. paradoxa near the southern shoreline of Lake Eyre South (A. Black, pers. comm., 2016). We suggest that A. m. indulkanna rarely occupies areas dominated by Z. paradoxa but Z. paradoxa may allow dispersal of A. m. indulkanna through this habitat type. However, it is possible that we have underestimated the occurrence of TBGWs in Z. paradoxa due to the low detectability of grasswrens. Asymmetric contact zones are common across avian populations with secondary contact (Dingle et al. 2010; Greig and Webster 2013; Beysard et al. 2015) and may be influenced by variable habitat distribution (Miller et al. 2014; Saura et al. 2014). Further studies that measure gene flow will reveal the extent and patterns of introgression between these two TBGW subspecies.

A large number of grasswren species and subspecies are in decline. Amytornis contains the highest proportion of threatened species and extinct subspecies within the Maluridae (reviewed in Skroblin and Murphy 2013). Habitat loss and fragmentation caused by overgrazing from livestock and feral herbivores are likely to be contributing factors (Schodde 1982; Reid and Fleming 1992; Garnett et al. 2011). Plant species respond differently to grazing, which can result in some plant species being more abundant and others less abundant (Fuhlendorf et al. 2001; Navarro et al. 2006). In this study, the region of parapatry and the core area for A. m. raglessi were associated with plant species (M. pyramidata) that prefer soil disturbed by grazers (Facelli and Springbett 2009). These ‘increaser’ plant species may have become more abundant since livestock grazing was introduced to the Stony Deserts over 100 years ago. The plant species A. n. omissa, which was found in the core area of A. m. indulkanna and in the region of parapatry, has low palatability to livestock and is likely to have become a dominant species after the introduction of feral grazers (Jessop 1995). Other plant species associated with grasswren habitat such as A. vesicaria (found in the core areas of both A. m. raglessi and A. m. indulkanna) and M. astrotricha (core area of A. m. raglessi) are likely to decrease in the presence of grazing because they are highly palatable or sensitive to grazing (Jessop 1995). Grasswrens may prefer habitat that has not been heavily grazed (Louter 2016) or habitat that contains ‘decreaser species’. A. m. raglessi is currently listed as vulnerable to extinction due to a reduction in available habitat (Garnett et al. 2011). Most of the arid rangelands of South Australia are used for pastoralism and are currently grazed extensively by domestic livestock. Conservation management for A. m. raglessi should focus on land-management practices and reducing the negative impacts of livestock grazing on long-lived slow-growing chenopod shrublands, which are critical for TBGW habitat.

Understanding associations between habitat parameters and gene flow across the landscape is important for conservation of populations and interpretation of evolutionary patterns (James et al. 1995). The TBGW is an important model for (1) planning conservation strategies because there are populations that are vulnerable to extinction and (2) understanding factors that affect hybridisation because there are subspecies that are parapatric. This study showed that gene flow across a region of parapatry may be affected by the distribution of different dominant plant species. This study offers insights into a model of ecological stepping-stones across fragmented habitat and ecological clines within a region of parapatry.

The authors declare no conflicts of interest.