When the constraints are not holonomic, why is it not possible to find such $q_i$s that $\delta q_i$s are independent of each other? In the book of Classical Mechanics by Goldstein, at page 20, it is given that

However, I cannoot understand from what has been presented so far that when the constraints are not holonomics, why is it not possible to find such $q_i$s that $\delta q_i$s are independent of each other ?
 A: OP is considering the very last step in the derivation of Lagrange equations (1.53) from d'Alembert's principle (1.45) in Goldstein. OP asks what would happen if the constraints are not holonomic?
Well, I disagree with Goldstein that thus far no restrictions has been on the nature of the constraints [...]. Frankly the holonomic conditions were used already in the very first step of the derivation, namely in the very definition of generalized  coordinates (1.38). So we would strictly speaking have to redo the whole derivation from scratch!
If we assume that Goldstein in the mentioned quote is not talking about the holonomic constraints used to define the generalized coordinates (1.38) in the first place, but that he is merely entertaining the possibility of some other (not necessarily holonomic) constraints (that don't perform virtual work), then he is basically saying that the presence of these extra constraints would imply that the virtual displacements $\delta q^j$ are not independent (because the system can not cover all generalized positions $q^j$ in the generalized configuration space), and that we therefore can not conclude Lagrange equations in the simplest form (1.53).
Now if these other constraints happen to also be holonomic, then we could pick a smaller set of independent generalized coordinates for the system.
Let us mention for completeness that a modified set of Lagrange equations can in principle be derived for so-called semi-holonomic constraints by introducing Lagrange multipliers $\lambda^i$. The resulting modified Lagrange equations are displayed in my Phys.SE answer here.
A: If a constraint on the motion of the system are holonomic, by definition of 'holonomic', the only thing it restricts is the position of the system to a subspace of lower dimension (subspace of the whole unrestricted space of the old coordinates). A set of new coordinates (which has fewer number of them than the old set) can be then introduced to describe position in this subset,  free of any remaining constraint. So, any continuous trajectory entirely in the subspace is possible. Think of how constraining a mass point in 3D space to a surface of a sphere allows us to introduce 2 new coordinates of the point, which we are then free to continuously change in any way (latitude and longitude).
However, if the constraints on old coordinates are not holonomic, the restriction on the possible motion is not simply a restriction of its position to a lower dimensional subset of the old space. The typical example is restricting things such as how fast and by which path the system may change its position. While the set of possible positions may be unrestricted, the possible trajectories connecting those positions in that set are restricted and so the differentials of any valid set of coordinates $q_i$ will be somehow restricted, either mutually or by some time-dependent condition. Think of a marble made of rubber, bound to always touch a desk, with no sliding possible. For any pre-selected contact point on the desk and any angular orientation of the marble, it can be put in there while obeying the non-sliding constraint at every step, just by rolling it cleverly. But it cannot be done "directly", by simply translating and rotating in the most obvious way, because that could require sliding. Instead, it has to be rolled back and forth appropriately to acquire the right angular orientation, while also moving the touch point to the destination. This means the possible trajectories connecting two possible positions are constrained, not every trajectory is possible. The constraint is that the speed of the lowest point always has to be zero, which translates to a constraint on the simultaneous values of variation of position on the plane and variation of the angular variables. 
If there was a complete set of independent coordinates for a system, there wouldn't be any restriction of the above kind on the trajectories, because, well, the coordinates would be completely independent of anything. Such system thus could not be having nonholonomic constraints.
