Rasterio convert vertical datum of NAVD88 to WSG84

Question

I've been asked to provide an elevation data for projected to GPS. To achieve that i've tried to convert the NAD83 geotiff with NAVD88 vertical datum (downloaded from USGS) into WSG84.

To do the conversion i'm using rasterio:

import rasterio
from rasterio.warp import Resampling, calculate_default_transform, reproject
dst_crs = "EPSG:4326"
base_filename = "test_raster"
path_to_original_geotiff= f"./{base_filename}.tif"
path_to_reprojected_geotiff= f"./reprojected_{base_filename}_4326.tif"
with rasterio.Env(
    CPL_DEBUG=True,
):
    with rasterio.open(path_to_original_geotiff) as src:
        print(src.crs)
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            "crs": dst_crs,
            "transform": transform,
            "width": width,
            "height": height
        })
    with rasterio.open(path_to_reprojected_geotiff, "w", **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

I've managed to do the horizontal transition correctly, however i do not see any data changes after the conversion. The min and max value of the file stays intact.

I have PROJ [7.2.1] and GDAL [3.2.2] installed in my env with PROJ_LIB variable set. In the directory where PROJ_LIB is pointing to I've put a g2012b geoid files downloaded from: https://download.osgeo.org/proj/vdatum/

What could i do to properly convert those files into WSG84 with rasterio?

--EDIT

My input CRS EPSG:4269. Acc. to the source link the USGS data relates to NAD83 (horizontal) and NAVD88 in vertical.

Here is a manual example of a transformation i expect to happen done via VDatum

Output of gdalinfo on my geotiff:

Driver: GTiff/GeoTIFF
Files: test_raster.tif
Size is 10812, 10812
Coordinate System is:
GEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101004,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4269]]
Data axis to CRS axis mapping: 2,1
Origin = (-99.000555556093445,28.000555556094923)
Pixel Size = (0.000092592592692,-0.000092592592692)
Metadata:
  AREA_OR_POINT=Area
  BandDefinitionKeyword=*
  DataType=*
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
  LAYOUT=COG
Corner Coordinates:
Upper Left  ( -99.0005556,  28.0005556) ( 99d 0' 2.00"W, 28d 0' 2.00"N)
Lower Left  ( -99.0005556,  26.9994444) ( 99d 0' 2.00"W, 26d59'58.00"N)
Upper Right ( -97.9994444,  28.0005556) ( 97d59'58.00"W, 28d 0' 2.00"N)
Lower Right ( -97.9994444,  26.9994444) ( 97d59'58.00"W, 26d59'58.00"N)
Center      ( -98.5000000,  27.5000000) ( 98d30' 0.00"W, 27d30' 0.00"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  NoData Value=-999999
  Overviews: 5406x5406, 2703x2703, 1352x1352, 676x676, 338x338
  Metadata:
    DESCRIPTION=*
    LAYER_TYPE=athematic
    RepresentationType=*