Subsoil Microbial Diversity and Stability in Rotational Cotton Systems

Microbial diversity has been well documented for the top 0–0.30 m of agricultural soils. However, spatio-temporal research into subsoil microbial diversity and the effects of agricultural management remains limited. Soil type may influence subsoil microbial diversity, particularly Vertosols. These soils lack distinct horizons and are known to facilitate the downward movement of organic matter, potentially supporting subsoil microbiota, removed from the crop root system (i.e., bulk soils). Our research used the MiSeq Illumina Platform to investigate microbial diversity down the profile of an agricultural Australian Vertosol to 1.0 m in bulk soils, as influenced by crop system (continuous cotton and cotton–maize) and sample time (preand in-crop samples collected over two seasons). Overall, both alpha(Chao1, Gini–Simpson Diversity and Evenness indices) and beta-diversity (nMDS and Sørensen’s Index of Similarity) metrics indicated that both bacterial (16S) diversity and fungal (ITS) diversity decreased with increasing soil depth. The addition of a maize rotation did not significantly influence alpha-diversity metrics until 0.70–1.0 m depth in the soil, where bacterial diversity was significantly higher in this system, with beta-diversity measures indicating this is likely due to root system differences influencing dissolved organic carbon. Sample time did not significantly affect bacterial or fungal diversity over the two seasons, regardless of the crop type and status (i.e., crop in ground and post crop). The relatively stable subsoil fungal and overall microbial diversity in bulk soils over two crop seasons suggests that these microbiota have developed a tolerance to prolonged agricultural management.


