Adaptation in U.S. Corn Belt increases resistance to soil carbon loss with climate change

Increasing the amount of soil organic carbon (SOC) has agronomic benefits and the potential to mitigate climate change. Previous regional predictions of SOC trends under climate change often ignore or do not explicitly consider the effect of crop adaptation (i.e., changing planting dates and varieties). We used the DayCent biogeochemical model to examine the effect of adaptation on SOC for corn and soybean production in the U.S. Corn Belt using climate data from three models. Without adaptation, yields of both corn and soybean tended to decrease and the decomposition of SOC tended to increase leading to a loss of SOC with climate change compared to a baseline scenario with no climate change. With adaptation, the model predicted a substantially higher crop yield. The increase in yields and associated carbon input to the SOC pool counteracted the increased decomposition in the adaptation scenarios, leading to similar SOC stocks under different climate change scenarios. Consequently, we found that crop management adaptation to changing climatic conditions strengthen agroecosystem resistance to SOC loss. However, there are differences spatially in SOC trends. The northern part of the region is likely to gain SOC while the southern part of the region is predicted to lose SOC.

. The maps of a) the number of NRI locations in corn/soybean rotation, b) the areaweighted average of clay content, and c) the area-weighted average of sand content in each sub-region. Maps were generated using the R "ggmap" package 1 (Version 2.6.1; http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf).  Figure S2. Changes in annual mean temperature and precipitation of the three GCMs (2041-2070 relative to 1981-2010). Maps were generated using the R "ggmap" package 1 (Version 2.6.1; http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf). Figure S3. Crop grain yield of climate change scenarios (ensemble mean of three GCMs) compared with the baseline scenario (2041-2070). Maps were generated using the R "ggmap" package 1 (Version 2.6.1; http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf). Figure S4. Predicted carbon input and decomposition factor (combined temperature and moisture effect) averaged for 3 GCMs for the baseline scenario (top row), the difference between no adaptation scenario and the baseline (second row), and the difference between the adaptation scenario and the baseline (third row). Maps were generated using the R "ggmap" package 1 (Version 2.6.1; http://journal.r-project.org/archive/2013-1/kahlewickham.pdf). Figure S5. Correlation between the initial prediction of planting dates (DOY) using temperature and NASS reported average DOY for corn and soybean.

Corn
Soybean a) b) Figure S6. The predicted maturity groups for corn and soybean for the 2061-2070 decade. Gray boundaries represent crop maturity distributions reported by the seed companies and the notation numbers were the group numbers in Supplementary Table 1. Maps were generated using the R "ggmap" package 1 (Version 2.6.1; http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf). Table S1. Parameters used for each maturity group. The spatial distribution of maturity groups can be found in Figure S4. Growing degree days (GDD) and biological days (BD) 2 were used for corn and soybean, respectively.

