Wow, that's really interesting. The same seems to be true for base 7; I haven't tried any other prime bases yet (and I don't know how it might extend to non-prime bases). Anyone have an idea why this seems to hold?

Edit: This basically works for base 10, too. I feel like the reason must either be very obvious or very deep.

I also tested in base 10, with primes under 1e9, got this:

             1       3       7       9
  1       0.0458  0.0746  0.0756  0.0540
  3       0.0599  0.0439  0.0707  0.0755
  7       0.0638  0.0677  0.0438  0.0747
  9       0.0805  0.0638  0.0599  0.0458

It appears to have the same transition probabilities as in base 5 (with the two center rows and columns swapped):

             1       2       3       4
  1       0.0458  0.0756  0.0746  0.0540
  2       0.0638  0.0438  0.0677  0.0747
  3       0.0599  0.0707  0.0439  0.0755
  4       0.0805  0.0599  0.0638  0.0458

I share your feeling about it being either obvious or deep.

Well, whether it's obvious or not, I think it's in the paper. Immediately under equation (1.1) of http://arxiv.org/abs/1603.03720, the authors are discussing the second correction term to the distribution of primes mod q, and they state:

"We can also show that c2(q; (a, b)) = c2(q; (−b,−a)) for any two reduced residue classes a and b (mod q)."

I'm not 100% certain this is responsible for the phenomenon that we're seeing, but it seems exceedingly likely. I think I'd need to stare at their formula for c2 for a long time to understand where this relation comes from, though.

I'm pretty solidly convinced of the pattern at this point: I've checked the "reflect across the anti-diagonal" pattern for the first 10 million primes, expressed in every base from 3 to 20, and it seems to hold up. (I haven't tried to establish any sort of bounds.)

As I've expressed in another comment, my preferred way of thinking about this symmetric pattern goes something like this:

The probability that (a prime congruent to x mod b is followed by a prime congruent to y mod b) seems to be equal to the probability that (a prime congruent to -y mod b is followed by a prime congruent to -x mod b).

I still haven't figured out whether it ought to be obvious, though if it is then I expect the language I've used above to be relevant. It's definitely not trivially obvious, because it's not an exact equality: if the pattern only shows up clearly after you've accumulated thousands or millions of primes, then it doesn't seem that it could be enforced by any sort of exact transformation. (For example, if the symmetric entries were somehow just counting the same pairs in two different ways, the numbers ought to be precisely equal rather than just increasingly close.)

There are a rather a lot of other patterns in the data; I expect that at least some of them must be accounted for in the original paper, but I haven't more than glanced through it yet.

> It appears to have the same transition probabilities as in base 5 (with the two center rows and columns swapped):

Well, yeah that's because 1,3,7,9 modulo 5 becomes 1,3,2,4 which is 1,2,3,4 with the centre two swapped.

I find this more interesting than the article and want to know more.

So a number ending in 'x' is as likely to be followed by 'y' as a number ending in 'n-x' is to be preceded by a number ending in 'n-y' (where n is the base).

Can't think of a trivial reason why that would be the case, something weird is happening.

Framing this differently, a prime congruent to x mod b is as likely to be followed by a prime congruent to y mod b as a prime congruent to -y mod b is to be followed by one congruent to -x mod b.