A Metapopulation Model to Assess Water Management Impacts on the Threatened Australian Lungfish and Neoceratodus forsteri

2. Materials and Methods

2.1. Study site

The Burnett River and tributaries account for approximately 40% of the known range of Australian lungfish, with an estimated total population of about 10,000 individuals [37]. Australian lungfish occur in the main channel and major tributaries of the Burnett River from the tidal barrage (Ben Anderson Barrage) upstream to Three Moon Creek near the township of Monto, 310 km from the mouth (Kind et al., 2008; Macdonald et al., 2012) (Figure 1). However, lungfish are much more common in the middle to lower reaches between the township of Gayndah and the tidal barrage, with the lower downstream section recognised as a key spawning area [34]. The river is heavily regulated and highly developed [40], with more than 40 structures (barrages, weirs and dams) [41] supporting water supply schemes [32,42]. Paradise Dam (37 m high), constructed in 2005, is the most recent infrastructure built and inundates approximately 24 km of riverine habitat when full [32]. The dam includes a stepped spillway on the downstream face where Australian lungfish mortalities have been reported to occur under a range of flow conditions [43], but especially flood events with 152 adult deaths reported in 2010 [43]. Likewise, an increase in population downstream of the dam was observed after two of the largest floods recorded in the Burnett River, suggesting a net downstream migration [37]. Overall, however, assessment of tag-recapture data suggested there was no significant short-term change in the population in the Burnett River since the dam was constructed [37]. Other barriers in the system include the tidal barrage dividing fresh from saline water which regularly leads to stranding of lungfish in the estuarine waters of the Burnett River [34,42].

To capture this structure in the distribution in the Burnett River catchment, we recognise three populations in this study:

• Population 3 (P3) – the riverine reach downstream of Claude Wharton Weir in Gayndah down to the upper impounded area of Paradise Dam (Adopted Middle Thread Distance (AMTD) 176–202.4 km).
• Population 2 (P2) – the water impounded area of Paradise Dam (AMTD 131.4–176 km).
• Population 1 (P1) – the riverine reach downstream of Paradise Dam through to Ben Anderson Barrage (AMTD 25.9–131.4 km)

Figure 1. Map of lungfish distribution in SE Queensland with the Burnett River study site and populations P1, P2, and P3 (highlighted).

Figure 1. Map of lungfish distribution in SE Queensland with the Burnett River study site and populations P1, P2, and P3 (highlighted).

2.2. The biology and ecology of Australian lungfish

An expert workshop with fish ecologists, managers and modellers, combined with knowledge from published and unpublished literature, provided an up-to-date understanding of Australian lungfish biology and ecology. The species is long-lived (77 years) [44] large-bodied (up to 1.5 m and 29 kg) [45] largely a sedentary species with restricted home ranges (1–1.5 km) and does not have an obligatory migratory life-history strategy, but may undertake larger scale (10s of kilometres) movements to seek out suitable spawning habitats between July and December or in response to elevated flow events [34].

The entire lifecycle of Australian lungfish is completed in freshwater [34] with a strong affinity for dense submerged macrophytes (>90% cover), considered critical habitat for spawning and recruitment [45]. Annual spawning occurs in late-winter (Austral) to late-spring (August to November) in shallow (still or flowing) waters [46] when water temperatures are 18–28°C [35,47], although may begin as early as July and extend to January [48]. Spawning occurs in pairs or groups with the large (3 mm diameter) loosely adhesive eggs encapsulated in a jelly capsule (total diameter of 1 cm after water hardening) [45], deposited singly, attached to plants and unlikely to survive if they fall to the substrate ) [45],. Spawning most often occurs at night, amongst dense aquatic macrophytes in still to fast-moving water [25,33,35,49,50]. Brooks and Kind [34] found that Australian lungfish egg densities increased with the density of macrophytes and the density of eggs was greatest with shorter macrophyte stands in shallower waters. Small fluctuations in the water level can impact on the viability of eggs, due to the species preference for egg deposition in shallow water. Conversely, eggs deposited in deeper calmer waters may be negatively impacted by deleterious effects like water quality, substrate composition, or light attenuation on macrophyte condition [36]. Flooding in river systems can remove aquatic macrophytes through scouring, leaving limited spawning habitat available [36] until recovery occurs some four or more years without further disturbances.

Recruitment is believed to follow a cyclical pattern of high or good spawning years interspersed with years of low or poor spawning. Pusey et al. [51] reported a one-in-five year ‘good’ spawning event occurring regardless of environmental conditions while Kemp [25] suggested that successful breeding seems to have occurred at intervals of around 20 years since 1900. Good spawning years have been linked with suitable antecedent water levels and baseflow conditions that benefit aquatic macrophytes, which are intricately linked with spawning [34]. However, breeding may be delayed or skipped if spawning habitat is disturbed [52] or unavailable [53]. There is evidence that successful spawning and recruitment within impoundment habitats is unlikely due to having steep bank profiles, fluctuating water levels, highly variable water quality conditions and anoxic soft substrates, resulting in conditions considered unsuitable for macrophytes persistence and creation of habitats that support egg development and the food and shelter for young fish [31,32,34].

2.3. Structure for Australian lungfish population model based on life-history analysis

Figure 2. A depiction of (a) the stage structure for female lungfish with associated estimated time transitions; and (b) the age network diagram of the structure adopted for the female only lungfish population model, where recruitment is a combination of fecundity and egg, larval and young-of-year survival (b).

Figure 2. A depiction of (a) the stage structure for female lungfish with associated estimated time transitions; and (b) the age network diagram of the structure adopted for the female only lungfish population model, where recruitment is a combination of fecundity and egg, larval and young-of-year survival (b).

