Let's simulate some data based on your fits:
shape.x <- 0.586219900
rate.x <- 0.012159695
shape.y <- 0.491757644
rate.y <- 0.006671180
set.seed(1) # for replicability
x <- rgamma(186,shape=shape.x,rate=rate.x)
y <- rgamma(50,shape=shape.y,rate=rate.y)
Now, if you are only interested in whether the means of the underlying distributions differ, the canonical approach is a standard t test on the data themselves, without fitting anything:
Welch Two Sample t-test
data: x and y
t = -2.8045, df = 60.679, p-value = 0.006761
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
mean of x mean of y
Note that this is (asymptotically) valid. t tests make pretty weak assumptions on the shape of the distributions of the underlying data. This is because by standard theorems, means of samples (which we are interested in here) are asymptotically normally distributed. You have sample sizes of 186 and 50. In such situations, I'll usually use a t test without a second thought.
That said, your estimated rate parameters indicate that you have some pretty degenerate gammas (if, as @Glen_b suggests, they are gamma at all):
foo.x <- seq(min(x),max(x),by=.1)
foo.y <- seq(min(y),max(y),by=.1)
(Compare what Wikipedia authors think is a "standard gamma".)
So the question pops up whether we have enough samples for asymptotics to kick in and make the t test valid, after all. One possibility would be to nonparametrically bootstrap the difference in means:
n.boot <- 10000
foo <- rbind(data.frame(X=x,group="A"),data.frame(X=y,group="B"))
set.seed(1) # for replicability
boot.nonparam <- boot(foo,
Now, I would say that $p=.0012$ from the nonparametric bootstrap is close enough to $p=.007$ from the t test, considering that the smaller group has 50 observations and no more.
Another alternative would be a parametric bootstrap, if you are really sure (as per theory) that your data actually should be gamma. Bootstrap your data, fit gammas to the resampled data, then either resample from the fitted gammas (my approach below), or directly calculate the expectations of the fitted gammas as shape/rate:
boot.param <- boot(foo,
gamma.A <- fitdistr(data[indices,"X"][data[indices,"group"]=="A"],"gamma")
resample.A <- rgamma(sum(data[indices,"group"]=="A"),
gamma.B <- fitdistr(data[indices,"X"][data[indices,"group"]=="B"],"gamma")
resample.B <- rgamma(sum(data[indices,"group"]=="B"),
Now we only have $p=.02$, quite a bit more than above.
Bottom line: you should really look at your data closely. If you fit a gamma, it will be pretty degenerate - so degenerate, in fact, that I'd be wary of making any inferences about whether the means of these fitted (!) gammas will be significantly different. I'd most trust the nonparametric bootstrap - but even that can be unstable for heavily skewed distributions (Good mentions a threshold of $n=100$ in his 2006 book). One other alternative would be a permutation test, where you assess the null distribution of means under random permutation of group labels.
And the bottom bottom line: of course p values are not "probabilities that populations (or means) are significantly different". See here.