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).
 
You have an extremely large number of polygons. Are these transformed raster data?