Regional Simulation Setup
The crop field locations were obtained from the National Resources Inventory (NRI) survey, which is a comprehensive long-term survey conducted by the U.S. Department of Agriculture's Natural Resources Conservation Service (NRCS) 11 . Management data such as crop rotation and irrigation are available in the survey. A total of 54,912 point locations were selected and corn/soybean rotation was modeled at each location. Locations were classified as predominantly corn-soybean if at least one-third of the years were cropped in corn and one third in soybeans over the recent history in the NRI survey.
In our simulation, we applied automatic irrigation in the DayCent simulations, if the fields were irrigated according to the survey. Tillage practices were modeled by grouping them into three categories: full tillage (FT; less than 15% residue remaining), reduced tillage (RT; 15-30% residue remaining), and no tillage (NT; above 30% residue remaining) using the data from the 2008 Crop Residue Management Survey conducted by the Conservation Technology Information Center (CTIC) 12 .
Nitrogen fertilizer application rates were calculated in two steps. We first simulated the C dynamics using a routine that automatically adds adequate fertilizer to avoid nitrogen limitations. The fertilizer rates for the second simulation were calculated using the potential crop yield for each decade from the first simulation and recommended application rates from DuPont Pioneer to meet the targeted yield (1.2 lb/bu or 0.021 kg/kg for corn and 1.5 lb/bu or 0.025 kg/kg for yield above 60 bu or 1633 kg for soybean; https://www.pioneer.com/home/site/us/agronomy/library/nitrogen-fertilizer-for-soybean/).
Simulations were executed with customized scripts in R to automatically build DayCent input files, simulate the histories and projections, and analyze the outputs. Simulations were performed using a computer cluster with 440 processors. Simulation results were aggregated to the county level for analysis by weighting NRI points using the expansion factors from the NRI survey, which are based on the two-stage sampling design.

DayCent model
The DayCent model 13 is a widely used biogeochemical model for simulations of soil carbon, soil nitrogen, and crop growth 5,14,15 . The model was used to simulate corn and soybean yields for the contiguous U.S. and found to compare well with NASS reported yield values 16 . Two methods are available for modeling crop phenology. One is the growing degree day (GDD) method which is commonly used by seed companies for corn development (base temperature of 10 C˚ and ceiling temperature of 30 C˚). Another newly added method is the biological days (BD) method 2 which simulates the photoperiod effect on soybean phenology. The model simulates the elevated CO2 effect on crop production, transpiration, root/shoot ratio, and carbon/nitrogen ratio [17][18][19][20] .

Maturity groups, planting dates and maturity dates.
Our study region was the U.S. Corn Belt. We analyzed the current crop management for the entire U.S. in order to determine the adaptation of crops under climate change. To simulate the varieties, we divided the crops into groups (Supplemental Table 1). The comparative relative maturity (CRM) is a common method of labeling the length of time to maturity for corn. Dupont Pioneer Seed Company grouped CRMs spatially (https://www.pioneer.com/home/site/us/agronomy/library/compare-maturitycorn-products/; accessed 05/06/2016). We adapted the map and associated the CRM with GDD 21 . For soybeans, we used the maturity group map from Monsanto (https://www.aganytime.com/asgrow/mgt/planting/Pages/Soybean-Maturity.aspx; accessed 05/06/2016); each maturity group was represented by the difference in their BD to maturity according to the sensitivity of the variety to temperature and photoperiod 2 .
In the adaptation scenario, planting dates and maturity groups were set to change every decade in our simulation. We created a simple regression model to predict the planting dates by analyzing the relationship between air temperature and the NASS reported survey data for corn as corn planting dates are affected by temperature 22 . We first estimated the dates when temperature was not a restriction for planting based on the average weekly air temperature (threshold of 15.4 C) and the average last day of frost (daily minimum temperature below 0 C) 22 for all the major corn-growing counties in U.S. Note that counties south of the Corn Belt were also included in the development of regression as the climate in these areas are similar to the projected climate in the study region. These dates were compared to state average planting dates from NASS 23 to derive the regression equations ( Supplementary Fig. 5a). A similar calculation was conducted for soybean planting dates. Even though the correlation for soybean was lower ( Supplementary Fig. 5b), the RMSE of 10.2 days based on the regression model predictions demonstrated that the effect on grain yields were small for regional simulations. Soybean planting dates are not as strongly correlated with temperature as corn because soybean growth also depends on day length.
We also used different methods for corn and soybean to predict maturity groups. i.e., length of the growing season, for the climate change projections in the adaptation scenario. For corn, which is not affected by day length, we used the first day of frost in the fall/winter (decadal average) to approximate the maturity date. To avoid heat stress at critical development stages, farmers use hybrids that mature long before the first frost dates in regions that are south of the Corn Belt 24 . Therefore, we matched the maximum GDDs from planting to frost date to the reported GDDs across the US (Supplementary Tab. 1) and then applied these regression relationships for future predictions ( Supplementary Fig. 6).
For soybeans, we developed a regression model to predict the physiological maturity (PM) dates. We used NASS reported planting dates and historical distribution of maturity group from data compiled by seed companies to calculate the PM dates for each weather grid cell (Supplementary Tab.1). These PM estimates were treated as observations and were compared to the predicted first day of frost in the fall/winter to develop the linear regression model for predicting PM dates (RMSE 9.3 days). For both the baseline and climate change scenarios, the BD values were calculated using the planting and PM dates. The suitable maturity group for a weather grid was determined for each decade by selecting the maturity group with BD (Supplementary Tab. 1) closest to the predicted BD in that weather grid ( Supplementary Fig. 6).

Model Calibration and Verification
Details of model calibration and verification can be found in our earlier study 25 . Based on our previous sensitivity analysis 26 , crop parameters were either derived from literature or calibrated using experimental sites located in different regions of the U.S. (Supplementary Tab. 2). Growing degree days for corn were calculated from a regression in a previous study 27 . Soybean photoperiod related parameters were taken directly from the work of Soltani and Sinclair 2 . The parameters associated with the effect of elevated CO2 on production and water use were from our research 28 . The SOC decomposition parameters were from the U.S. greenhouse gas (GHG) inventory 29 . DayCent simulated corn and soybean yields in the contiguous U.S. with more accuracy than most other published studies using process-based models, with an overall R 2 of 0.54 and 0.54 for corn and soybean, respectively 25 .
We did not assume any yield potential increase in the scenario simulations. Yields are likely to increase further due to technology development and other factors in the future 33 . However, it was beyond the scope of our study to predict the potential for future crop breeding and genetic improvements to increase the yield potentials for corn and soybeans. Our objective was to determine the role of crop adaptation alone on SOC dynamics with future climate change by parameterizing the model to use current crop varieties. Therefore, this analysis addresses the potential for adaptation through selection of varieties that currently exist in the United States. There may be even greater opportunity to enhance C input and increases in SOC in the future with crop breeding and other genetic improvementsin combination with evolving technologies and management strategies to maintain soil health and meet food security goals.