Warming inhibits increases in vegetation net primary productivity despite greening in India

India is the second-highest contributor to the post-2000 global greening. However, with satellite data, here we show that this 18.51% increase in Leaf Area Index (LAI) during 2001–2019 fails to translate into increased carbon uptake due to warming constraints. Our analysis further shows 6.19% decrease in Net Primary Productivity (NPP) during 2001–2019 over the temporally consistent forests in India despite 6.75% increase in LAI. We identify hotspots of statistically significant decreasing trends in NPP over the key forested regions of Northeast India, Peninsular India, and the Western Ghats. Together, these areas contribute to more than 31% of the NPP of India (1274.8 TgC.year−1). These three regions are also the warming hotspots in India. Granger Causality analysis confirms that temperature causes the changes in net-photosynthesis of vegetation. Decreasing photosynthesis and stable respiration, above a threshold temperature, over these regions, as seen in observations, are the key reasons behind the declining NPP. Our analysis shows that warming has already started affecting carbon uptake in Indian forests and calls for improved climate resilient forest management practices in a warming world.


Introduction
It is unequivocal that the anthropogenic Green House Gas (GHG) emissions are the key drivers of global warming since the pre-industrial period. Measurements show that during the last decade, the global average of annual anthropogenic CO2 emissions has reached "the highest levels in human history" at 10.9 ± 0.9 PgC year -1 , 31% of which was stored by the terrestrial vegetation 1 . Terrestrial vegetation acts as a sink in the global carbon cycle through photosynthesis by uptaking the atmospheric CO2 and recycling to its higher energy form, carbohydrate, and storing it as biomass 2 . Recent satellite observations have concluded that the Earth's green cover has increased significantly in the last two decades [3][4][5][6] . The 'greening' is generally measured with the Leaf Area Index (LAI), computed with one-sided leaf area for broadleaf canopies and half of the needle surface area for coniferous canopies 7 . The CO2 fertilization effect 3,4 is the key contributor to the global greening trend. The other contributors include climatic factors 3,4,8 (example: the rising temperature in the high latitudes) and the land use, land cover change 9,10 .
This global greening intuitively enhances vegetation productivity and hence, increases the carbon sink potential. However, the recently published observed 11 and model-driven 12 studies have not found a proportional increase in vegetation productivity at a global scale.
Gross Primary Productivity (GPP) is one of the key indicator for measuring vegetation productivity. GPP is defined as the amount of CO2 captured by the plants in unit time during the time of photosynthesis. GPP has increased only 0.08 %, whereas the global green cover has increased about 0.23 % during the period 2000-2015. The lower increase in GPP is attributed to the atmospheric moisture stress globally 12 . Net Primary Productivity (NPP) is another important indicator for determining vegetation productivity.NPP is the difference between GPP and autotrophic respiration. There exists inconsistency in the global trends between NPP and GPP arising due to increasing autotrophic respiration 11,13 . Climate extremes like droughts and heatwaves are reported to be one of the key reasons behind the decline in terrestrial vegetation productivity [14][15][16] .Temperature and precipitation are reported to be the two most important climatic factors that have control over vegetation productivity 8,[16][17][18] . Zhang et al. 8 hypothesized that increasing temperature drives global NPP decrease, despite a stable GPP for the terrestrial ecosystem. Temperature 19 , total water storage 20 and soil moisture 21  been attributed to CO2 fertilization 22 . Other studies 23,24 with observational data and model output showed that precipitation is the major driver behind NPP variations. Nayak et al. 24 also estimated an increase of NPP at 0.005 PgC year -

