1. Factors that determine C4 grass distributions have been well documented, with evidence in the literature for C₄ photosynthetic subtypes displaying varying levels of drought susceptibility. However, the interactions between C4 photosynthetic subtype and phylogeny add complexity and are relatively under-studied.  
2. We use species distribution modelling to determine the influence of rainfall on distribution patterns of representative C₄ grass families and subtypes. Select C₄ grass species, representing different photosynthetic subtypes (NADP-Me and NAD-Me) and lineages (Panicoideae and Aristidoideae), were then subjected to a progressive 58-day drought period and recovery phase, to explore drought responses through leaf water relations, gas exchange and chlorophyll fluorescence.
3. We show Panicoideae NADP-Me species to be more susceptible to drought than both Panicoideae NAD-Me and Aristidoideae NADP-Me species due to apparent greater metabolic impairment. The differences between groups were related to how rapidly photosynthesis declines with exposure to drought and the rate of recovery post-drought, rather than the maximum extent of photosynthetic decline. The mechanisms for the relative maintenance of plant water status differed between the Panicoideae NAD-Me species, which utilised greater stomatal control, and the Aristidoideae NADP-Me species, which maintained water uptake through osmotic adjustment. 
4. Synthesis: We show here that drought susceptibility differs both phylogenetically and according to photosynthetic subtype, but that the role of phylogeny may outweigh physiological control. This research adds novel insight into the physiological differences behind observed rainfall-related differences in C₄ grass distribution patterns.

Materials and methods

Species distribution modeling

