SARIMA modelling results. Choosing the right lag for seasonal data After differentiating a monthly climatic data with a lag of 12, and being sure that, at least, one more differentiation will turn my series into white noise (ndiffs only results in 1)...I used auto.arima to model my series and I got a result like $ARIMA(p,1,q)(0,0,Q)_{12}$. 
I am hesitating why $d=1$ if data has been previously differentiated.
Answer: My data was differentiated without taking logarithms. Variance got better with that! 
How can I make sure I choose a good lag for seasonal data?
Answer: try with a subseries plot ggsubseriesplot().
 A: I took your 456 monthly values into AUTOBOX which examines series by ITERATIVELY considering anomalies , power transformations , transience in EITHER model parameters or model error variance over time .... https://autobox.com/pdfs/ARIMA%20FLOW%20CHART.pdf
It developed a model which broadly speaking which does not use a regular differencing as yours did BUT it did introduce the need for a logarithmic transformation.
Following is the identified model  and here  ... (2,0,0)(1,1,0)12 in logs 
The ambiguity that you encountered was to apply regular differencing or not . The acf of the seasonally differenced data is here suggesting an ar(2) for short term memory NOT regular differencing . The message here is that the conclusion not to regularly difference and ar/ma structure should be made AFTER the seasonal differencing is in place.
The two pulses are of minor importance but might be of interest as they reflect values that were identified as extraordinary. AUTOBOX detected an ar(12) factor while you used an ma(12) .
The recent actuals and forecasts are here with the residual acf here suggesting model sufficiency 
