Strange effective degrees of freedom (smoothness) selected for smooth component in GAM model with mgcv Consider the following very simple example:
set.seed( 1 )
SimData <- data.frame( X = runif( 1000, 0, 1 ) )
SimData$Y <- rnbinom( nrow( SimData ), mu = 100*sin( SimData$X*2*pi )+100, size = 10 )
plot( gam( Y ~ s( X ), data = SimData, family = nb( link = log ) ) )

Here is the result:

The effective degrees of freedom is very-very close to 9, giving rise to the suspicion that the default basis dimension is not large enough. Let's increase it:
plot( gam( Y ~ s( X, k = 20 ), data = SimData, family = nb( link = log ) ) )


Hm. The overall picture is essentially the same, yet, the EDF suggests that the basis dimension is still not large enough!
Even if we increase k to 30, the EDF still gets larger (24.8), so, essentially EDF seems to simply follow the k limit, which is pretty bizarre... (especially for such a simple functional form, and especially that is was already well captured by the default model).
EDIT (07 Sep, 2017): According to an answer to the original question, the application of adaptive splines (bs="ad") can be the solution to this problem. However...
Let's take another simple example:
set.seed( 1 )
SimData <- data.frame( X = runif( 1000, 0, 1 ) )
SimData$Y <- rnbinom( nrow( SimData ), mu = 100*sin( SimData$X*2*pi*3 )+1000, size = 10 )
plot( gam( Y ~ s( X ), data = SimData, family = nb( link = log ), method = "REML" ) )


Seems perfect! Let's now "spoil" it:
SimData$Y[ SimData$X<=0.05|SimData$X>=0.95 ] <- 0

This gives rise to the original problem: EDF is 8.92 with default k, 18.91 if k=20, 46.74 if k=50 etc. As an illustration, here is the k=50 case:

So indeed, we have the original problem. Let's try therefore bs="ad":
plot( gam( Y ~ s( X, bs = "ad" ), data = SimData, family = nb( link = log ), method = "REML" ) )


So unfortunately the problem remained, even with bs="ad"! (It's really the same situation: increasing k to 50 gives an EDF of 44.72, k=100 gives 78.03. That's why I decided to edit this question instead of starting a new one: this seems to be somehow the same story...)
 A: Fitting on the log scale is causing the problem; on that scale, as shown in your plot, the degree of smoothness is not really the same across the function; at low values of x there is little curvature in the true function, but there is a lot of curvature around the minimum at ~ x = 7.5. As, by default, the spline wiggliness is controlled by a single smoothness parameter, this has the effect of assuming the same degree of smoothness or curvature everywhere. As you increase k, the model is able to adapt to the data around the rapid rise and fall, but in doing so it also adapts, because it has too given the smoothness parameter/penalty, to sampling noise elsewhere in the data. 
With these data I was able to achieve an acceptable fit by changing the link function to the sqrt link via family = nb(link = "sqrt"). On this scale the rise and fall is not as sharp or constrained to a narrow range of x as it is on the log scale.
m <-  gam(Y ~ s(X, k = 20), data = SimData, family = nb(link = "sqrt"), method = "REML")

The fitted spline is:

There is still some evidence of adaptation to sampling noise in the region where x is low, but this is quite trivial and by any reasonable metric would be indistinguishable from the true function I would guess.
The basis dimension test implemented in gam.check() also seems reasonably happy with this — note that I increased k to 20 but this appears sufficient for this model, whereas k = 10 wasn't:
> gam.check(m)

Method: REML   Optimizer: outer newton
full convergence after 4 iterations.
Gradient range [2.520842e-06,0.0003287577]
(score 4439.901 & scale 1).
Hessian positive definite, eigenvalue range [5.53589,345.0859].
Model rank =  20 / 20 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

       k'  edf k-index p-value
s(X) 19.0 14.9    1.01    0.76

The diagnostic plots are also pretty reasonable:

If you want to fit on the log scale, an adaptive spline can be used, which also appears to give good results for these data. An adaptive spline uses multiple smoothness penalties (parameters) for each spline and in effect is like fitting several separate splines to parts of data, each with their own smoothness parameter. Such splines require a lot more data to achieve good results however.
mad <-  gam(Y ~ s(X, k = 20, bs = "ad"), data = SimData, family = nb, method = "REML")

The fitted spline is now much smoother for low values of x, but there is perhaps a little overfitting to sampling noise at large x. On the scale of the response this is going to be trivial however.

and gam.check() is also happy that k = 20 is sufficient for this basis/model/data:
> gam.check(mad)

Method: REML   Optimizer: outer newton
full convergence after 11 iterations.
Gradient range [-0.0002868044,0.001512152]
(score 4461.377 & scale 1).
Hessian positive definite, eigenvalue range [0.0005701379,344.5224].
Model rank =  20 / 20 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

       k'  edf k-index p-value
s(X) 19.0 12.9    0.97    0.24

and again, the diagnostics don't look that bad either

although there is a hint of non-constant variance in residuals in the plot to the top right.
I'll also add that the test in gam.check() for the dimension of the basis is not infallible. It is a guide at best. If the model looks good in other respects, I would place less emphasis on the dimension test results. However, for the default k = 10 model you fitted, the model isn't really sufficiently wiggly to capture the low counts around the dip and it is quite badly biased there.

I'd argue that whilst just increasing k wasn't helping here, the initial model was also not adequate for these data.