Introduction
Soil microbiota mediate a range of functions within the soil environment, including the decomposition of organic matter (OM) and subsequent release and cycling of essential nutrients-carbon (C), nitrogen (N) and phosphorus (P)-into the terrestrial ecosystem [1][2][3][4][5][6]. These mechanisms are not only vital to a healthy soil, but are key in managing and maintaining sustainable agricultural practice [7][8][9]. As such, an understanding of the influence of disturbance, such as agricultural management and climate change on soil microbiota, is of utmost importance.
In cropping systems, microbiota are subject to disturbance from agricultural practice including minimum tillage, irrigation and crop rotations-all of which impose physical changes in the soil structure and nutrient availability [10]. These disturbances to the soil environment may have ramifications for microbial communities, exhibiting sensitivity (the extent to which a microbial community is altered in response to a disturbance), resistance (the extent to which a microbial community can resist change (DOC)/1.10% (SOC) at 0-0.15 m to 48 mg/kg (DOC)/0.7% (SOC) at 0.3-0.5 m and 32 mg/kg (DOC)/0.58% (SOC) at 0.7-1 m depth. Other soil characteristics (i.e., P, K, and Ca) for this site have been published and can be found in Polain, Guppy, Knox, Lisle, Wilson, Osanai and Siebers [26], Nachimuthu, Watkins, Hulugalle, Weaver, Finlay and McCorkell [29] and Hulugalle, Nachimuthu, Kirkby, Lonergan, Heimoana, Watkins and Finlay [27]. The field site is part of a long-term (approximately 30 years) study, where historical tillage and rotations constituted the main plot treatments, with maize being added as a rotation in 2011 as subplots. This site is managed by the NSW Department of Primary Industries (NSW DPI) and undergoes minimum tillage in the top 0.1 m of the soil profile, post picking, as part of an insect management regime to eliminate Helicoverpa spp. (Cotton bollworm or Native budworm) larvae in a process known as 'pupae busting' [44]. Crops are irrigated at a rate of 1 ML ha −1 (equivalent to 100 mm rain) at 7-20 day intervals with approximately 4-8 irrigations per season depending on the amount of seasonal rainfall. Both systems received 260 kg N ha −1 per season, with no P fertiliser applied due to critical Colwell P levels being <6 mg kg −1 for cotton [45]. Samples were collected at the midpoint between the head ditch and tail drainage ends (the site map can be viewed in Osanai, Knox, Nachimuthu and Wilson [25]) in continuous cotton (CC) and cotton-maize (CM) plots within the same field [27,29] over two seasons and at different crop status, these being post pupae busting and prior to planting in October 2015 and 2016 (pre-crop) and at flowering status (in-crop) January 2016 (CC-cotton; CM-maize) and 2017 (CC-cotton; CM-cotton). During season 1 (2015/2016), the minimum and maximum soil temperatures were 6.7 and 36.5 • C, respectively, with a total rainfall of 537 mm and irrigation volumes of 720 mm (CC) and 564 mm (CM) [28,46]. During season 2 (2016/2017), the minimum and maximum soil temperatures were 8.2 and 40.2 • C, respectively, with a total rainfall of 318 mm and irrigation volumes of 1008 mm for both systems [28,46].
Sampling was established as a split-plot design, based on the main plot treatments, where there were four replicate plots each for CC and CM, with one core extracted per replicate plot at each time of sampling in the plant line between cotton plants (in-crop), or from the reformed planting hill (pre-crop), which resulted in a total of 32 cores extracted over the two seasons. These soil cores were extracted using a powered soil corer, with a 42 mm inner diameter, without the use of lubricants [47] to the depth of 1 m and placed into polyvinyl chloride (PVC) half-pipes, wrapped with plastic wrap to store and support the core. Cores were transported to the laboratory in cooled, insulated containers, where subsamples were removed in triplicate from three depths (0-0.15, 0.3-0.5 and 0.7-1 m), avoiding any root material and stored at −20 • C until required.
DNA was extracted from each subsample using PowerSoil ® DNA Isolation Kits (MO BIO Laboratories Inc. Carlsbad, CA, USA) according to the manufacturer's instructions. Briefly, 0.25 g of soil was weighed and subjected to a series of steps to promote cell lysis, removal of inhibitors (i.e., non-DNA material), and washing of DNA, before being eluted in 65 µL elution buffer. The triplicate samples were pooled, with DNA quantity (1 ng/µL to 50 ng/µL required) and purity (A260/A280 ratio range 1.6 to 1.9) assessed using a NanoDrop ™ 8000 Spectrophotometer (ThermoFisher Scientific, Waltham, MA, USA).
Three technical replicates of the pooled soil DNA soils were sent to the Australian Genome Research Facility (AGRF) in Brisbane (Queensland) for Diversity Profiling using the MiSeq Illumina Platform. This service identifies the relative proportion of organisms within mixed microbial communities [48]. AGRF-validated primers 341F and 806R were used to target 16S variable regions V3-V4 for bacteria and archaea, whilst primers 1F and 2R were used to target ITS regions for fungi (Table 1).
Bioinformatics methods were completed by the AGRF as part of quality control [48]. Briefly, these methods involved the initial assembly of the paired ends of the forward and reverse reads using PEAR (version 0.9.5), followed by primer trimming [49]. The trimmed sequences were processed using Quantitative Insights into Microbial Ecology (QIIME 1.8) [50], USEARCH (version 8.0.1623) [51,52] and UPARSE [53] software. Sequences were quality filtered, with duplicate sequences, singletons and unique sequences removed and data sorted by abundance [48]. Sequences were clustered and chimera filtered using "rdp_gold" for 16S [52,54,55] and "Unite" for ITS [56] databases as references.
The number of reads in each OTU were mapped back to OTUs with a minimum identity of 97%. QIIME was used to assign taxonomy using the Greengenes databases for 16S sequences [54] and the Unite database for ITS sequences [56]. Prior to the construction of mathematics functions (indices) to determine soil microbial diversity, the raw and quality control reads were statistically compared using IBM SPSS Statistics 25 program for statistical interrogation. Repeated-measures Analysis of Variance (ANOVA) was used to assess significant differences (p ≤ 0.05) between the experimental factors of time (season and crop status), depth, system, time by depth, time by system and depth by system. Table 1. Forward and reverse primers used to target 16S and ITS regions [48].

Target
Primer Sequence Alpha-diversity (α-diversity), the diversity within a sample (i.e., the diversity within a sample obtained from a particular depth), was assessed using the Chao1 (S Chao1 ) estimator for species richness [57,58], Gini-Simpson Diversity (D GS ) and Gini-Simpson Evenness (E GS ) indices (Table 2) [59][60][61]. Chao1 (Equation (1)) is a non-parametric method to estimate species richness, with the premise that rarer species infer more information about the number of species missing from a population [57,58].
where S is the number of species, F 1 the number of singletons and F 2 the number of doubletons in an assemblage [62]. The Gini-Simpson Diversity Index (Equation (2)) takes into consideration species richness and abundance, measuring the probability that two randomly selected individuals belong to two different species [58,60].
where S is the total number of species in a community and p i is the proportion of individuals of a particular species (OTU i ), divided by the total number of species found (S) [60,63]. Ranging from 0 (low diversity) to 1 (high diversity), the D GS index favours abundant species, highlighting the dominant species in a microbial community [60]. The Gini-Simpson Evenness Index (Equation (3)) assesses the abundance (evenness) of the dominant species within an assemblage, with an increased evenness value indicating that the assemblage is more diverse [64].
where D GS is the Gini-Simpson Diversity value determined (Equation (2)) and S is the number of species observed [59]. The E GS value is considered meaningful when evaluating the microbial community's ability to resist or recover from disturbances, where a value closer to 0 indicates a lower species abundance and, therefore, lower level of diversity [64]. These aforementioned indices were calculated and constructed in Microsoft Excel, then imported into IBM SPSS Statistics 25 program for statistical interrogation. Repeated-measures Analysis of Variance (ANOVA) was used to assess significant differences (p ≤ 0.05) between the experimental factors. Where data were not normally distributed, log 10 (Log 10 ) or square root ( √ ) transformations were performed prior to ANOVA.
Beta-diversity (β-diversity), the diversity observed between two samples (i.e., the diversity between samples recovered from CC and CM at 0-0.15 m depth) was determined using PRIMER-e v7 (Quest Research Limited, Albany, New Zealand) software, where system by depth absolute abundance OTU data at the phyla level were used for non-metric multidimensional scaling (nMDS) based on Bray-Curtis dissimilarity [65]. Permutational Multivariate Analysis of Variance (PERMANOVA) was used in the analysis of the Bray-Curtis distance data for depth, depth by system and time. Assessed and available soil environmental data (Table S1) were obtained and used with the β-diversity data to perform distance-based linear modelling (DISTLM) sequential tests, which were used to plot distance-based redundancy analysis (dbRDA), in an effort to identify some of the main drivers in our system ( Figures S1 and S2). The results of these data interrogations have been used to support theories for system differences or similarities where relevant.
Finally, Sørensen's Index of Similarity (I S ) was used to observe microbial community similarities at the Phylum level [66,67]. Where With C being the number of species two samples have in common, S 1 is the total number of species in sample 1 and S 2 is the total number of species in sample 2 [66].

Quality Control and Microbial OTUs
The number of OTU reads significantly differed (p < 0.001) between the raw and quality control data, with a difference of 14,286 ± 358 for 16S and 39,879 ± 1302 for ITS data. Statistical analysis of the data set type (raw versus quality control) against depth (16S p = 0.965; ITS p = 0.320), crop system (16S p = 0.845; ITS p = 0.664) and sampling point (16S p = 0.942; ITS p = 743), revealed no significant differences. Statistical analysis of the data sets + depth + crop system (16S p = 1.000; ITS p = 0.742) and data sets + depth + sample point (16S p = 0.999; ITS p = 0.829) also revealed no significant differences.

Bacteria (16S)
The average species richness (S Chao1 ) estimate for bacteria (16S) was significantly influenced  (Table 2). Average bacterial species richness estimates were not significantly affected by the remaining factors of crop system (p = 0.638) and sample time (p = 0.093). No significant difference was observed for the average species richness for bacteria at depth as influenced by crop system (p = 0.727) or for species richness by depth as influenced by crop system and sample time (p = 0.114) ( Table 2). However, there was a significant difference Depth significantly (p < 0.001) influenced the diversity of bacterial (16S) communities, with the average bacterial Gini-Simpson Diversity Index (D GS ) decreasing from 1.00 ± 2 × 10 −4 (0-0.15 m) to 0.99 ± 1 × 10 −3 (0.70-1.0 m) ( Table 2). Crop system significantly (p = 0.005) influenced D GS , with CM (0.99 ± 5 × 10 −4 ) having increased diversity, in comparison to CC (0.99 ± 9 × 10 −4 ). Crop system also significantly (p = 0.011) influenced bacterial diversity at depth (0.70-1.0 m), with a greater decrease in bacterial diversity in CC (average D GS = 0.99 ± 2 × 10 −3 ) than CM (average D GS = 0.99 ± 9 × 10 −4 ) ( Table 2). No significant effect (p = 0.247) was detected for depth, as influenced by sample time and crop system on bacterial diversity, with D GS values averaging 0.99 ± 5 × 10 −4 . However, bacterial community diversity decreased with depth at each sample time over the two seasons, independent of crop system, except for the notable decreases in diversity in CC at 0.70-1.0 m, pre-crop for both seasons 1 and 2.
Depth did not significantly influence the evenness (E GS ) of bacterial communities (p = 0.391) down the soil profile (Table 2). However, the average E GS value for bacterial communities increased between depths 0.30-0.50 m (4 × 10 −4 ± 2 × 10 −5 ) and 0.70-1.0 m (5 × 10 −4 ± 2 × 10 −5 ). Crop system did not significantly influence the evenness of bacterial (p = 0.388) communities down the soil profile, however, with the average bacterial E GS value for both CC and CM being 10 −4 ± 6 × 10 −5 ( Table 2). The CM system trended in higher evenness of microbial communities at 0-0.15 and 0.30-0.50 m, whilst evenness increased in the CC system at 0.70-1.0 m. Sample time and crop system did not significantly influence the evenness of bacterial (p = 0.349) species down the soil profile (Table 2), with an average E GS of 4 × 10 −4 ± 4 × 10 −5 . Whilst the evenness of bacterial communities was relatively consistent across systems, there was an increase in evenness during season 1 pre-crop at 0-0.15 m under CM, with an average of 1 × 10 −3 ± 9 × 10 −4 , which was >2-fold greater than the remaining average values.
There was no apparent influence of depth on fungal (ITS) community diversity (p = 0.637), with average fungal D GS differing by 0.01 from 0-0.15 to 0.70-1.0 m, with an increase of 0.04 at 0.30-0.50 m (Table 3). Fungal diversity down the soil profile was not significantly (p = 0.938) influenced by crop system. However, the average D GS was 0.03 ± 0.02 lower under the CM system (Table 3). No significant effect (p = 0.486) was detected for depth, as influenced by sample time and system on fungal diversity, with D GS values averaging 0.90 ± 1 × 10 −2 . Fungal community diversity did not demonstrate a clear change with soil depth over the two growing seasons although there was a notable, but not significant difference between crop system diversity at 0-0.15 m pre-crop, season 1, with a decrease in diversity in the CM system.
Depth did not significantly influence the evenness of fungal communities (p = 0.715) down the soil profile (Table 3), with average E GS increasing from 0-0.15 and 0.30-0.50 m (3 × 10 −3 ± 7 × 10 −4 ) were lower in comparison to 0.70-1.0 m (8 × 10 −3 ± 7 × 10 −4 ). Crop system did not significantly influence the evenness of fungal (p = 0.556) communities down the soil profile, with the average E GS values differing by less than 0.001 between each system, at each depth ( Table 3). The CM system trended in higher evenness of fungal communities at 0-0.15 and 0.70-1.0 m, whilst evenness was increased in the CC system at 0.30-0.50 and 0.70-1.0 m. Sample time and crop system did not significantly influence the evenness fungal (p = 0.744) species down the soil profile (Table 3), with an average E GS of 6 × 10 −3 ± 4 × 10 −4 . Overall, the evenness of fungal communities fluctuated over the two growing seasons, with the E GS index particularly variable at 0.70-1.0 m.

Dissimilarities between Microbial Communities
PERMANOVA analysis of absolute abundance OTUs for bacterial (16S) found a significant difference between samples collected at different depths (p = 0.001) and samples collected at depth under the different crop rotations (p = 0.001). However, there were no significant differences between the diversity of samples collected from different crop systems (p = 0.139) or sample time (p = 0.087). Non-metric multidimensional scaling (nMDS) was used to visualise the differences in bacterial communities down the soil profile and between the two systems. Clear separation of clusters with depth and to a lesser extent, separation by system, was evident particularly between CC and CM at 0-0.15 m and 0.70-1.0 m (Figure 1). Soil Syst. 2020, 4, x 9 of 18

Dissimilarities between Microbial Communities
PERMANOVA analysis of absolute abundance OTUs for bacterial (16S) found a significant difference between samples collected at different depths (p = 0.001) and samples collected at depth under the different crop rotations (p = 0.001). However, there were no significant differences between the diversity of samples collected from different crop systems (p = 0.139) or sample time (p = 0.087). Non-metric multidimensional scaling (nMDS) was used to visualise the differences in bacterial communities down the soil profile and between the two systems. Clear separation of clusters with depth and to a lesser extent, separation by system, was evident particularly between CC and CM at 0-0.15 m and 0.70-1.0 m (Figure 1).   Dissimilarities between fungal (ITS) communities down the soil profile, as influenced by continuous cotton (CC) and cotton-maize (CM) crop systems, using nonmetric multidimensional scaling (nMDS) of absolute abundance OTUs based on the square root transformation of Bray-Curtis dissimilarity, with a stress score of 0.14.

Similarities between Microbial Communities
Similarities between bacteria and fungi community profiles significantly decreased (p < 0.001) between the depths of 0-0.15 cm and 0.70-1.0 m, with Sørensen's Index of Similarity (IS) values of 0.77 and 0.53 for bacterial and fungal communities, respectively (Figure 3). At each depth, Acidobacteria (14-19%), Actinobacteria (29-49%) and Proteobacteria (14-27%) were the prevalent bacterial phyla (Figure 3a). Both Actinobacteria and Nitrospira increased with each depth increment sampled, being more prevalent in the 0.30-0.50 and 0.70-1.0 m samples than they were in the 0-0.15 m samples. Ascomycota (47-61%) and Basidiomycota (19-35%) were amongst the most prevalent phyla for fungi (Figure 3b). Under the CC system, the fungal groups were relatively constant, with the exception of Zygomycota becoming less abundant at 0.30-0.50 cm, before returning to a higher abundance at 0.70-1.0 m. In contrast, under the CM system, the relative abundance of Agaricomycetes, belonging to the Basidiomycota phylum, increased with depth, whilst Ascomycota, namely the Sordariomycetes, decreased in relative abundance. Whilst crop system did not significantly influence bacteria and fungi profiles with depth, IS values between CC and CM systems decreased by 0.08 (bacteria) and 0. 16

Similarities between Microbial Communities
Similarities between bacteria and fungi community profiles significantly decreased (p < 0.001) between the depths of 0-0.15 cm and 0.70-1.0 m, with Sørensen's Index of Similarity (I S ) values of 0.77 and 0.53 for bacterial and fungal communities, respectively (Figure 3). At each depth, Acidobacteria (14-19%), Actinobacteria (29-49%) and Proteobacteria (14-27%) were the prevalent bacterial phyla (Figure 3a). Both Actinobacteria and Nitrospira increased with each depth increment sampled, being more prevalent in the 0.30-0.50 and 0.70-1.0 m samples than they were in the 0-0.15 m samples. Ascomycota (47-61%) and Basidiomycota (19-35%) were amongst the most prevalent phyla for fungi (Figure 3b). Under the CC system, the fungal groups were relatively constant, with the exception of Zygomycota becoming less abundant at 0.30-0.50 cm, before returning to a higher abundance at 0.70-1.0 m. In contrast, under the CM system, the relative abundance of Agaricomycetes, belonging to the Basidiomycota phylum, increased with depth, whilst Ascomycota, namely the Sordariomycetes, decreased in relative abundance. Whilst crop system did not significantly influence bacteria and fungi profiles with depth, I S values between CC and CM systems decreased by 0.08 (bacteria) and 0. 16

Discussion
To our knowledge, this is one of the first studies investigating subsoil and bulk soil microbial diversity using high-throughput sequencing in an agriculturally managed Vertosol, in the southern hemisphere [31,68,69]. Whilst our data have followed trends observed in studies from the northern hemisphere (i.e., decline in diversity with increasing depth), there have been a few surprises [30,70]. The first is that depth appears to have a lower level of influence on fungal (ITS) diversity for both alpha-and beta-diversity [33,34,71]. The second is that the imposition of maize into cotton rotations had a limited effect on bacterial α-diversity (significantly noticeable at 0.70-1.0 m) and no effect on fungal α-diversity. However, significant differences between systems was evident in β-diversity measures for bacteria (0-0.15 and 0.30-0.50 m) and fungi (0.30-0.50 m). Finally, sample time did not significantly influence microbial alpha-or beta-diversity, with the exception of bacterial α-diversity with depth, with increased Chao1 values during season 2 at 0.30-0.50 m.
Depth was a significant factor influencing the diversity of bacteria, both within a sample (αdiversity) and between samples (β-diversity), whilst depth was only a significant factor for fungal βdiversity. Despite a significantly higher species richness (SChao1) at 0.30-0.50 m, bacterial α-diversity (DGS) was overall lower as depth increased, with the low evenness values (EGS) indicating low abundance of the different bacterial species present with increasing depth. A decrease in bacterial diversity with increasing depth is not a surprising result, as availability and accessibility to C, N and P sources becomes increasingly difficult with depth (Table S1) in both natural and managed ecosystems [30,[70][71][72][73]. The trend of decreasing bacterial diversity with increasing depth was also apparent in β-diversity analyses, with nMDS showing pronounced separation between each depth assessed. We performed distance-based linear modelling (DISTLM) sequential tests and plotted

Discussion
To our knowledge, this is one of the first studies investigating subsoil and bulk soil microbial diversity using high-throughput sequencing in an agriculturally managed Vertosol, in the southern hemisphere [31,68,69]. Whilst our data have followed trends observed in studies from the northern hemisphere (i.e., decline in diversity with increasing depth), there have been a few surprises [30,70]. The first is that depth appears to have a lower level of influence on fungal (ITS) diversity for both alphaand beta-diversity [33,34,71]. The second is that the imposition of maize into cotton rotations had a limited effect on bacterial α-diversity (significantly noticeable at 0.70-1.0 m) and no effect on fungal α-diversity. However, significant differences between systems was evident in β-diversity measures for bacteria (0-0.15 and 0.30-0.50 m) and fungi (0.30-0.50 m). Finally, sample time did not significantly influence microbial alpha-or beta-diversity, with the exception of bacterial α-diversity with depth, with increased Chao1 values during season 2 at 0.30-0.50 m.
Depth was a significant factor influencing the diversity of bacteria, both within a sample (α-diversity) and between samples (β-diversity), whilst depth was only a significant factor for fungal β-diversity. Despite a significantly higher species richness (S Chao1 ) at 0.30-0.50 m, bacterial α-diversity (D GS ) was overall lower as depth increased, with the low evenness values (E GS ) indicating low abundance of the different bacterial species present with increasing depth. A decrease in bacterial diversity with increasing depth is not a surprising result, as availability and accessibility to C, N and P sources becomes increasingly difficult with depth (Table S1) in both natural and managed ecosystems [30,[70][71][72][73]. The trend of decreasing bacterial diversity with increasing depth was also apparent in β-diversity analyses, with nMDS showing pronounced separation between each depth assessed. We performed distance-based linear modelling (DISTLM) sequential tests and plotted distance-based redundancy analysis (dbRDA) on a subset of diversity data against environmental factors, and both DOC and gravimetric water content (GWC) were equally proportional (29.9% each) in explaining 73.6% of total variation ( Figure S1) in bacterial diversity. It appears that bacterial communities in the top 0-0.15 m of the soil profile are more positively correlated to DOC, whilst subsoil bacterial diversity (i.e., 0.30-0.50 and 0.70-1.0 m) was negatively correlated to DOC, despite subsoils having a higher GWC. Sørensen's Index of Similarity (I S ) indicated that both Proteobacteria and Actinobacteria dominated the top 0-0.15 m of the soil profile, an observation consistent with these copiotrophs preferences for nutrient-rich environments [74][75][76]. Interestingly, the relative abundance of Actinobacteria increased with depth, whilst the relative abundance of the oligotrophic (adaptable to low-nutrient habitats) Acidobacteria [75,76] decreased with depth at 0.7-1.0 m, in an apparent negative correlation. There have been reports that water-extractable organic carbon (i.e., DOC) has decreased availability in soils with increased clay content, and therefore the DOC content in moist soil conditions does not indicate the availability of this C to microbiota [77,78]. It is, therefore, plausible that the decrease in DOC accessibility down the soil profile resulted in overall lower bacterial diversity in subsoils, which was monopolized by the highly competitive Actinobacteria [72,79].
The lack of significant differences in fungal α-diversity with depth (with the exception of S Chao1 ) is surprising, as fungal diversity has been reported as being limited to surface soils, where oxygen and OM are abundant [70,71,80]. Granted, S Chao1 indicates a significant decrease in species richness below 0.15 m but, overall, the abundance of these fungal species is low. It is plausible that the high clay content and subsequent cracking in the Vertosol profile facilitates the movement of O 2 , OM and H 2 O down the soil profile, providing micropockets of nutrients to be accessed by fungal hyphae networks [19,81,82]. β-diversity analyses using nMDS highlight the effect of depth on fungal diversity, with communities in the top 0-0.15 m separate from the remaining depths, which were clustered together. In contrast to our observation with bacterial diversity, DISTLM sequential analyses indicated that GWC (48.9%) was the highest contributor to the 75.8% total variation in diversity data on the dbRDA plot ( Figure S2), when environmental factors were taken into consideration. In this context, GWC was negatively correlated with DOC, which in turn was driving the separation of fungal diversity in the top 0-0.15 m of the soil profile. Sørensen's Index of Similarity (I S ) indicated the dominance of Ascomycota in these top soils, largely comprised of members belonging to Sordariomycetes, which are known to prefer warm, moist and well-oxygenated environments [83]. As depth increased, the prevalence of Basidiomycota, (predominantly members of Agaricomycetes) increased, possibly reflecting differences in responses to water and carbon between these two classes of fungi [83].
Crop system influence on α-diversity down the soil profile was only apparent in subsoil (0.70-1.0 m), with a decrease in bacterial diversity (D GS ) under CC prior to planting (pre-crop) in both seasons. However, time was not a significant factor. The decrease in bacterial diversity in the CC system at 0.70-1.0 m depth is mostly likely attributed to the prevalence of maize crop roots in the CM system after crop removal, as observed by Rasse and Smucker [84]. Although not measured as part of these experiments, the presence and distribution of roots have been examined in this rotation and field previously [15,85]. The prevalence of maize roots would create diverse biopore networks, with OM accumulating and being flushed through the soil profile upon irrigation [19,43,86], thus becoming available to the subsoil microorganisms in this system. Another possible contributor to consider is the increased P in the 0.70-1.0 m soils under CM (Table S1). While these soil measurements were taken during season 1, the major P pool in these soils is calcium phosphate, which has a long turnover time and is thus unlikely to have changed [26]. Therefore, increased accessibility of C and increased P could both be co-contributors to increased diversity in CM at 0.70-1.0 m Crop system influence down the soil profile was apparent in β-diversity measures for both bacteria and fungi. The crop system differences were more pronounced for bacteria at 0-0.15 m, but a clear dispersion between system mean values at 0.30-0.50 m was also observed. In contrast, system differences occurred only at 0.30-0.50 m for fungi. The drivers behind microbial diversity patterns between the two investigated systems remain unclear. The DISTLM analysis and dbRDA plots demonstrated that DOC and GWC were key players in microbial diversity, but there was little separation between systems. However, it seems that fungi diversity at 0.30-0.50 m in CC was more positively correlated with DOC, compared to fungal diversity in CM at the same depth ( Figure S2). A further observation that is worth mentioning is the distinct increase in Basidiomycota at 0.70-1.0 m under CM, as indicated by Sørensen's Index of Similarity. Whilst not a significant increase, it is highly possible that the Basidiomycota have proliferated in this region due to the persistence of lignin containing maize roots in subsoils, acting as a C source for members of this phyla with lignin degrading metabolic pathways [15,84,85,87,88].
Sample time (i.e., prior to planting, mid-crop development, post crop removal, and rotational phase) has been linked to belowground changes in soil microbiota diversity [14,[89][90][91][92]. However, this was not evident down the soil profile in the systems we studied, with the exception of bacterial α-diversity at 0.30-0.50 m during season 2. Season 2 was characterised by temperatures falling outside the ideal temperature range (11-36 • C) for cotton growth and drought setting into the region [40]. This resulted in poor crop establishment, which, in turn, would have had a knock on effect for root growth and soil microbiota, thus disrupting the structure of microbial communities [14,89]. The increased S Chao1 values at 0.30-0.50 m could be reflective of previously suppressed bacterial species proliferating in favourable conditions, which would contribute to the increased D GS values, although these were not significant increases [70,93]. Whilst we cannot definitively state that microbial communities are resistant to the influence of agricultural disturbance, without post crop diversity comparisons, there does appear to be a level of resilience in the microbial communities studied [10,94,95].
This study represents an initial exploration into Australian subsoil microbial diversity in cotton systems, establishing a baseline for understanding these vital biota and their resilience in highly dynamic soils. This work provides a framework against which to explore the functional diversity of subsoil microbiota in future research, which may assist in defining how the addition of maize has influenced changes in microbial diversity, which remains unclear [39,96,97]. Our research investigating both long-and short-term microbial activity in these soils suggest that the differences between systems is likely a result of how microbiota respond and allocate C and other nutrients for metabolic processes [26,71,98].

Conclusions
The majority of agricultural soil microbial diversity studies have focused on the top 0-0.30 m, with exploration into subsoil (>0.3 m) limited and microbial diversity in rotational cotton systems cultivated in Vertosols non-existent. Vertosols are highly dynamic, with both vertical and lateral soil movements facilitating the downward movement of organic matter, oxygen and water, with the potential to support microbial life beyond the top 0.3 m. Despite fungal diversity patterns deviating from published research, our first hypothesis, relating to decreasing microbial diversity with increasing depth, was upheld. Our second hypothesis that microbial diversity would increase down the soil profile under the cotton-maize system was partially refuted in relation to α-diversity, with significant differences only occurring at 0.70-1.0 m for bacteria. Significant differences in β-diversity between systems was evident. However, the exact mechanisms behind diversity differences were masked by the strong depth effect. Our final hypothesis, that microbial diversity would increase in similarity when both systems were under the same crop rotation (cotton), was neither refuted nor supported, as there were no significant differences between the systems at any given time point, indicating a level of resilience in microbial communities to agricultural management.
The dynamic nature of these soils is likely to be responsible for the general lack of difference observed in the α-diversity of microbial communities between the continuous cotton and cotton-maize rotations. Differences with depth, within the bacterial and fungal β-diversity, appear to be attributable to changes in DOC and available water, both resulting from changes in plant rooting patterns within this dynamic soil. The distinct lack of diversity between sample times as influenced by season and crop status is reflective of the resilience in the microbial communities studied. The data presented in this paper are, to our knowledge, the first investigation into Vertosol subsoils in agricultural crops, using high-throughput sequencing. As such, we present a baseline for future investigations into subsoil diversity of Australian agricultural soils and indeed Vertosols.