Thank you for the interesting question!
First, let's consider the easy part. Given that you want the distances from every vertex for all distances, you are asking for all-pairs shortest paths.
Therefore you need to find $|V|(|V|-1)/2$ distances between pairs of vertices, and you cannot do better than essentially $c|V|^2/2$ steps for some small constant $c$.
Your restriction to low bounded degree $d$ implies that $|E| \le d|V|/2$. For $d\le 8$ this means $|E| \le 4|V|$.
As Shir points out in a comment, running BFS starting from each vertex therefore takes at most $C(|V|+|E|)|V| \le C(|V| + 4|V|)|V| = 5C|V|^2$ steps for some small constant $C$. The actual running time will probably be less, since your $d\le 8$ is an upper bound; depending on the graph the multiplier may be closer to the average degree.
So the solution will have an unavoidable $|V|^2$ factor, multiplied by a constant between $c/2$ and $(1+(d/2))C$.
That is as far as a purely asymptotic analysis goes.
However, you indicate you need to do this millions of times (presumably for different graphs), so even a small improvement is likely to be worthwhile, and that it would be worth improving the constants.
Assuming that $c$ and $C$ are similar, the question is then whether it is possible to reduce the multiplier from the $(d/2)+1$ upper bound (which seems to be around $3$ in your setting), to get it closer to the optimal $1/2$.
Depending on the structure of the graph, a modified BFS may give a lot of information about intermediate shortest paths. If you start a BFS at $u$, then see vertex $v$ in the search at distance $d$, then see vertex $w$ at distance $d+e$ in the subtree rooted at $v$, you can conclude that $d(u,w)=d+e$ and $d(u,v)=d$, but also that $d(v,w)=e$. In short, one can extract some of the entries $d(v,w)$ in the distance matrix (where $v$ and $w$ are different from $u$) from the shortest path tree rooted at $u$. This seems to be essentially the idea of this SO answer to a related question.
Whether the cost of extracting additional information out of the shortest path trees is worth the cost depends on the structure of your graphs. For this to be worthwhile, it should allow a lot of information to be extracted from the first few BFS searches, so that many of the later vertices from which a BFS would be started already have their entries in the distance matrix filled in and can be skipped.
The following question then seems interesting:
How should a "good" set of starting vertices be chosen?
With "good" I mean a small set $U$ of vertices so that after modified BFS has been run with each $u\in U$ as the starting vertex, then all (or nearly all) of the remaining vertices all have their shortest path vectors fully computed. Examining some low-degree graphs this seems to be an interesting idea: for a path, either endpoint allows a single BFS to compute all the distance vectors; for a cycle, taking two vertices far apart works. For trees, a single BFS starting from an internal vertex $u$ can be used to derive the distance matrix, since if $d(v,w)$ is not known at the end of the BFS, then $v$ and $w$ must be descended from different children of $u$, and then $d(v,w) = d(v,u)+d(u,w)$ which are both known.
The Nash-Williams decomposition of a graph into edge-disjoint forests might then be relevant. The arboricity of a graph with maximum degree $8$ seems to be at most $4$ (I might be wrong here; however, David Eppstein confirms that arboricity of bounded degree graphs is bounded). So it might be possible to find a small number of vertices from which BFS search trees yield all relevant shortest path information, i.e. so that the distance matrix can always be computed using these trees, and so that the edges of the BFS search trees cover the graph. I do not know how large "small" is in this statement, or even if it can always be made to be $4$.
I didn't find anything solving this particular version of the problem when searching, but expect it has been thought about already.
There is a recent paper which uses a distributed Nash-Williams decomposition to obtain maximal independent sets.
Perhaps someone else has a pointer to a more directly relevant paper.
- Leonid Barenboim and Michael Elkin, Sublogarithmic Distributed MIS Algorithm for Sparse Graphs using Nash-Williams Decomposition, Distributed Computing 22, 363–379, 2010. doi:10.1007/s00446-009-0088-2 (preprint)