Trends of LAI, GPP and NPP
Figure 1 (a1-c1) presents the annual time series of LAI, GPP, and NPP over India. We find that the LAI over India has increased steadily (statistically significant) in the 21 st century, which is consistent with the findings from Chen et al. 5 , where they attributed the greening in India to the increased crop area. The annual time series of GPP and NPP have not shown any statistically significant trend despite this greening; the NPP showed a slight (not statistically significant) decreasing trend. MODIS NPP, GPP, and LAI products are obtained by MODIS radiometric information. There is a high chance that these date sets are affected by covariability. Extensive ground truth verifications are also impossible for India due to the limited number of flux tower locations having long-term datasets. Furthermore, the MODIS products do not consider increasing productivity due to carbon fertilization 27 .
Hence, we used two more datasets: FLUXCOM GPP and NEE products and Contiguous Solar-Induced chlorophyll Fluorescence (CSIF) data products for more robust analysis.
Supplementary Figure 1 shows the climatologies of MODIS GPP, FLUXCOM GPP and CSIF over the Indian landmass. All three products show high photosynthesis during the last half of the monsoon, August and September. However, the values from MODIS maintain high GPP during the post-monsoon till February, while the other two climatologies drop significantly. It is noteworthy that none of these data products are observations and have their own limitations. It is expected that the GPP in India during post-monsoon should be maintained due to high radiations and high soil moisture from monsoon recharge. The same may continue in the winter, a cropping season in India known as Rabi season. MODIS GPP and CSIF dip during pre-monsoon summer because it is a dry, non-cropping season. The FLUXCOM GPP climatology shows an increase from February-March, which is unexpected. We plotted trends for CSIF, FLUXCOM GPP, and FLUXCOM NEE  We visually observed that in regions 1,4 and 5, the LAI trends diverge from NPP trend; to demonstrate this statistically, we generated a scatter plot of standardized LAI and NPP trends (Supplementary Figure 7). We also produced probability density functions (pdf) of standardized LAI and NPP trends using a bivariate kernel density estimate for each region (Supplementary Figure 8). Each point in the scatter plot and in the sample for generating pdfs represents grid. We kept the LAI (Standardized) trend on the x axis and the trend of NPP (Standardized) on the y axis. The density of LAI trend pixels is higher in the first and fourth quadrants in all five regions, signifying widespread greening over these five regions. We also found a trend divergence between LAI and We observed a steady increase in forest cover area in these regions. Cropland has also increased these regions (Supplementary Figure 12 b and f), except for region 4 (Supplementary Figure 12d). Therefore, land use and land cover changes cannot account for the divergence in trend between LAI and NPP over these regions.
We linearly regressed NPP as a function of LAI (NPP = f (LAI)) and also looked into the correlation between NPP and LAI to perform a more detailed analysis of the annual  There is a spatial resemblance between India's warming hotspots and the regions with declining NPP, indicating a strong temperature control on the vegetation productivity.
For a quantitiative understanding of the role of temperature on diverging trend of NPP and LAI, we linearly regressed NPP as a function of LAI (NPP = f (LAI)) and calculated the trend in the residuals for each grid (Supplementary Figure 16). The residual trend is decreasing in nature in regions (1, 4, and 5). In these regions (1,4 and 5), we also observed a significant warming trend (Figure 3a). We also generated the scatter plots between temperature and the residuals of linear regression (NPP against LAI, Supplementary Figure 17 (a,b,c,d,e)).We found that the unexplained (by LAI) components of NPP are negatively correlated with the temperature for all regions.
These plots demonstrate the negative impacts of warming on the NPP.  14 The GPP generally increases with temperature till an optimum value. When the temperature crosses the optimum value, GPP declines 32 while respiration rises 33 , leading to decreasing NPP 34,35 . Tropical forests typically have the aforementioned characteristics, with a decrease of 9.1 megagrams of carbon uptake per hectare per degree Celsius warming in the mean daily maximum temperature in the warmest month 36 . We plot the variations of GPP in the 5 selected regions with temperature in Figure 4. The curves in the subplots are the smoothed variations as obtained using nonlinear kernel regression 37 . From these curves, we obtained the temperature value where GPP reaches the maximum, and from there, it starts dropping. That temperature is termed as optimum temperature. As MODIS provides only Leaf and Fine root Maintenance respiration (LFrMR) , we also used them in Figure   4. For regions 1, 4, and 5, the LFrMR are either stable or increase beyond the optimum temperature. Hence, with warming above the optimum temperature, the net Photosynthesis   Figure. 4 c, e, g,i and k . For regions 1, 4, and 5, there are statistically significant increasing trends in the number of days above optimum temperature (T GPP opt). These trends, along with low photosynthesis above optimum temperature, results in low net carbon uptake by the vegetation. Interestingly, for regions 2 and 3, the LFrMR drops after the optimum temperature; however, the drops are not as high as the drops in GPP. Region 2 is India's warmest and arid region, and probably the vegetation is already adapted to a warm environment. The decrease in LFrMR could be because of this long-term adaptation of vegetation. For regions 2 and 3, the numbers of days exceeding the optimum temperature do not have a trend. Hence, the decline in NPP is not observed in these regions.  Granger Causality, we may identify the causal factors; however, quantifying contributions from a causal variable to an impact variable is difficult. Furthermore, all the meteorological variables are connected, with the possibility of multiple confounding factor(s). Hence, we did not attempt to quantify the exact causal contribution from temperature to productivity. This is one of the limitation of the present causal analysis and probably a generalized limitation of any widely used causal discovery methods.
Carbon fertilization plays a role in improving water use efficiency and primary productivity and shoule be considered as a potential causal factor in Granger Causality analysis.
However, we do not have any available long-term dataset in India that considers carbon fertilization. This is a major limitation of the study. Further, the atmospheric CO2 concentration does not change with a high-frequency variability (monthly scale) that was used in the Granger Causality analysis. Considering annual CO2 will result in a small annual sample insufficient to perform any causal analysis.

