I would like to do a reverse clip (or erase in ArcMap) from shapefiles (polygons) then to merge the two shapefiles into a single shapefile. From this post Reverse clipping (erasing ) in R?, I tested the function gDifference from the rgeos package but it's very long (several days). How can I improve the speed or is there another way to do a reverse clip in R?

Here are some information on the two shapefiles (I kept the ArcMap terminology http://pro.arcgis.com/en/pro-app/tool-reference/analysis/erase.htm)

## Input feature
> shp1
class       : SpatialPolygonsDataFrame 
features    : 61378 
extent      : -74, -71.5, 45, 45.875  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
variables   : 55

## Erase feature
> shp2
class       : SpatialPolygonsDataFrame 
features    : 429418 
extent      : -79.53649, -61.75487, 44.99248, 50.25604  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
variables   : 4

I would like to do a reverse clip (or erase in ArcMap) from shapefiles (polygons) then to merge the two shapefiles into a single shapefile. From this post Reverse clipping (erasing ) in R?, I tested the function gDifference from the rgeos package but it's very long (several days). How can I improve the speed or is there another way to do a reverse clip in R?

Here are some information on the two shapefiles (I kept the ArcMap terminology http://pro.arcgis.com/en/pro-app/tool-reference/analysis/erase.htm)

## Input feature
> shp1
class       : SpatialPolygonsDataFrame 
features    : 61378 
extent      : -74, -71.5, 45, 45.875  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
variables   : 55

## Erase feature
> shp2
class       : SpatialPolygonsDataFrame 
features    : 429418 
extent      : -79.53649, -61.75487, 44.99248, 50.25604  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
variables   : 4