MaxEnt Species Distribution Models (SDMs) were developed using the maximum entropy algorithm implemented in MaxEnt v.3.3.3k to display the current SDMs for an array of C₄ grass species from the Panicoideae and Aristidoideae families. Maximum Entropy Models (MaxEnt) use machine-learning principles to predict the presence/absence of species across geographic space based on the relationship between site-specific environmental conditions and species occurrence records (Elith et al., 2011). From this, MaxEnt models can predict a measure of habitat suitability for a particular species in a particular location. Habitat suitability is a factor of both the density of individuals of a species and the likelihood of the species being present (Elith et al., 2011). Although there are limitations to the chosen method (e.g., using solely presence data as opposed to presence-absence records; Phillips et al., 2009), MaxEnt SDMs provide a good base from which to guide our physiological hypotheses. SDMs were run for the study species as well as for multiple southern African Panicoideae and Aristidoideae species to confirm that patterns held throughout the groups (see species listed in Table S1). Models were separated into (a) ‘study species models’ which includes three species of Panicoideae NADP-Me, two species of Panicoideae NAD-Me, and three species of Aristidoideae NADP-Me (see Table 2) and (b) ‘multi-species models’ which includes multiple Panicoideae (38 species) and Aristiodoideae (6 species) species that occur in southern Africa. Species occurrence records were obtained from the Global Biodiversity Information Facility (GBIF, httpp://www.gbif.org) accessed on 20 November 2021. Environmental spatial data was obtained from the WorldClim—Global Climate Data database (http://www.worldclim.org; Fick & Hijmans, 2017), providing 19 environmental layers for MaxEnt modeling (see Table S2). MaxEnt models were trained with the whole extent of southern Africa as the species maximum potential ranges. The MaxEnt models produced a spatial ‘habitat suitability’ (H.S.) based on a scale from 0 to 1 for southern Africa (including South Africa, Lesotho, eSwatini, Botswana, Namibia, Zimbabwe, Mozambique, Angola, Zambia, Malawi, Tanzania, and the Democratic Republic of Congo) as well as probability curves relating to each environmental variable. Of most interest for this study was the correlation between precipitation and probability of occurrence, comparing areas of peak occurrence along a rainfall gradient between plant groups. Area under the curve (AUC) values were used to assess model performance using 10 000 random pseudoabsences during model evaluation. AUC values for all models were > 0.85, and mean AUC was 0.94 (Table S1).

Plant collection, growth conditions and experimental set-up

The species classified as ‘study species’ in the models were experimentally grown and studied in a drought experiment (see Table S3 for more detailed information of each species). Tristachya leucothrix, Heteropogon contortus, Aristida diffusa and Aristida congesta were collected from the Makhanda area (Eastern Cape, South Africa), while Aristida junciformis was collected at Port Alfred (Eastern Cape, South Africa). Whole plants were dug up in the field, trimmed and potted such that each pot represented an individual plant. Panicum coloratum, P. stapfianum and Alloteropsis semialata were grown from existing potted plants that were trimmed and re-potted. Alloteropsis semialata was collected at Middelburg, Gauteng Province, South Africa, (S25°50’, E29°24’), and P. coloratum was grown from seed sourced by Taylor et al. (2010). Table 2 gives a summary of the species used in the experiment. All plants were potted in 10 L pots containing 6.7 kg of a homogenous soil mixture made from locally obtained top-soil, representative of the soil the grasses grow in naturally, and kept in a clear polythene tunnel at the Department of Botany, Rhodes University. Average min/max temperatures (± SD) of the tunnel for the duration of the experiment were 17.6 ± 1.2 and 34.1 ± 4.6°C respectively and the average tunnel temperature was 25.1 ± 8.8°C. Diurnal photosynthetically active radiation (PAR) in the tunnel ranged from 64 to 1014 µmol s^-1 m^-2, with midday values of 1012 ± 2 µmol s^-1 m^-2. Plants were well watered (using field capacity of the soil as a guide) and hydroponic fertilizer (Chemicult - 1g L^-1) was added twice in the month leading up to the experiments. Six treatment and six control replicates (except P. coloratum and A. diffusa which had five replicates due to mortality) of each species were used in all the experiments. 

Progressive drought was imposed by starting experiments with potted plants watered to field capacity (±20% SWC), and then allowing them to decrease soil water content (SWC) by ±0.3% each day over the subsequent 58 days (see Fig. 1). This daily reduction is representative of natural soil drying (Ripley et al., 2010).  On day 58 plants were re-watered and maintained at field capacity over the remaining 11 days (recovery phase). During the dry-down phase of the experiment, potted plants were weighed every second day, and supplementary water added where necessary to ensure that all plants dehydrated at similar rates. Field capacity of the soil was determined by soaking pots in water for 24 hours and then allowing the soil to drain to constant mass under gravity. During this period, the evaporation from the soil surface was minimized by covering the pots with plastic lids. To estimate SWC it was necessary to determine the dry weight of the soil added to each pot and to estimate the weight of the plants. Soil dry weights were determined by oven-drying soil at 70ºC for 72 hours and representative plants' weights were determined by harvesting a subset of plants from each species. Evaporation from the soil during the experiment was minimized by adding 1 kg of fine stone (< 1 cm diameter) to the soil surface. Hence as plant, soil, pot and stone weights were accounted for, the % SWC for the potted plants could be calculated as  SWC(%) =  [(soil wet mass - soil dry mass)/soil dry mass] x 100.

Leaf gas exchange, chlorophyll fluorescence and plant water relations

Gas exchange, chlorophyll fluorescence, and leaf water relations (Ψleaf and RLWC) were measured on various occasions during the dehydration and re-watering phase of the experiment (Fig. 1). Instantaneous measures of net CO₂ assimilation rates (A), stomatal conductance (G_st), intrinsic water-use efficiency (A/ G_st) and intercellular CO₂ concentration (Ci) were conducted on the youngest fully expanded leaf (first down from the apical bud) of the control and treatment plants. These parameters were measured on the days indicated in Fig. 1, except for the Aristidoideae species which were not measured on day 10. Measurements were made using a Li-6400-40 LCF photosynthesis system (Li-Cor Inc., Lincoln, NE, USA) between 10:30 am and 3:30 pm under laboratory conditions. Plants were acclimated under a sodium vapour light at a photosynthetic photon flux density (PPFD) similar to that used in the leaf chamber. Cuvette conditions were maintained as follows: incoming (reference) ambient CO₂ concentration (Ca) was supplied at 400 μmol mol⁻¹, a PPFD of 1200 μmol m^-2 s^-1 was supplied by a blue-red LED light source, leaf temperature was set at 29°C and vapour pressure deficits (VPD) ranged between 1 – 2.5 kPa. Leaf areas were measured and entered manually, and gas exchange parameters were calculated according to the equations of von Caemmerer and Farquhar (1981). To ensure that measurements were conducted at near saturating light intensities, photosynthetic response to incident light intensity was measured on control plants according to Long and Bernacchi (2003). 

Chlorophyll fluorescence measurements were made immediately following each gas exchange measurement as not to disrupt the steady state photosynthesis. Leaves were acclimated until steady-state fluorescence (F_s) was achieved. A multiphase flash (MPF) protocol was used to ensure maximum reduction of QA. The following MPF settings were used: 30% ramp, 250ms for phase 1 and 3 and 500ms for phase 2. The light intensity required to ensure QA reduction was experimentally determined. Chlorophyll fluorescence parameters measured are defined (Baker, 2008) and where necessary their calculations and units are shown. PSII maximum efficiency, Fv’/Fm’ = (Fm’ – Fo’) / Fm’. At a given photosynthetic photon flux density (PPFD), this estimates the maximum PSII photochemistry (efficiency of oxidised (QA) PSII reaction centers). Fm’ is the maximal fluorescence during the saturating light phase (PPFD > 7000 μmol m^-2 s^-1) (QA maximally reduced) and Fo’ is the minimal fluorescence of a briefly darkened (6 seconds at 740nm), light-adapted leaf (QA maximally oxidised). PSII operating efficiency, ΦPSII = (Fm’ – Fs) / Fm’. At a given PPFD, this estimates the efficiency at which light absorbed by PSII is used for QA reduction, (steady state photosynthesis). Photochemical quenching, qP = (Fm’ – Fs) / (Fm’ – Fo’). At a given PPFD, this estimates the PSII reaction centers (QA) that are oxidised. This includes photosynthesis and photorespiration. Electron transport rate, ETR = ΦPSII x f x I x αleaf (μmol electrons m⁻² s⁻¹).

Midday leaf water potentials, relative leaf water contents and osmotic adjustment

The leaves used for gas exchange measurements were either excised at midday on the same day as the gas exchange measurements or on the following day. The excised leaves were immediately weighed, and the leaf water potential (Ψleaf) was measured using a Schőlander pressure chamber. Following Ψleaf measurements, the leaves were placed upright in a glass vial which contained enough water to cover the first 10mm of the excised end of the leaf. The leaves were left in the dark overnight to regain full turgor pressure. The following morning the leaves were blotted, weighed, and then placed in a drier at 70 ºC for 48 hours, after which they were weighed again. This method allowed the Ψleaf and relative leaf water content (RLWC) to be obtained for the same leaf. Trial experiments were conducted to determine if the measurement of Ψleaf with a pressure chamber affected the rehydration of leaves and it was found to have no significant effect. . Ψleaf could not accurately be measured at day 56 (± 3.5% SWC) because of the extreme leaf dehydration. Ghannoum et al. (2003) showed that the relationship of Ψleaf to RLWC was mostly linear. Models were fitted to the mean Ψ_leaf and corresponding RLWC data for each species (control and treatment plants) during the drought and recovery phase. All the species showed a linear relationship of Ψleaf to RLWC, thus a straight-line function (y = mx – c) was used to describe Ψleaf at day 56 (see Fig S3).

Pressure volume (PV) curves were constructed by determining the relationship between RLWC and Ψleaf (see Fig. S3). This was done by sequentially dehydrating leaves and determining Ψleaf and RLWC at regular intervals. Initially well-watered potted rooted plants were bagged overnight to ensure the leaves reached full turgor potential. For 11 replicate plants, the fully expanded second or third leaf produced after the cotyledon was excised, weighed, and the corresponding Ψleaf was obtained by using the Schőlander pressure chamber. Subsequent leaves were allowed to slowly dehydrate in a humidified bell jar and Ψleaf and RLWC were measured at repeated intervals. PV curves constructed from control leaves were used to determine individual components of leaf water potential (Ψleaf = ΨP + Ψπ; see Fig. S4). The reciprocal of the Ψleaf was plotted against RLWC and a model was fitted to the data using the equations of Schulte and Hinckley (1985). The TLP (ΨP= 0 MPa) was defined as the RLWC at which Ψleaf equaled Ψπ (osmotic potential). Ψπ at 100% RLWC was calculated as the y-intercept of the straight line. To determine the osmotic adjustment (OA) of treated plants, their Ψleaf and RLWC were measured at various intervals during the dry-down experiment, once RLWC had declined sufficiently to ensure that ΨP= 0 and that changes in 1/ Ψleaf were solely dependent on 1/Ψπ. This response of 1/ Ψleaf to RLWC was fitted with a straight line and 1/Ψπ at 100% RLWC (y-intercept) was calculated. OA was calculated from the change in Ψπ with drought stress. This was expressed in absolute terms for well-watered individuals or relative to the values for well-watered controls for individuals exposed to drought. Osmotic adjustment was defined as the difference between Ψπ of control and drought-treated plants at 100% RLWC. As different species showed inherently different Ψπ, OA was expressed as a percentage of the control value (relative OA) to allow the comparison of OA between species. 

Statistics

The data obtained were analysed using the R language and platform (R Core Team 2020). To determine the influence of the treatments (both plant type and days of exposure to drought) on plant ecophysiology and water relations, linear mixed-effects models were fitted in R using the ‘nlme’ package (Pinheiro et al., 2021).  Species were treated as a random factor in the models. The use of planting treatment as a random factor in the statistical models was initially included to avoid planting differences between species influencing the results but was removed once it was determined that it had no effect on the statical results. ANOVAs were conducted on the mixed-effect model outputs to get appropriate pairwise and global model statistics. Statistics presented in the results show the interaction between plant type and days of exposure to drought (i.e., differences between plant types in responses to drought over time). Comparisons including watered (control) plants are presented in Fig S2. Model evaluation was conducted for all MaxEnt models run, ensuring the Average Area Under the Curve (AUC) values were greater than 0.7. To determine the difference in habitat suitability relative to mean annual precipitation between plant types, Chi² tests were conducted on Generalised Additive Model (GAM) distributions using the ‘gam’ package in R (Hastie, 2020).

The data obtained were analysed using the R language and platform (R Core Team 2020).