How can I combine USGS GeoTIFF quadrangles into a single GeoTIFF using `rasterio`?

There are several similar questions about this but most as about QGIS and are non-specific about the data involved. I have a specific data source/format and I need to use Python, and preferably rasterio.merge to do the heavy lifting.

The files in question are USGS 2020 quadrangles, downloaded from https://ngmdb.usgs.gov/topoview/viewer/#4/40.01/-100.06 in GeoTiff format

I need to combine 4 of them (with a common corner) into a larger, seamless GeoTiff. I have succeeded at this in the past using rasterio.merge similar to the following

import rasterio
from rasterio.merge import merge
from rasterio.plot import show

quadrangles = [rasterio.open(f, 'r') for f in files]
combined, new_transform = merge(quadrangles)
show(combined, cmap='terrain')

But in this case it results in one image (complete with legend, scale, header, footer) stacked on top of another. When I try this with four images, two of them seem to get completely lost.