The description of the -a_nodata
parameter from the documentation of gdal_translate https://gdal.org/programs/gdal_translate.html
-a_nodata [value]
Assign a specified nodata value to output bands. It can be set to none to avoid setting a nodata value to the
output file if one exists for the source file. Note that, if the input
dataset has a nodata value, this does not cause pixel values that are
equal to that nodata value to be changed to the value specified with
this option.
I believe that in your original image value -3.4e+38 is used to mean nodata. Then you fired a gdal_translate command
gdal_translate -a_nodata 0.0 ...
What happens is that all pixel values of 0.0 has been marked as nodata in the target image. If -3.4e+38 used to mean nodata before it does not mean it anymore but all those pixels contain now real data for QGIS and SAGA. The -a_nodata parameter does not change the pixel values, it is just touching the metadata.
Without having test data it is hard to know but it feels like you should have -3.4e+38 as nodata value. If it is not so in the original image (check it with gdalinfo) you can do set if by running
gdal_translate -a_nodata -3.4e+38 ...
EDIT
I downloaded the image from your link and I am correcting myself. The original image has no nodata set.
gdalinfo NitrogenBalanceOnLandscape_140Crops.tif
…
Band 1 Block=432x224 Type=Float32, ColorInterp=Gray
Overviews: 2160x1080, 1080x540, 540x270, 270x135
I set the nodata to 0 with the similar command that you used and now I have
gdalinfo with_nodata.tif
…
Band 1 Block=4320x1 Type=Float32, ColorInterp=Gray
NoData Value=0
I digitized two polygons with QGIS, one over the ocean and one partly over Africa. Then I used the zonal statistics tool from the processing toolbox. The results in the table show zero and null for the all nodata area, and reasonable values for the African triangle. For me it seems that nodata pixels are correctly excluded from the zonal analysis.

EDIT 2
By following the workflow that is described in the question I am getting the same results. For some reason that I can't explain the workflow original image -> QGIS raster calculator (layer1 > 0) * layer1 -> assign nodata=0 with gdal_tranlate results an image that has 4320 pixels with value -3.4028234663853e+38
Min=-340282346638529993179660072199368212480.000 Max=2752376.250
Minimum=-340282346638529993179660072199368212480.000, Maximum=2752376.250, Mean=-2090745287322099897642461654474031104.000, StdDev=26590834825873999547824535346249465856.000
256 buckets from -3.4095e+38 to 6.6722e+35:
4320 0 0 0 0 0 0 0 0 0 0 0 0 0
I can get correct result with gdal_calc https://gdal.org/programs/gdal_calc.html but even here I can't explain why.
The command is
gdal_calc -A nitro.tif --outfile nitro_calc.tif --calc="(A>0)*A" --NoDataValue=0
Compare the resuls from the original workflow and the gdal_calc workflow (columns _ccount, _csum, and _cmean)

Understanding what happens requires probably very deep knowledge on how QGIS, GDAL, and numpy work. Perhaps you can reach best experts from qgis-users or qgis-developers mailing list.
EDIT 3
The origin of the issue is probably in the QGIS raster calculator. The lowest pixel row of the original image contains just 0 values and raster calculator with (layer1 > 0) * layer1
should not do anything for the lowest row. But actually it considers that the lowest row is nodata and changes the values which used to be 0 into -3.4e+38. Later GDAL is told to label 0 as nodata and then -3.4e+38 turns into real data and pixels are used in the zonal statistics.
