Multiply integrand with -1, and the precision changes?

Short answer:

One workaround is to use Method -> "LevinRule" instead.

---

Long answer:

As mentioned by Xavier in a comment above, changing BesselJ[0, x] to BesselJ[0, Re[x]] resolves the issue:

NIntegrate[{1, -1} BesselJ[0, Re@x], {x, 0, ∞}, WorkingPrecision -> 32, 
 Method -> "ExtrapolatingOscillatory"]
    Precision /@ %

{0.99999999999999999999999999999979, -0.99999999999999999999999999999979}
{32., 32.}

But why it works? Codes in the UPDATE of this answer solve the mystery. The truth is, NIntegrate has internally switched to "LevinRule":

Pictured by Simon Wood's shadow.

The output is indeed the same as that with option Method -> "LevinRule:

NIntegrate[{1, -1} BesselJ[0, x], {x, 0, ∞}, WorkingPrecision -> 32, 
 Method -> "LevinRule"]
Precision /@ %

{0.99999999999999999999999999999979, -0.99999999999999999999999999999979}
{32., 32.}

Notice in this case you'll be unable to adjust the MaxRecursion option because of a bug mentioned here, but this seems not to be a big deal:

pmhankelLevin[p_, sign_: 1, prec_: 32] := 
 NIntegrate[sign ξ BesselJ[0, ξ] f[p, ξ], {ξ, 0, ∞}, 
  WorkingPrecision -> prec, Method -> "LevinRule"]

lst1 = pmhankel[#, -1] & /@ Range@32;
(* You'll see warning NIntegrate::ncvb and be unable to eliminate it 
   by adjusting MaxRecursion option because of the bug mentioned above 
   when generating lst2 and lst3 *)
lst2 = pmhankelLevin /@ Range@32;
lst3 = pmhankelLevin[#, -1] & /@ Range@32;
(* But the difference between lst2, lst3 and lst1 is negligible: *)
Max /@ {(lst1 + lst2)/lst1, (lst1 - lst3)/lst1}

{2.27320119973*10^-6, 2.27320119973*10^-6}

Though a workaround is found, the reason why the precision is influenced by the sign when "ExtrapolatingOscillatory" method is used still remains unclear. I'm looking forward to answer(s) addressing the issue.

The oddity in this case comes from NSum which is being called in a certain way from NIntegrate. This is a simple example that has roughly the same behavior (note in this case the exact result is known to be $\mp \ln 2$):

NSum[(-1)^n/n, {n, 1, Infinity}, 
          Method -> {"AlternatingSigns", Method -> "WynnEpsilon"}, WorkingPrecision -> 32]

(* -0.6931471805599453094172318803247 *)

NSum[-(-1)^n/n, {n, 1, Infinity}, 
          Method -> {"AlternatingSigns", Method -> "WynnEpsilon"}, WorkingPrecision -> 32]

(* 0.693147180559945309417232 *)

where the second result has several digits fewer than the first.

Is that a bug? Not necessarily, because both results have at least 16 correct digits which certainly attains the default PrecisionGoal, which is WorkingPrecision/2.

Still, I agree the consistency could be improved in this case and I have filed a report for the developers to take a look.

---

Update

This has been improved in the just released Mathematica 11.0.