We transformed stage-based life-history information of Australian lungfish into an age-based life-history (Figure 2), which allowed us to construct an age-based model that better describes the ecology of the species, and better aligns with mainly annual timescales for planning and implementation of management actions. The use of a stage- and age-based model was essential to capture the interaction of key threatening processes for this species that act upon different stages such as flooding (removes habitat for eggs/larvae/young-of-the-year), loss of fish passage (adult spawning movement) and spillway mortality (adults moving downstream with flows) [36,37]. To ensure complete capture of the birth and death process, it was decided to construct an 80-age class model where fish in the final 80-year-old age class do not survive. Adopting an 80-age class model was considered sufficient to capture the longevity of the species based on the maximum age measured in the wild (77 years) [54]. The 80 × 80 matrix (Figure S1), with a pre-breeding census [55], represents this construct, allowing the population growth rate to be calculated from the associated female only characteristic equation (1), where the sex ratio is 1:1 [45], and estimated vital rate parameters with a meaning and description of the parameters explained in the following figures, tables and text.

$$1=S_E{S}_L{S}_0{S}_1{S}_2{S}_3{S}_4{S}_5{S}_6{S}_7{S}_8{\displaystyle \sum _{j=10}^{j=80}{\lambda }^{-j}0.5F_j{\displaystyle \prod _{i=10}^{i=j}{S}_{i-1}}}$$

(1)

The matrix can also be used to perform sensitivity analysis and calculate reproductive values, and other mathematically derived analysis such as generation time and the net reproductive rate [56,57].

2.4. Parameter estimation

2.4.1. Age-specific survival

The age-structured matrix requires estimates of age-based survival rates and age-based fecundity as a function of recruitment to one-year-old fish. Age data is essential to generate estimates of age-specific survival [18,58,59]. The data from Fallon et al. [54] provided 158 aged Australian lungfish, with a maximum age of 77 (Figure S2). An age-frequency data set was compiled with a regression of the form $freq\sim {a/i}^c-b$ fitted to the data. This form produces an idealised age-frequency for all age-classes (not restricted to the data age-classes) from which survival rates were estimated as the ratio between idealised age classes, taken as the ‘average’ or mean survival rate [17,60]. Age-specific survival rates were estimated from the proportion of fish in one age divided by the proportion of fish in the previous age class to estimate $S_i$, i = 1…79 (see Table 1), age-specific survival rates, that ranged between 0.2331–0.9997 (Figure S3).