Discussion
Recent decades have experienced increasing atmospheric CO2 concentration and widespread warming around the globe. India is not an exception. Increased warming resulted in an increase in water vapor carrying capacity and subsequent increase in the vapor pressure deficit. Stomatal conductance reduces under increased VPD and transpiration increases till a VPD threshold 38,39 . Both of these processes lead to reduced photosynthesis and these impacts have been observed globally 38,39 Global studies showed that the VPD increased by 0.0017 ± 0.0001 kPa year−1 (Data: CRU) the GPP decreased by 0.23 ± 0.09 Pg C year−1 and 0.31 ± 0.11 Pg C year−1 since the 1999, as per the EC-LUE and MODIS models, respectively 40 . On the contrary, model-driven studies show improved WUE [41][42][43] and vegetation productivity under increased CO2 concentration due to carbon fertilization 41,42,44 . The science question remains if Carbon fertilization can make up for the reduced productivity of the vegetation due to warming and the increased VPD.
A synthesis study performed in the IPCC Assessment Report 6, Working Group 1, Chapter 5 1 and Chapter 11 45 showed with a high confidence that the vegetation carbon sink will be less efficient at a very high warming level. Major limitations in such studies are: 1) the models are still inadequate to consider the physiological response of plants to the increased VPD, and 2) the satellite products like MODIS do not consider Carbon fertilization. The opposing impacts of both the factors and the resulting uncertainty in projections are more pronounced at a lower warming level. The literature agrees that high VPD impacts will dominate globally at a higher warming level over Carbon fertilization due to multiple direct and indirect effects 1 .
The lack of ground-truth data availability in India is one of the major bottlenecks in performing a detailed analysis of vegetation productivity trends and responses to increased VPD and carbon fertilization. The primary contributor to greening in India is agricultural expansion 5 . Recent studies show decreased crop yield due to warming 46 . This finding is probably one of the proxies that confirm the warming-induced declining vegetation productivity in India and is in agreement with our conclusions from the present analysis.
In the last two decades, the green cover of India has increased [MoEFCC, 2021; https://unfccc.int/sites/default/files/resource/INDIA_%20BUR-3_20.02.2021_High.pdf"] and Indian forests are still a net sink of carbon 47 ; however, our study finds that this sink may be weakening due to warming. The NPP from temporally consistent Indian forests shows a significant declining trend of -0.75±0.58 Tg C year -2 and has declined by 6.19% were the good samples per grid, averaged across the grids, considered during the study period. We also presented a spatial map of overall good pixel count at the grid level in the Supplementary Figure 21. We discovered that, on average, more than 85% of the pixels in each grid were taken from 2001 to 2019. Except for a very few grids, almost all the grids have more than 70% good pixels. Since good pixel count never reaches low (<50%) for any grids, we have not performed any sensitivity analysis. The data is processed as per the user guidelines 7 . These 8-day data sets of LAI are converted into monthly and annual datasets.

MODIS GPP and PSNnet Product:
In this study, we used MOD17A2HGF Version 6 Gross Primary Productivity (GPP) product and Net photosynthesis (PSNnet) product, with the spatial and temporal resolutions are 500 meters (m) and 8 days respectively, for calculating vegetation productivity. MOD17A2HGF Version 6 achieved stage 3 validation 48 , this current version has not validate globally but previous versions of MOD17A2H GPP datasets are validate globally across different biomes with Fluextower GPP in many studies 49,50 (R 2 > 0.50). Net photosynthesis (PSNnet) is the difference between GPP and leaf and root respiration. We used the 8-day cumulative GPP and PSNnet product. These datasets also contain a Data quality control layer which we applied for removing the cloud-covered and bad pixels. We have not performed any sensitivity analysis as mentioned earlier.We processed the data as per the user guidelines 48 and converted these 8-day data sets of GPP and PSNnet into monthly and annual datasets.

MODIS NPP Product:
We obtained the Net Primary productivity from MOD17A3HGF Version 6 product; it gives us annual net primary productivity at 500 meters (m) spatial resolution. Previous version of MOD17A3H dataset are consistent with Fluxtower NPP reflecting R 2 value of 0.56 51 . Annual NPP is calculated by taking the difference between GPP and Autotrophic respiration (RA) which is the sum of Growth respiration (RG) and Maintenance respiration (RM). This data is processed as per the user guidelines 52 . Forests are considered as the forest cover. For Cropland, we considered the grids with 60% area as cultivated 56 .
Precipitation and Temperature Data : Daily gridded precipitation data at 0.25°× 0.25° resolution provided by India Meteorological Department (IMD) is used in this study 57 .
Daily average surface temperature is also obtained from IMD at 1°×1° resolution 58 Monthly average and annual average precipitation and temperature are then obtained from these daily datasets.
Vapour Pressure Deficit (VPD) Calculation : Vapour Pressure Deficit is the difference between saturation water vapour pressure and actual water vapour pressure. We calculated VPD by using Teten formula (Equation 1) 59

Analytical solution of MODIS Leaf and Fine root Maintenance Respiration equation:
We have looked into the equations for calculating MODIS NPP 1 and solved them analytically to investigate the effect of climate variables (mainly temperature) on net primary productivity. We have mainly focused on maintenance respiration since it is highly controlled by temperature. Putting the values of a, b in Eq. 6 and Q10 in Eq. 9, we have plotted the temperature response of 1 × ( ) and 1 × ( ) separately in the Supplementary Fig 18 (a) and (b) respectively.
We discovered that initially, Leaf MR rate increased with increasing temperature (Supplementary Fig 18a), but after a certain temperature limit, it began to fall. The temperature limit where Leaf MR rates fall will vary depending on biome type. We have taken cl as constant since it is not a function of temperature. However, the Froot MR rate keeps on increasing with an increase in temperature (Fig S18b). The rate of increase will vary across biome types but will always be positive since cr > 0. We have taken cr as constant since it is not a function of temperature.
For analyzing the collective effect of both leaf and froot MR rate changes with temperature, we have taken We have plotted the temperature response of ( ) in the Supplementary Fig 18c. We can observe that, initially, leaf and respiration rates increase with an increase in temperature, but after a certain temperature, they remain stable. The temperature at which the rate of change will be stable will also vary across different biome type.