Searching multiple geodatabases for feature classes and merging the results

I’ve got some Sea Level Rise data from NOAA (https://coast.noaa.gov/slrdata/) stored in 66 different geodatabases. Each geodatabase represents a state or district and contains a set of feature classes of sea level rise at different depth (0 feet, 1 feet,…). As is shown in these pictures:

My aim is to merge all the features that have the same water depth, for example, 0 feet, to a single feature class to create a sea level rise map of the US at different water depths. So I used Arcpy to find polygon feature classes that end with "slr_0ft" through these geodatabases. And then merge the results to one single feature class. Here’s the code I used:

import arcpy, os
workspace = "G:DirectStudySLRgdb" #My workspace that contains 66 GDBs
search = "slr_0ft"
feature_classes = []
for dirpath, dirnames, filenames in arcpy.da.Walk(workspace,
                                                  datatype="FeatureClass",
                                                  type="ANY"):
    for fname in filenames:
        # search for string in string to eliminate the need for exact filenames
        if search.upper() in fname.upper():
            feature_classes.append(os.path.join(dirpath, fname))
    for fc in feature_classes:
         print fc

arcpy.Merge_management(feature_classes, "G:DirectStudySLRMerge_0ft")  
print "Done!!"

The code was successful in printing paths and filenames (see below). But the output file from the merge function is an empty shapefile. I checked it with ArcMap and found it contains no polygon. Also, the script didn’t give me any kind of error or warning. I guess the mistake should be in the arcpy.Merge_management() function. How can I fix my script?

I’m not a native speaker so I tried best to explain my question clearly.