This is how you can do that:
library(raster)
# erase
x <- shp1 -shp2
# union (append)
y <- x + shp2
But it won't be faster than gDifference
(which is used under the hood --- the benefit of the method shown here is that attributes are not lost).
Perhaps this speeds things up a bit
agg <- aggregate(shp2)
x <- shp1 -agg
y <- x + shp2
You have an extremely large number of polygons. Are these transformed raster data?