Here is a short example. I assume that by overlay you are looking for intersects; namely to dissolve all polygons that intersect from both layers.

    library(sp)
    library(rgeos)
    # Create a dataset
    poly <- SpatialPolygons(list(
      Polygons(list(Polygon(coords = matrix(c(1, 1, 4, 3, 4, 2, 1, 1), ncol =  2, byrow = TRUE))), ID = "1"),
      Polygons(list(Polygon(coords = matrix(c(3, 1.5, 3.5, 2, 4, 2, 3, 1.5), ncol = 2, byrow = TRUE))), ID = "2"),
      Polygons(list(Polygon(coords = matrix(c(4, 1, 5, 2, 5, 1, 4, 1), ncol = 2, byrow = TRUE))), ID = "3")
    ))

    # Split dataset to two SpatialPolygon objects
    polyA <- poly[1, ]
    polyB <- poly[2:3, ]

    # Show dataset
    plot(poly, axes = TRUE)
    plot(polyA, axes = TRUE, add = TRUE)
    plot(polyB, axes = TRUE, add =TRUE, col = "red")


You can see the sample data set below. `PolyB` has two polygons from which one intersects the polygon in `PolyA`.

[![Dataset basic plot][1]][1]

Using `gUnion` will result in one polygon of all polygons in both objects, as you have suggested: `plot(gUnion(polyA, polyB))`.

Yet, if you select only those polygons that intersects `polyB[polyA, ]`, dissolve will give you the expected result:

`plot(gUnion(polyA, polyB[polyA, ]))`

[![Union overlay][2]][2]

You can and should subset the first layer as well, `gUnion(polyA[polyB, ], polyB[polyA, ])`, if it has more than one feature.

    

    


  [1]: https://i.sstatic.net/LfFhz.png
  [2]: https://i.sstatic.net/ACyzu.png