Early life-history survival, the survival of eggs $\left(S_{eggs}\right),$ larvae $\left(S_{larvae}\right)$ and young-of-the-year (YOY) juvenile fish $\left(S_0\right)$ in the wild are unknown. Egg survival was estimated to be 0.2 based on observations and repeated sampling of eggs over time (Steven Brooks, unpubl. data.), larval survival was assumed to be one fifth of egg survival, i.e. $S_{larvae}=0.04,$ which is similar to the estimate by Houde [61] of larval survival for freshwater fish of 0.05 and consistent with larval survival in other Australian freshwater fish [18,58,59]. Survival of YOY juvenile fish $\left(S_0\right)$, over the 10 month period prior to becoming a one-year-old fish, was estimated from the age specific survival equation in Table 1 with i = 3/12 (the first month as a YOY) and i + 1 = 12/12 (the final month as a YOY: note that setting i + 1 as the final month is equivalent to taking the product of each individual month to be one fifth of first year survival, i.e. $S_0=0.05209$ The combined survival of early life-history $S_0\times S_{larvae}\times S_{eggs}\approx 0.00042$provides approximately one in 2400 eggs hatch and survive to become one year old fish, similar to other long-lived Australian native fish [17]. The high survival rates from $S_{10} \mathrm{to} S_{79}$ accords with the expected high survival rates for adults [62].

Table 1. Life history and key parameters for the Australian lungfish metapopulation model.

Survival rates		Values pre-burn in
$S_{eggs}$	0.2	P1 upstream initial pop size	10,000
$S_{larvae}$	0.04	P1 upstream LTA	10,000
$S_0$	0.0521	P2 midstream initial pop size	12,500
$S_i=\left[a{\left(i+1\right)}^{-c}-b\right]/\left[a{i}^{-c}-b\right]$	i:1...79	P2 midstream LTA	12,500
a	208.97	P3 downstream initial pop size	7,500
b	-1.55	P3 downstream LTA	7,500
c	2.14	Burn in time steps	200
Coefficient of variation in survival rates		Post burn in (time step 200)
All survival rates	10%	P1 average initial pop size	10,083.79
Paradise dam spillway survival reduction	50%	P1 min initial pop size	819
(only occurs for WRD flow)		P1 max initial pop size	22,508
Fecundity		P2 average initial pop size	12,488.16
$Fe{c}_i=MF/\left(1+100\times \exp \left(-i+10\right)\right)$	i:10…80	P2 min initial pop size	1,116
MF	75,000	P2 max initial pop size	32,508
Coefficient of variation in fecundity	10%	P3 average initia7 pop size	7,552.16
Proportion of females breeding (PFB)	95%	P3 min initial pop size	398
PFB in Paradise Dam	2%	P3 max initial pop size	17,530
		Maximum movement rates
Single population model summary		P1 to estuary	0.005
Maximum age	80 years	P1 to P2 WRD flow	0.001
Annual growth rate (λ)	1.0587	P1 to P2 Natural flow	0.025
Generation time (G) (years)	46.24	P2 to P1	0.05
Net reproductive rate	8.13	P2 to P3	0.025
Average no. of eggs survives to adulthood	~8*10^-6	P3 to P2	0.025

Table 1. Life history and key parameters for the Australian lungfish metapopulation model.

Survival rates		Values pre-burn in
$S_{eggs}$	0.2	P1 upstream initial pop size	10,000
$S_{larvae}$	0.04	P1 upstream LTA	10,000
$S_0$	0.0521	P2 midstream initial pop size	12,500
$S_i=\left[a{\left(i+1\right)}^{-c}-b\right]/\left[a{i}^{-c}-b\right]$	i:1...79	P2 midstream LTA	12,500
a	208.97	P3 downstream initial pop size	7,500
b	-1.55	P3 downstream LTA	7,500
c	2.14	Burn in time steps	200
Coefficient of variation in survival rates		Post burn in (time step 200)
All survival rates	10%	P1 average initial pop size	10,083.79
Paradise dam spillway survival reduction	50%	P1 min initial pop size	819
(only occurs for WRD flow)		P1 max initial pop size	22,508
Fecundity		P2 average initial pop size	12,488.16
$Fe{c}_i=MF/\left(1+100\times \exp \left(-i+10\right)\right)$	i:10…80	P2 min initial pop size	1,116
MF	75,000	P2 max initial pop size	32,508
Coefficient of variation in fecundity	10%	P3 average initia7 pop size	7,552.16
Proportion of females breeding (PFB)	95%	P3 min initial pop size	398
PFB in Paradise Dam	2%	P3 max initial pop size	17,530
		Maximum movement rates
Single population model summary		P1 to estuary	0.005
Maximum age	80 years	P1 to P2 WRD flow	0.001
Annual growth rate (λ)	1.0587	P1 to P2 Natural flow	0.025
Generation time (G) (years)	46.24	P2 to P1	0.05
Net reproductive rate	8.13	P2 to P3	0.025
Average no. of eggs survives to adulthood	~8*10^-6	P3 to P2	0.025

2.4.2. Age-specific fecundity

The matrix construct requires age-specific survival and recruitment to one-year-old fish to be estimated using the equation: $R_i=S_0\times S_{larvae}\times S_{eggs}\times Fe{c}_i$, where $Fe{c}_i$ is age-specific fecundity and stage-specific survival estimates outlined above. It is likely that some females may mature early, and accordingly a sigmoidal function was used to calculate age-specific fecundity increasing from 10year olds to 80-year-olds, $Fe{c}_i=MF/\left(1+100\times \exp \left(-i+10\right)\right)$ (see Figure S4). Australian lungfish lay 200–600 eggs in a single event, it is assumed that fish are serial spawners laying for several days in a row and repeated during the season with the maximum total potential fecundity of 77,289 based on the average of three large running ripe female lungfish from the Brisbane River in September 2019 (David Roberts, Seqwater, unpubl. data). We assumed a maximum fecundity of 75,000.

2.5. Mathematical and sensitivity analysis

The characteristic equation (1) when parameterised with the mean fecundity and survival rates provides several summary population statistics on the status of the species. Stossel et al. [17] identify these statistics as population growth rate (λ) which summarises the capacity of the population to grow (or decline); the net reproduction rate $\left(R_0\right)$ which is the average number of female lungfish produced per female; and generation time (G) which measures the average age of parents of the current cohort [63]. Generation time is one of the measures used in categorising level of threat in the Guidelines for Using the IUCN Red List Categories and Criteria [63], and directly influences the time period of assessment and the period over which conservation efforts need to be considered [63,64].

Alongside these estimates, the use of sensitivity analysis allows for the exploration of uncertainties in the deterministic model and to assess this model for any disproportionate changes to the population growth rate due to uncertainty or changes in parameter estimation [65] (Todd et al., 2024). This approach reveals the importance of each parameter and its influence on model outcomes. Two types of perturbation analysis are used: (1) elasticity analysis, measuring the proportional rate of change in the population growth rate for a small change in a vital rate; (2) sensitivity analysis, measuring the absolute rate of change in growth rate given a small change in a vital rate; and we also calculate the age-specific reproductive values, measuring the contribution of an age-class to future generations, and summarising reproduction, survival, and timing [56].

2.6. Modelled flow effects on Australian lungfish

Conservation of Australian lungfish relies on maintaining sufficient macrophyte habitat for spawning, and for the survival of eggs and larvae. Smaller floods have stream velocities sufficient to remove macrophyte leaves, but allow for rapid regeneration, whereas higher flood flows can scour the bed of the stream and completely remove macrophytes including all propagules. Large flows can cause total loss of macrophyte cover, with adverse consequences for Australian lungfish spawning and recruitment success [36]. Recovery of macrophytes from large floods takes multiple years, regenerating from seed stock and propagules and can take up to four or more years to recover in minimally impacted catchments, longer in catchments with significant water impoundments.

In the Burnett River, flows ≥ 40,000 megalitres (ML)/day result in a proportional loss of macrophyte habitat based on observed losses of macrophytes following a range of flow magnitudes (Figure S5). This proportional loss of macrophyte habitat will remain for 365 days until natural recovery (or 3 years if the flow ≥ 155,000 ML/day – due to observed loss of all macrophytes and substrate) [36], unless one of the following happens:

• A larger flow occurs. In this case, the level of macrophyte loss will now be the corresponding level as per Figure S5, with the recovery period reset to 365 days (or 3 years, depending on size).
• A flow ≥ 100,000 ML/day, but smaller than the last flow to reduce cover. In this case the recovery period will be extended back to 365 days.
• flow ≥ 40,000 ML/day, but smaller than the last flow to reduce cover. In this case the recovery period will be extended to $floor\left(365\times \left(flo{w}_d-40000\right)/60000\right)$

, unless this is less than the current recovery period (in which case, the recovery period won’t change), where $flo{w}_d$ is daily flow data. It is thought that these flows would affect the recovery time of macrophytes as recovering macrophytes would be more sensitive to flows due to shallower root systems [66].

After the recovery period has ended, macrophyte habitat is reset to 100%.

In the model, the maximum potential recruitment from a breeding season (as a proportion) is based on the level of macrophyte habitat and the level of macrophyte habitat was determined by daily flow data as above. During the spawning season (August to November), in the model Australian lungfish have a daily probability of spawning as a parabolic curve ($P_d\left(spawn\right)$: Figure S6). Hatching occurs after 23–30 days and larvae take another 14–21 days while the yolk sack is consumed before becoming free swimming juveniles, all the while it is assumed they remain within the dense macrophyte habitat used for spawning [25,67]. Therefore, we assumed macrophyte habitat is critical for the 240-day period following spawning, the period that spans the most habitat sensitive life stage (eggs to YOY). To calculate the daily contribution to seasonal recruitment, the proportion of eggs spawned on day (j), is multiplied by the minimum level of macrophyte habitat from j until j+240. Summing this over each day (d) in the breeding season (L) gives the maximum potential recruitment based on macrophyte cover:

$$MP{R}_t={\displaystyle \sum _{d=1}^L{P}_d\left(spawn\right)\times {\displaystyle \prod _{j=0}^{j=240}MacrophyteCove{r}_{d+j,t}}}$$

(2)

Macrophyte habitat also influences survival of YOY fish, with declining effects on 1–4-year-old survival. The proportion of available habitat (macrophyte cover) was used to scale survival of YOY fish and reduced scale of 0.5, 0.25, 0.125, 0.0625 on survival of fish aged 1–4.

2.7. Density dependence

The ecology of Australian lungfish provides some direction for a suitable density-dependence mechanism. Juveniles have been rarely encountered in the wild and the recruitment of juveniles to the adult population is not well understood [51,68]. The relatively high survival rate of freshwater larvae, in comparison with marine larvae, supports the likelihood that the juvenile stage controls the level of recruitment [61]. Australian lungfish less than one-year-old are vulnerable to interspecific and intraspecific competition and changes in resource availability [45]; juveniles are thought to remain in their natal habitat for long periods [45]; and, in above ground tanks, yearling lungfish shelter in complex shaded habitats [45]. Kemp [68] suggests that low recruitment frequency, or episodic recruitment, is a normal strategy for such a long-lived species. These insights indicate that competition for resources may be high when there is strong recruitment, and that the competition is highest for YOY lungfish as they transition to one-year-old lungfish and this boundary is a plausible location for the impacts of density dependence (increased food demands and habitat requirements as fish grow). Young-of-year lungfish are much smaller compared with other age classes and are likely to be able to utilise habitat (food and refuge) not available to larger lungfish. As the total population would likely influence the strength of density dependence, density dependence was applied to YOY with their number proportionally decreased when the total population was above a weighted average population size (or long-term average in the absence of flow variation). The density dependence ratio of the average population size to total population size was weighted by a factor of 5 to account for size disparity (1 year old ~300 mm and adults up to 1.5 m [34,45]), similar to density dependence mechanisms used for other Australian freshwater fish [58,59].

2.8. Individual population stochastic population model construct

Using the 80 age-class projection matrix (Figure S1) as the structure, the pre-birth census female-only population model, where males were thought not to be limiting, is represented by the following equations:

$$\begin{array}{l}{N}_{i+1,t+1}=Bin\left(N_{i,t},S_{i,t}\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=1\dots 79\\ N_{1,t+1}=Poisson\left(Den{s}_t\times MR{P}_t\times S_{0,t}\times S_{larvae,t}\times S_{eggs,t}\times 0.5\times Egg{s}_t\right)\\ Egg{s}_t={\displaystyle \sum _{i=10}^{i=80}{N}_{i,t}\times Fe{c}_{i,t}}\\ Fe{c}_{i,t}=\mathrm{round}\left(MF\times \left(1+0.1\times fecsn{d}_t\right)/\left(1+100\times \exp \left(-i+10\right)\right)\right)\hspace{1em}\hspace{1em}i\ge 10\\ Den{s}_t=\{\begin{array}{cc}5\times LTA/T{N}_t,\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{for} T{N}_t>5\times LTA\hfill \\ 1,\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{for} T{N}_t\le 5\times LTA\hfill \end{array}\\ T{N}_t={\displaystyle \sum _{i=1}^{80}{N}_{i,t},}\end{array}$$

where t is an annual time interval. $N_{i,t}$ is the number of female Australian lungfish in the ith age-class at time t; $S_{i,t}$ is a random variate describing environmental variation in the survival rates in the ith age class drawn from normal distributions transformed to the unit interval [69] with no loss of information from the specified means and standard deviations; $S_{0,t}$, $S_{larvae,t}$ and $S_{eggs,t}$ are random variates describing environmental variation in survival rates of lungfish less than one-year-old, larvae, and eggs respectively, similarly to age-specific survival at time t, where survival rates across all ages and stages are perfectly correlated; $Den{s}_t$ is the density-dependence factor for one-year-old lungfish based on the total number of female lungfish at time t; $MR{P}_t$ is the maximum recruitment potential governed by flow at time t; $Egg{s}_t$ is the total number of eggs produced at time t from females aged 10 and older with full maturity from age 20 and older; $Fe{c}_{i,t}$ is the fecundity for lungfish in the ith age-class at time t, where age-fecundity is perfectly correlated but not with survival; MF is the maximum fecundity for lungfish and $fecsn{d}_t$ is a standard normal deviate at time t for variation in fecundity due to environmental variation; LTA is the long-term average total female population size (including juveniles); $T{N}_t$ is the total female population size (including juveniles) at time t; $Bin\left(n,s\right)$ is a random variate representing demographic variation in the transition from one age-class to the next drawn from a binomial distribution $Bin\left(n,s\right)=X\sim \mathrm{Binom}\left(n,s\right),$ and $Poisson\left(m\right)$ is a random variate representing demographic variation in recruitment drawn from a Poisson distribution $Poisson\left(m\right)=X\sim \mathrm{Poi}\left(m\right)$. Variation in survival rates is unknown, we assumed a coefficient of variation for all survival rates of 10%.

2.9. Burnett River metapopulation model

In the Burnett River, Paradise Dam poses a barrier to Australian lungfish movement and has fragmented what was once a continuous population into three partially connected populations [37]. We constructed a metapopulation model containing three populations that are allowed to move between populations (Figure 1) based on specific rules related to flow. Flow data used in the model for the Burnett River was based on modelled flow data from 1^st July 1890 to 30^th June 2009 extracted from the SOURCE program that is used by the Queensland Government to guide management of water resources particularly through the use of water plans [70]. The reason modelled streamflow data were used instead of measured streamflow is two-part: (1) the gauged flow data from each site is for a shorter time period (GS136007A Burnett River at Fig Tree since 1997; 1/3 the lifespan) and included missing data; and (2) the use of modelled data allowed the separation of flow alteration impacts for lungfish. Separation of the flow-alteration impacts was achieved by comparison of two specific modelled scenarios. The first modelled flow sequence was a Water Resource Development (WRD) flow which represents current conditions including all water storages in the Burnett water plan area, including Paradise Dam, and the full use of existing water entitlements, as well as improved environmental flow release rules from dams implemented in the 2014 water plan [70] (Figures S10, S11). The second modelled flow sequence was a Predevelopment (Natural) flow which represents a counterfactual flow and was modelled without water extraction or the impact of dams and weirs, including the absence of Paradise Dam. The length of both WRD and Natural flow sequences was 118 years, long enough for population-level impacts to become apparent.

2.9.1. Recruitment and movement rules

Recruitment

Recruitment in Population 1 and Population 3 were calculated using modelled streamflow data at two river gauge sites (GS136002D Burnett River at Mount Lawless: Figure S7; and GS136007A Burnett River at Fig Tree: Figure S8; respectively) [71], based on equation (2). For the WRD flow, it was assumed minimal recruitment occurred in Population 2 (MRP=0.02, based on failed spawning and recruitment in storages) [35] and for the Natural flow it was assumed the recruitment in Population 2 followed the same pattern as Population 1.

WRD flow movement

All downstream movement in the Burnett River model is a proportional movement of the population based on daily flow data, calculated as the sum across the year of daily movement, based on distributions for each population (Figure S9), multiplied by the proportion that hadn’t moved before that day, up to a maximum proportion dependent on the population. The maximum annual movement from population P3 to population P2 is 0.025. There is no movement between these populations during flows under 10,000 ML/day, with flows above this resulting in a daily movement between 0.005 and 0.025 (Figure S9). Because there are no physical barriers between populations P3 and P2, we assumed a minimum annual movement of 0.005, with the sum of daily movement added to this amount. The maximum annual movement from population P2 to population P1 is 0.05. There is no movement between these populations during flows under 40,000 ML/day, with flows above this resulting in a daily movement between 0.002 and 0.05 as in Figure S9 [72]. The maximum annual movement from population P1 into the estuary (past the Ben Anderson Barrage) is 0.005. No movement into the estuary occurs during flows under 8,600 ML/day, with flows above this resulting in a daily movement between 0.0001 and 0.005 (Figure S9). Australian lungfish rarely survive in the estuary [45] and are assumed to die following movement into the estuary.

Upstream movement from Paradise Dam was assumed to be negatively correlated with flow. The annual movement from population P2 to population P3 is calculated as 0.03 less the annual movement from population P3 to population P2 (giving between 0.025 and 0.005). There is limited upstream movement from population P1 as Australian lungfish were rarely reported using the Paradise Dam fish lift [43]. As such, we set movement from population P1 to population P2 at 0.001 for the WRD flow.

Natural flow movement

For Natural flow movement between population P3 and population P2 was used for movement between population P2 and population P1 (see Figure S9). For movement between population P1 and the estuary we used the downstream rule for WRD flow and an 80% return rate for juveniles and adults from the estuary to population P1 assuming no tidal barrage in place but some loss due to predation [45].

2.9.2. Initial conditions

The initial population size and age distribution represents the antecedent conditions of a population model. Age data presented in Fallon et al. [54] shows the age data to have significant gaps and it is likely that the antecedent conditions will have a mix of years where successful recruitment has occurred and years where little or no recruitment occurred [25,51]. We generated a random distribution of proportional age classes that sum to one, with a mix of zero (years of no recruitment) and non-zero years (Figure S10) and multiplied this distribution by the initial population size to generate an initial population distribution. The long-term average female population size (LTA) was set at a level to return an average female adult population size equivalent to the initial female adult population size. The model was run for 200 time-steps (years), a burn in period, to allow the model to escape the influence of the initial conditions, before the flow sequences were implemented. We used a 200 year burn-in as the species is long-lived. This approach has two additional benefits over escaping the influence of the initial conditions before the flow sequences were implemented from time step 201: 1) a variety of population structures are tested after time step 200; and 2) a variety of population sizes are tested after time step 200, creating a range of initial population sizes when the flow sequence is implemented. We set the initial female adult population size for P1 at 7,500, for P2 at 12,500, and P3 at 10,000. The total female adult population size at time step 200 was an average population size of 30,017.23 (see Table 1 for P1, P2 and P3 populations sizes at time step 200). From time step 201 we implemented the flow sequence and applied the associated rules for macrophyte loss and movement.

2.9.3. Other assumptions

Due to infertility and uncertainty about fish near the maximum age being reproductive, it was assumed that only 95% of adult fish produce fertile eggs [25]. Paradise Dam is a large impoundment with a stepped spillway and in high flows some lungfish moving with the flowing water are damaged by impacting on the spillway in overtopping events causing mortality. It was assumed that in high flows, the lungfish moving downstream from population P2 to population P1 would have a 50% reduction in survival due to high turbulence of spillway flows causing mortal injuries.

2.10. Model outputs – expressions of risk

All scenarios were run for 118 time-steps (the length of the WRD and Natural flow sequences) following the 200 time-step burn-in period, with 1000 iterations or trajectories to explore the underlying distributions of the stochastic parameters’ influence on population outcomes. One method to compare stochastic scenarios is to record the minimum population size from each trajectory and present these minimum populations sizes as a cumulative frequency distribution known as risk curves [55]. Risk curves represent the likelihood of extinction (probability of falling to zero) but also the likelihood of falling below some non–zero population threshold known as quasi-extinction risk [55] and has been used successfully on multiple species to document differences between the effectiveness of modelled scenarios on population risk [13,59,73]. Risk curves can be easily compared, a risk curve closer to zero represents a greater likelihood of lower populations sizes in comparison with a risk curve further from zero [73]. Associated with risk curves is the concept of quasi-extinction. Quasi-extinction risk is the likelihood of a population falling to, or below, some non-zero threshold that may be detrimental to the population if below or increase the likelihood the population is exposed to additional threats that it may not be able to recover from, such as behavioural, genetic or ecological [57]. Quasi-extinction risk is a key consideration in conservation management, with those species with reduced distribution range, reduced numbers or limited genetic diversity at a higher risk of extinction and quasi-extinction. Furthermore, as risk curves are a probability distribution we can also calculate the expected value of the risk curves, the expected value of the minimum population size distribution (EV) provides a quantitative comparison between risk curves [13,74]. The EV was recorded in the risk curve figures with the percentage difference between the Natural flow and the WRD flow. We also present plots of all 1000 trajectories, with summary statistics as well as plots of mean trajectories with shaded areas of standard deviation from the mean.

2.11. Model coding

The Australian lungfish model for the Burnett River was constructed in the software package ESSENTIAL [75]. This software and approach have been used extensively to model a variety of other aquatic species in Australia to aid in determining appropriate future management options [13,17,65,73,76]. The platform allows running of a metapopulation model and in this application to Australian Lungfish in the Burnett River includes the options of (a) Water storages scenario including the WRD flow sequence which includes Paradise Dam, limits fish movement upstream of the dam, includes a percentage of spillway mortality on flood flows as well as severely limits the macrophyte growth potential in the dam storage; and (b) running a Natural scenario including the Natural flow sequence and no water storages to impact movement rates. Each scenario was run for 1000 iterations over 318 time-steps with a burn in of 200 time-steps.

3. Results

3.1. Mathematical and sensitivity analysis

The Australian lungfish population growth rate was calculated as $\lambda =1.0587$ for the deterministic (matrix) model, the average number of female fish produced per female in her lifetime was calculated as $R_0=8.13$ and the generation time was calculated to be G = 46.24. Perturbation analysis shows that the population growth rate is relatively insensitive to proportional change (elasticity analysis) in vital rates, except for the aggregated adult survival rate which is not an individual age class ($S_{79}$: Figure S11a). The population growth rate is sensitive to additive changes (sensitivity analysis) in larval and YOY survival (Figure S11a). Survival of larvae and YOY fish are vital rates that are very hard to measure in the wild. Our best estimate for YOY survival rate comes from the age specific survival equation in Table 1 and may change if additional age data becomes available, however this would change all age survival estimates not just the YOY survival rate. While the model is most sensitive to the estimation of larval survival, doubling larval survival only produces a 2.6% change in the population growth rate. The reproductive value indicates that the highest value was achieved between the ages of 20 and 40 and declines thereafter (Figure S11b), a common trait for reproductive values [56], and is still relatively high around the generation time age, indicating the contribution that these older fish provide to future generations.

3.2. Model outcomes

For the Burnett River, the 118-year WRD flow data was used to analyse the impact of flow on recruitment and to model the response of population dynamics to changes in macrophyte habitat based on the flow rules identified in the methods. The results highlight that there are some differences in overall adult population trajectories between the three populations within the Burnett River using the WRD flow data (Figure 3). All populations exhibit a decline through the first 80 years, reflecting the lower recruitment opportunities through this period. Populations P3 and P1 (Figure 3a and 3c) recover after this period. However, population P2 (Figure 3b) continues to decline and remains low. Reduced spawning habitat is the key driver for the decline in population P2 (Paradise Dam) coupled with a net migration of lungfish downstream to population P1, albeit low level, due to limitations in the upstream fishway at Paradise Dam. The combined total population trajectories (Figure 3d) have a similar pattern to population P3 (Figure 3a) though recovery is proportional less due to population P2.

Comparing WRD flow with Natural flow produces different responses depending on the population (Figure 4). The mean trajectories and ± 1 standard deviation in populations P3 and P1 show little difference between the two flows (Figure 4a and 4e) which is reflected in the associated risk curves and expected value of the minimum population size EV (Figure 4b and 4f). In contrast, the mean trajectories and ± 1 standard deviation in population P2 show sharp differences between the two flows (Figure 4c) which is reflected in the associated risk curves and EV (Figure 4d) indicating a higher risk of a small population size under the WRD flow compared with the Natural scenario flow. At the metapopulation level the clear differences in population P2 between the flow scenarios are reduced (Figure 4g), though the risk curves indicate a higher risk of a smaller population size for the WRD flow compared with the Natural scenario flow (Figure 4h). This is an example of the benefits of constructing a metapopulation model to demonstrate the effects of flow changes at the local population level.

The expected value of the minimum population size distribution between WRD flow and the Natural flow: produces in population P3 a Natural flow EV 11% higher than the WRD flow (Figure 4b); produces in population P2 a Natural flow EV 179% higher than the WRD flow (Figure 4d); produces in population P1 a Natural flow EV 1% lower than the WRD flow (Figure 4f); and the metapopulation produces a Natural flow EV 29% higher than the WRD flow (Figure 4h).

Figure 4. Plots of mean female adult population size paired with the risk curves: population P3 (a) mean and standard deviation and (b) risk curve; population P2 (c) and (d); population P1 (e) and (f); and total population TP (g) and (h). Orange – modelled WRD flow and blue – modelled Natural flow: shaded areas are ± 1 standard deviation.

Figure 4. Plots of mean female adult population size paired with the risk curves: population P3 (a) mean and standard deviation and (b) risk curve; population P2 (c) and (d); population P1 (e) and (f); and total population TP (g) and (h). Orange – modelled WRD flow and blue – modelled Natural flow: shaded areas are ± 1 standard deviation.

The quasi-extinction risks for Population P2 highlight the difference between the WRD and Natural flows (Table 2) consistent with the risk curves (Figure 4d). For the Natural flow the risk of falling below 500 female adults was 0.015 and 0.123 for the WRD flow, and the risk of falling below 2,000 female adults was 0.227 for the Natural flow and 0.981 for the WRD flow.

Figure 5. Plots of modelled recruitment events paired with the proportion of available aquatic macrophyte for each population: population P3 (a) mean abundance and (b) macrophyte habitat proportion; population P2 (c) and (d); and population P1 (e) and (f). Orange – modelled WRD flow and Blue – modelled Natural flow: shaded areas are ± 1 standard deviation.

Figure 5. Plots of modelled recruitment events paired with the proportion of available aquatic macrophyte for each population: population P3 (a) mean abundance and (b) macrophyte habitat proportion; population P2 (c) and (d); and population P1 (e) and (f). Orange – modelled WRD flow and Blue – modelled Natural flow: shaded areas are ± 1 standard deviation.

The modelled effects of flow on macrophyte habitat resulted in episodic recruitment to 1 year olds (Figure 5). In populations P3 and P1, changes in macrophyte habitat and abundance of recruits were broadly consistent between the two flows (Figure 5a, b and 5e, f), whereas in comparison, in population P2, there were almost no available macrophytes (Figure 5d) and thus no recruitment under the WRD flow (Figure 5c).

4. Discussion

We developed a metapopulation model to conduct a population viability analysis (PVA) for Australian lungfish in the Burnett River which considered both viability and uncertainty to inform conservation management. The population model was constructed for this specific river for this long-lived, threatened species using the best available estimates of survival and fecundity rates, density dependence constraints and habitat loss from variable river flows impacting on the early life-history stages (eggs; larvae; and young-of-year fish) to ensure that the resultant model was both realistic and appropriate for this system [12].

The model considered the main identified threats to the species being the loss of habitat (naturally occurring and induced through water infrastructure) and the consequences of water infrastructure being the cause of other threats: altered hydrology and water regulation, spillway mortality and barriers to movement where for example fish experience mortality through exposure to excess salinity and predation by stranding in estuaries below tidal barrages [30]. The ability to model these influences on the long-term population viability of a long-lived, threatened species cannot be easily achieved through field-based monitoring techniques, and enables the assessment of various priority conservation measures, whilst being explicit about the uncertainties and minimising subjective judgements [77,78]. The lungfish PVA model is an important tool towards longer-term management of this species as it allows biologists and managers to test specific scenarios for the distinct sub-populations that may influence their viability as well as understanding consequences for the whole metapopulation.

4.1. The role of aquatic macrophytes as habitat for Australian lungfish

The Australian lungfish life-history traits align to that of periodic strategists, fish that are larger, long-lived and thereby reduce the risks of periods of unfavourable recruitment conditions [79]. Species with this strategy tend to have a lower degree of variability in abundance and high variability in recruitment, with the longevity of these species offsetting the period between strong recruitment events, which can be up to several decades [79]. As there is at least a 10-year lag from spawning to reproductive maturity in Australian lungfish, recruitment success may not be evident for decades using traditional field monitoring programs alone. The Burnett River lungfish PVA model has demonstrated that when the influence of flood flows on macrophyte development are included as inputs to spawning and recruitment success, that there are significant impacts to the resultant average adult population size over time (Figure 4, 5). This repeated loss of macrophytes due to flooding was considered the key reason for a decadal lack of recruits of Australian lungfish in the Mary River observed during the high flood period of the 1970s [54].

Whilst large dams like Paradise Dam on the Burnett River may reduce the peak discharge of flood flows and allow macrophytes downstream of the dam to persist during some flow events, this is compensated for by the significantly reduced macrophyte habitat availability and suitability throughout the inundated upstream reach (population P2) where up to 45 km of river channel now contains suboptimal habitat [32,36]. Additionally, dams can prevent the reestablishment of valuable macrophyte habitat through sequestration of seedbanks in deep sediments. Furthermore, large dams have the capacity to capture smaller flows disrupting natural hydrochory processes essential to aid recovery after destructive flood events that scour macrophytes [80]. This reduction in propagule dispersal required for recolonisation by macrophytes after large floods, could extend the period of reduced spawning potential for Australian lungfish, exacerbating the population impacts of habitat loss from inundation.

4.2. Impacts of water management

While inundation of aquatic macrophytes and the fluctuating water levels in storages have a direct effect on aquatic habitat for Australian lungfish in storages [66], interruption of movement patterns is also a key threat [30]. Direct monitoring has documented mortality of lungfish wounded as a result of impact with the stepped spillway of Paradise Dam following downstream displacement during high overtopping flows [36]. Translation of this spillway effect to population-levels may not be evident for many decades using traditional monitoring. Additional mortality of individuals occurs after flooding events due to being stranded downstream of the tidal barrage in the estuary as the fishway is ineffective and lungfish are unable to utilise it to return to freshwater environments. Both of these additional sources of mortality add to the loss of adults from the breeding population and are relevant to the conservation of Australian lungfish where adult survival was estimated to be very high, typical of long-lived species [81]. While the destructive impact of large floods has been included in the model, the influence of small stable flows during the breeding season necessary for spawning stimulus and egg development should be considered in future model development. Improved environmental flow management rules from storages in the Burnett River have, however, been shown to be associated with successful lungfish spawning downstream of storages [82]. These environmental management rules have been implemented throughout the core range of lungfish in the Burnett River and may offset some of the negative impacts of dams by providing a more natural low flow regime downstream.

4.3. Population viability and risk

The Australian lungfish population growth rate for the deterministic (matrix) model was calculated to be relatively low, $\lambda$ = 1.0587. The average number of female fish produced per female in her lifetime was relatively high $R_0$ = 8.13. Generation time (the average age of parents of the current cohort) was calculated to be very long G = 46.24. We do not consider these summary statistics unusual for a long-lived fish with a slow life-history, for example Bjørkvoll et al. [83] found that generally, marine species with a slow life-history had low reproduction and population growth, high survival rates and long generation times.

Our modelling of Australian lungfish in the Burnett River over a period of 118 years is the first attempt to understand population viability and extinction risk for this species. The results suggest that although the current hydrological threats may reduce the average adult population, it is likely that this results in only a modest increase in risk of the population being smaller for the total population when compared with the modelled Natural flow scenario (Figure 4). However, it is likely that there is a considerable increase in risk for the population within the Paradise Dam impoundment with a relatively high likelihood of the population being small under the WRD flow scenario. This highlights the benefits of developing a metapopulation model to be able to examine the impacts of spatially variable habitat alterations and effects on individual populations within a large river system with significant regulation. Conservation of Australian lungfish needs to consider the multiple genetic stocks across the current range and hence the importance of managing threats unique to each [20]. The original conservation advice listed Australian lungfish as ‘Vulnerable’ highlighted that there was ‘sufficient evidence to suspect that the adult breeding population will undergo a substantial decline over the next three generations’ [30]. This listing preceded significant research into age estimation, population estimations from tag-recapture and genetic differentiation. The generation time for Australian Lungfish in the Burnett River was calculated to be very long (G = 46.24), with three generations equating to 138 years. This is similar to the modelled simulation of 118 years where the model suggested that there may be a 28% reduction in the abundance of adult female lungfish, assuming water infrastructure has been in place for that whole period. Whilst the outcome of this scenario is of concern and aligns with the original conservation listing, there is no evidence to suggest these pressures will cause the extinction of the species. The Burnett River PVA model now allows us to test a range of scenarios to understand which management actions may be taken to reduce the risk to the metapopulation [10].

Population viability models are essential tools for prioritising management actions for long-lived species where generation time often exceeds even the longest long-term monitoring programs in well researched species [84]. For example, population modelling of the US Atlantic loggerhead sea turtle (Caretta caretta), a species with a ~70-year lifespan [85], highlighted that mortality of adult turtles by drowning in prawn trawls was a major threat that could be mitigated through implementation of turtle excluder devices [81]. Even though these devices became mandatory in 2002, due to the long generation time for marine turtles, it has taken nearly two decades for the numbers of nesting turtles within the local area to improve [86].

4.4. Future uses of the model

The existing Burnett River PVA model for Australian lungfish has highlighted the trajectory of the species in the river system with the existing modelled threats. Due to construction issues with Paradise Dam, the spillway height has been reduced and a new dam is currently being investigated [87]. In addition, the Burnett Basin water plan is being reviewed which means that hydrological modelling for the river system will be updated with new stream flow sequences that include the new dam design and an additional 14 years of streamflow data [40]. Updating the flow sequences in the PVA model will improve the confidence in understanding the longer-term impact of spillway mortality, macrophyte losses, spawning successes and losses to the estuary as all of these effects have been measured for the past decade [36]. This highlights the importance of updating models and modelling as new data, knowledge and management scenarios arise. Whilst the PVA model can be used as purely to understand extinction risks, its real power is to compare scenarios with altered parameters to evaluate outcomes of future conservation actions. Such actions could include reducing the lag period between flood scouring and macrophyte recovery downstream of large dams by actively rehabilitating lost aquatic macrophytes, improved fishway passage effectiveness for lungfish on modelled barriers, or reduced mortality rates of altered spillway design. Population models can also be used to understand where the vulnerability of the species (especially for long-lived species), lies in relation to the life cycle, whether that be the early life-history, breeding, adult mortality or at all stages [81,88].

Many of the knowledge gaps originally identified within the Draft National Recovery Plan for Australian lungfish have been advanced, including understanding the age structure of the population [36,54]. Such information has provided an essential basis for this modelling and the subsequent conservation management of this threatened species. The populations in separate drainages need to be managed as genetically distinct units, as they are experiencing differing threatening processes. The increased research on the species in the past decade has enabled realistic and precise construction of a species population model to represent geographically distinct populations [12]. The metapopulation model for the Burnett River Australian lungfish, comparing the two flow scenarios, identified a significant increase in risk of the population being small for one of the populations most affected by habitat modifications and a modest increase in risk for the total population within the river system. The model has the flexibility to test future conservation strategies to manage impacts, explore options to promote sustainable population levels and in the longer-term, may be used to help remove Australian lungfish’s threatened species status.

5. Conclusions

We developed a metapopulation model for a very long-lived threatened species using interventions that target different elements of the species’ life-history that drives subsequent population dynamics. The study of long-lived species may be restricted by a lack of high-quality long -term data, and particularly in the case of threatened species, restricting our capacity to predict systemic decline or recovery in response to threat amelioration. As such, it is imperative that tools such as this metapopulation model are used to guide and develop future management strategies, particularly for species of high conservation value, even where data may be limited. Alongside targeted research and monitoring that address knowledge gaps, metapopulation models can inform species conservation efforts. The PVA model presented here can be updated as knowledge on life-history or key threats become better understood to re-test alternate management scenarios into the future. Moreover, the approach described here could be adapted to inform management of other long-lived threatened species.