The combined data object has only the 2000 or 3000 most variable features, which were chosen at the beginning of the workflow. All the other features are lost in the integrated data. Is there any option to integrate the data having all features in the integrated data object?
I think you are missing the concept here. As igor says the integration runs on the reduced dimensions, most commonly PCA, not on a per-gene basis. And for this you only want to include the variable genes as all other genes do not any information. The integrated values are only useful for clustering and visualization which, again, is commonly performed on a reduced dimension level. For something like differential gene expression one would not use the integrated values as the integration precedure creates dependencies between the data points and might change magnitude and directions. There are many threads on this already, e.g. issues at Seurat Github or at section 13.8 in OSCA:
At this point, it is also tempting to use the corrected expression values for gene-based analyses like DE-based marker gene detection. This is not generally recommended as an arbitrary correction algorithm is not obliged to preserve the magnitude (or even direction) of differences in per-gene expression when attempting to align multiple batches. For example, cosine normalization in fastMNN() shrinks the magnitude of the expression values so that the computed log-fold changes have no obvious interpretation. Of greater concern is the possibility that the correction introduces artificial agreement across batches.
In summary, for integration you should use the most variable genes which is the reason why this is the default in most workflows.