Yes sorry, merge actually makes no sense, I'm not sure why I thought it would work.
The below code should keep writing the outputs into your workspace, with the last feature class being the one you need e.g. "SP_Test20". You may want to adapt it to process these in memory and then return only the final feature.
Also as Blake say's - remove fc1 from your workspace and put it into a separate geodatabase, as it will be iterating over itself (or put some logic i.e. if fc !=... : , but I would say the best thing to do is just take it out of there).
This is very simple and won't account for removing shape_length, shape_area etc. or any field mappings and join rules, it simply returns the first matching feature.
import arcpy, os
from arcpy import env
from datetime import datetime as d
startTime = d.now()
arcpy.env.workspace = r"C:\ArcGISHome\multipleSpatialJoin\b.gdb"
fc1 = r"C:\ArcGISHome\multipleSpatialJoin\a.gdb\a"
arcpy.env.overwriteOutput = True
fcs = arcpy.ListFeatureClasses() #gets the feature classes
counter = 1
for fc in fcs:
print(fc)
out_fc = "SP_Test" + str(counter)
arcpy.SpatialJoin_analysis(fc1,fc, out_fc, "JOIN_ONE_TO_ONE", "KEEP_ALL", "", "INTERSECT", "", "")
counter += 1
fc1 = out_fc
print ('(Elapsed time: ' + str(d.now() - startTime)[:-3] + ')')