gdal_merge.py is not working with NODATA values

Question

I am trying to merge a couple of GeoTiff files with gdal_merge.py. Example input:

In the input data the nodata value is -32678. Hence I tried:

gdal_merge.py -n -32678 -o DEM.tif -of GTiff na_DEM.tif ca_DEM.tif

and

gdal_merge.py -n -32678 -a_nodata -32678 -o DEM.tif -of GTiff na_DEM.tif ca_DEM.tif

However, for some reason the nodata values are not recognized as the output file giving me erroneous values where the two files overlap. Below I provided the output from gdalinfo on one of the input files.

Driver: GTiff/GeoTIFF
Files: ca_DEM.tif
       ca_DEM.tif.aux.xml
Size is 14160, 6000
Coordinate System is:
GEOGCS["WGS 84",
  DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
          AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
  PRIMEM["Greenwich",0],
  UNIT["degree",0.0174532925199433],
  AUTHORITY["EPSG","4326"]]
Origin = (-119.000000000000000,29.999999999999002)
Pixel Size = (0.004166666666667,-0.004166666666667)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
 Upper Left  (-119.0000000,  30.0000000) (119d 0' 0.00"W, 30d 0' 0.00"N)
Lower Left  (-119.0000000,   5.0000000) (119d 0' 0.00"W,  5d 0' 0.00"N)
Upper Right ( -60.0000000,  30.0000000) ( 60d 0' 0.00"W, 30d 0' 0.00"N)
Lower Right ( -60.0000000,   5.0000000) ( 60d 0' 0.00"W,  5d 0' 0.00"N)
Center      ( -89.5000000,  17.5000000) ( 89d30' 0.00"W, 17d30' 0.00"N)
Band 1 Block=14160x1 Type=Int16, ColorInterp=Gray
  Min=-42.000 Max=4065.000 
  Minimum=-42.000, Maximum=4065.000, Mean=836.674, StdDev=815.941
  NoData Value=-32768
  Metadata:
    STATISTICS_MAXIMUM=4065
    STATISTICS_MEAN=836,67402450978
    STATISTICS_MINIMUM=-42
    STATISTICS_STDDEV=815,94103514887