Inconsistent results with gammainc for values -0.6 < a < -0.2 and x < 0.01

I had been getting some unexpected results with the gammainc() function and I have isolated an issue within the function. The following code compares results from expint::gamminc(a,x) to pracma::incgam(x,a) using -0.6 < a < -0.2 and x < 0.1.

library(ggplot2)

# create values for x and a for sensitivity analysis
x <- 10^(-5:-1)
a <- seq(-0.6, -0.2, 0.01)

# get all combinations of a and x
data <- expand.grid(a = a, x = x)

# calculate results from both functions
data$exGamma <- mapply(expint::gammainc, a=data$a, x=data$x)
data$pmGamma <- mapply(pracma::incgam, x=data$x, a=data$a)

# calculate the difference between the functions.
data$funCompare <- data$exGamma - data$pmGamma
data$x <- factor(data$x)

ggplot(data) + geom_line(aes(a, funCompare, color = x)) +
  ylab("expint::gammainc(a,x) - pracma::incgam(x,a)")

As you can see, from -0.6 < a <= -0.5, the functions return the same value. For values x < 0.01, an abrupt discontinuity appears in one of the functions as soon as a > -0.5. The discontinuity gets larger as x approaches zero. It dissipates as a approaches -0.3.

The question remains, which function has the discontinuity? The graph below shows that the issue is with expint::gammainc().

a <- seq(-0.501, -0.499, length.out = 1000)
exGamma <- sapply(a, expint::gammainc, x = 1e-05)
pmGamma <- sapply(a, pracma::incgam, x = 1e-05)
plot(a, exGamma, ty = "l", ylab = "gamma")
lines(a, pmGamma, lty = 2, lwd =2)
legend(-0.4999, 635, c("expint::gammainc(a,1e-05)", "pracma::incgam(1e-05, a)"), lty = c(1,2))