you can use reduce() to perform union, for example we first make a GRanges object with your coordinates:

sep = GRanges( seqnames = "chr1", ranges = IRanges(start=c(58710367,58710476),end=c(58719391, 58716325)), strand = "+")

resulting:

GRanges object with 2 ranges and 0 metadata columns:
  seqnames               ranges strand
     <Rle>            <IRanges>  <Rle>
[1]     chr1 [58710367, 58719391]      +
[2]     chr1 [58710476, 58716325]      +

and then using reduce() to perform union:

U = reduce( sep )

resulting:

GRanges object with 1 range and 0 metadata columns:
  seqnames               ranges strand
     <Rle>            <IRanges>  <Rle>
[1]     chr1 [58710367, 58719391]      +

Remember that this function is sensitive to strand direction and won't union intervals from "+" and "-" strands, but if strand sets to "*" it works