This part of the code

    x2subject <- as(hits, "IntegerList")
    x2d <- relist(mcols(hits)$distance, x2subject)
    hits[unlist(min(x2d) == x2d)]

was to me highly idiosyncratic (in retrospect I guess there was only one thing that I didn't recognize immediately, which was the coercion of hits to IntegerList, but that was enough to send me off on the follow...). I started with a 'base R' solution (using 'query' instead of 'x' for consistency with the argument names of distanceToNearest()). Here's the distance between each 'nearest' query  / subject pair

    distance <- distance(query[queryHits(hits)], subject[subjectHits(hits)])

The challenge is to apply a function (which values are the minimum?) to the distances, grouped by query hit. If I didn't have any groups, I could use fun1() to return a logical vector indicating the values that are minima 

    fun1 <- function(x)
        min(x) == x

To apply this function to grouped values, a base R solution is to use the function ave()

    keep <- as.logical(ave(distance, queryHits(hits), FUN = fun1)

and then to clean up results

    ## clean-up
    mcols(hits)$distance <- distance
    hits[keep]

distanceToNearest() doesn't actually return ties, so I modified the function used in ave() to

    fun2 <- function(x)
        min(x) == x & !duplicated(x)

I think the performance of ave() wouldn't be too bad for a typical script, and the overall code seems to me easier to understand (ave() is probably relatively specialized knowledge, but the problem of operating on grouped values in base R comes up quite frequently on StackOverflow, so once one recognizes the problem being faced -- apply a vectorized function to a vector split into groups -- a solution isn't too far away).

One bag of tricks I do know about S4Vectors is that a list-of-integers can be efficiently created, manipulated 'group-wise', and turned into a vector. So I replace the use of ave() with

    lst <- splitAsList(distance, queryHits(hits)
    keep <- unsplit(fun1(lst), queryHits(hits))

S4Vectors applies fun1() to each element of lst, i.e., group-wise. Maybe this isn't that much less idiosyncratic than Herve's -- my first line is Herve's lines 1 & 2, and the second line is the subset in Herve's line 3. But knowing the splitAsList() / unlist(fun()) trick applies to many different situations, whereas as(hits, "IntegerList") is very specific. So my revision is

    distanceToNearest3 <- function(query, subject, fun = fun2) {
        hits <- union(
            precede(query, subject, select="all"),
            follow(query, subject, select="all")
        )
        ## minimum distance by query
        distance <- distance(query[queryHits(hits)], subject[subjectHits(hits)])
        lst <- splitAsList(distance, queryHits(hits))
        keep <- unsplit(fun(lst), queryHits(hits))
        ## clean-up
        mcols(hits)$distance <- distance
        hits[keep]
    }

Yes, that's right, if a query overlaps a subject range, no hits for it are returned. Depending on exactly what is needed here, there may not be a good solution. The nearest algorithm, as I recall, is greedy and based on positions, so it would need to be rewritten to support this mode. Perhaps the asker could describe their concrete use case. The bedtools closest -io command is documented to behave like Val's solution, where subject ranges with any overlap in the query are ignored.